导读
prokka注释细菌基因组,获得该基因组的CDS、Gene、RNA、EC、COG等注释信息。
一、基本信息
文章:Prokka: rapid prokaryotic genome annotation
杂志:bioinformatics
时间:2014
github:https://github.com/tseemann/prokka
二、prokka注释bin
####################################
##### prokka:type gene EC COG #####
####################################
echo -e "\033[32mprokka基因组注释: \033[0m"
for i in Bin_all/Bin_pick/*.fa; do
file=${i##*/}
base=${file%.fa}
prokka $i --outdir Bin_all/Bin_prokka/$base --prefix $base --locustag $base --metagenome --cpus $thread --kingdom Bacteria
echo -e "\033[32m$i prokka Done...\033[0m"
done
prokka参数:
--outdir 输出文件夹,自动创建
--prefix 输出文件名前缀
--locustag 预测基因编号前缀
--metagenome 高度片段化基因组,宏基因组
--kingdom Archaea|Bacteria|Mitochondria|Viruses 注释模式
三、结果汇总
echo -e "\033[32m集合prokka注释tsv文件:\033[0m"
for i in Bin_all/Bin_prokka/bin*; do
fold=${i##*/}
mv $i/$fold.tsv Bin_all/Bin_prokka/prokka_out_table # prokka注释结果
done
echo -e "\033[32m统计prokka注释tsv文件:\033[0m"
R CMD BATCH --args /home/cheng/huty/softwares/script_bin/prokka_fun_count.R
- 拙劣的R脚本 (′д` )…彡…彡
setwd("Bin_all/Bin_prokka/prokka_out_table")
files=list.files(pattern="bin.*.tsv") # 读取所有文件名
Bin_ID=vector()
for(i in 1:length(files))
{
Bin_ID[i]=as.character(strsplit(files[i], split=".tsv"))
# 提取所有文件名
}
ml=list()
for(i in 1:length(files))
{
ml[[i]]=read.table(files[i], sep='\t', na.string="", stringsAsFactors=F, header=T, quote="", comment.char="")
# 读取所有数据框
}
EC_num=vector()
for(i in 1:length(files))
{
EC_num[i]=lengths(ml[[i]]["EC_number"])-sum(is.na(ml[[i]]["EC_number"]))
# EC总数/BIN
}
COG_num=vector()
for(i in 1:length(files))
{
COG_num[i]=lengths(ml[[i]]["COG"])-sum(is.na(ml[[i]]["COG"]))
# COG总数/BIN
}
GENE_num=vector()
for(i in 1:length(files))
{
GENE_num[i]=lengths(ml[[i]]["gene"])-sum(is.na(ml[[i]]["gene"]))
# 基因总数/BIN
}
Total_num=vector()
for(i in 1:length(files))
{
Total_num[i]=lengths(ml[[i]]["ftype"])
# 功能预测总数/BIN
}
CDS_num=vector()
rRNA_num=vector()
tRNA_num=vector()
tmRNA_num=vector()
for(i in 1:length(files))
{
CDS_num[i]=0
rRNA_num[i]=0
tRNA_num[i]=0
tmRNA_num[i]=0
for(j in 1:length(ml[[i]][,"ftype"]))
{
if(ml[[i]][j, "ftype"]=="CDS")
{
CDS_num[i]=CDS_num[i]+1
}
else if(ml[[i]][j, "ftype"]=="rRNA")
{
rRNA_num[i]=rRNA_num[i]+1
}
else if(ml[[i]][j, "ftype"]=="tRNA")
{
tRNA_num[i]=tRNA_num[i]+1
}
else if(ml[[i]][j, "ftype"]=="tmRNA")
{
tmRNA_num[i]=tmRNA_num[i]+1
}
}
}
prokka_result=data.frame(Bin_ID, Total_num, EC_num, COG_num, GENE_num, CDS_num, tRNA_num, rRNA_num, tmRNA_num)
write.table(prokka_result, file="../prokka_out_table.txt", sep="\t", quote=F, row.names=F)
四、结果
Bin_ID Total_num EC_num COG_num GENE_num CDS_num tRNA_num rRNA_num tmRNA_num
bin.102 2816 1062 1151 1459 2786 27 2 1
bin.104 3184 1173 1580 1892 3130 51 2 1
bin.105 2625 974 1166 1422 2585 37 2 1
bin.106 3610 1196 1281 1674 3564 42 3 1
bin.107 1414 459 466 595 1389 22 3 0
bin.108 1695 646 942 1081 1661 31 2 1
bin.109 1568 637 970 1080 1543 22 2 1
bin.11 1205 483 623 740 1176 28 0 1
更多:
宏基因组功能注释(以COG为例
每日一生信--COG注释(终结版)
数据库地址:https://ftp.ncbi.nih.gov/pub/COG/
网友评论