宏基因组分箱（六）Prokka基因预测，注释COG、EC、RNA

导读

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

$\color{green}{😀😀KOOK😀😀}$

最后编辑于：2021.05.06 10:38:45