点击蓝字↑↑↑“微生态”，轻松关注不迷路

本文由阿童木根据实践经验而整理，希望对大家有帮助。

原创微文，欢迎转发转载。

只做宏基因组太单调？为什么不试试宏基因组Binning呢？一次测序，“宏基因组”+“Binning”两种分析，微生太帮您一站式处理宏基因组难题。现在，微生太免费向所有人分享Binning的整套分析流程，包含：生信分析代码和R语言绘图代码。我们一共设计了7个课时，每周一次，课表(进度)如下。

对Binning分析、R语言绘图感兴趣的朋友千万别错过。错过也没关系，每次课程不仅有回放，还有技术贴带您回顾课程内容。

图1

下面我们一起来回顾第四节课的主体内容吧。

 一、分析内容

COG，即Clusters of OrthologousGroups of proteins(同源蛋白簇)。COG是由NCBI创建并维护的蛋白数据库，根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成。COG分为两类，一类是原核生物的(一般称COG)，另一类是真核生物(一般称KOG)。通过比对可以将某个蛋白序列注释到某一个COG中，每一簇COG由直系同源序列构成，从而可以推测该序列的功能。ENZYME收录了7大类酶的四级分类信息。EC编号或EC号是酶学委员会(EnzymeCommission)为酶所制作的一套编号分类法，每一个酶的编号都以字母“EC”起头，接着以四个号码来表示，这些号码代表逐步更细致的为酶作出分类。Ribosomal RNA genes (rRNA)、Transfer RNA genes (tRNA)、Non-coding RNA(ncRNA)分别使用RNAmmer、Aragorn、Infernal进行预测。

COG官网：https://www.ncbi.nlm.nih.gov/COG/

COG注释信息：

http://eggnogdb.embl.de/download/eggnog_4.5/data/NOG/NOG.annotations.tsv.gz

EC官网：https://www.qmul.ac.uk/sbcs/iubmb/enzyme/

 二、分析方法

使用Prokka软件对每个Bin进行功能注释，获得每个Bin中的rRNA、tRNA、tmRNA、基因、直系同源蛋白簇(COG)、EC注释信息。使用Kofamscan进行KEGG功能注释。使用eggnog-mapper进行GO注释。使用Diamond软件和CAZyme数据库对每个Bin进行碳水化合物酶注释，获取每个Bin中的碳水化合物酶信息。

三、分析过程

1 输入数据

输入数据是经质检筛选后的clean bins，存放于Bin_pick文件夹，如下：

图2

2 分析流程

# 1. prokka基因组注释for i in Bin_all/Bin_pick/*.fa; do    file=${i##*/}    base=${file%.fa}    prokka $i --outdir Bin_all/Bin_prokka/$base --prefix $base --metagenome --cpus $thread --kingdom Bacteria    echo -e "\033[32m$i prokka Done...\033[0m"done

# 2. 统计prokka注释tsv文件R CMD BATCH --args /home/cheng/huty/softwares/script_bin/prokka_fun_count.R

# 3. 抽提COG信息for i in Bin_all/Bin_prokka/prokka_out_table/bin.*.tsv; do    file=${i##*/}    fold=${file%.tsv}    mkdir Bin_all/cog/$fold    cat $i | awk -F"\t" '{if($6 != "") printf("%s\t%s\n", $1, $6)}' > Bin_all/cog/$fold/${fold}_cog.txt    echo -e "\033[32m\t $i cog extract Done...\033[0m"done# 4. COG分类注释level 1 level 2 (func)for i in Bin_all/cog/bin.*; do    fold=${i##*/}    Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_annotation.R $i/${fold}_cog.txt /home/cheng/huty/databases/cog/cog_anno.txt ${fold}_cog_annotation.txt    echo -e "\033[32m\t $i annotate cog Done...\033[0m"done# 5. 统计func数量for i in Bin_all/cog/bin.*; do    fold=${i##*/}    cat $i/${fold}_cog_annotation.txt | awk -F"\t" 'BEGIN{OFS="\t"}{if($6 != "NA") print $3}' | sed '1d' | sort | uniq -c | awk 'BEGIN{OFS="\t"}{print $2,$1}' | sed '1 ifunc\tCount' > $i/${fold}_cog_sum.txt    echo -e "\033[32m\t $i summary annotated cog Done...\033[0m"done

# 6. 统计结果注释level 1for i in Bin_all/cog/bin.*; do    fold=${i##*/}    Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_count_annotation.R $i/${fold}_cog_sum.txt /home/cheng/huty/databases/cog/cog_category.txt ${fold}_cog_sum_annotation.txt    echo -e "\033[32m\t $i 分类 annotated cog Done...\033[0m"done# 7. 可视化统计注释结果for i in Bin_all/cog/bin.*; do    fold=${i##*/}    Rscript /home/cheng/huty/softwares/script_genome/prokka_cog_count_annotation_bar.R $i/${fold}_cog_sum_annotation.txt ${fold}_cog_sum_annotation.pdf    convert -density 400 -quality 200 $i/${fold}_cog_sum_annotation.pdf $i/${fold}_cog_sum_annotation.png    echo -e "\033[32m\t $i 分类可视化 annotated cog Done...\033[0m"done

3 统计：prokka功能注释

从prokka注释结果问价夹中提取bin.tsv文件到一个文件夹，打开其中一个bin.tsv查看文件格式，文件格式(左到右列依次是)：预测基因编号、functiontype、基因长度、基因名称、EC编号、COG编号、蛋白名称。

图3

用自写的R脚本prokka_fun_count.R统计汇总所有BIN的prokka注释结果，格式如下(左到右列依次为)：BIN编号、注释总数、EC数、COG总数、GENE总数、CDS总数、各类RNA总数。

图4

4 统计：COG分类注释

提取bin.tsv文件中的第一列和COG列，剔除未得到COG注释的行，接着用R语言merge函数(prokka_cog_annotation.R)和自制COG注释文件(cog_anno.txt)进行COG二级注释，获得function、level_1、level_2注释信息，然后统计function数量获得丰度表，接着用类似的R merge方法(prokka_cog_count_annotation.R)做注释(cog_category.txt)。

过程如下：

图5

结果如下：

图6

5 COG分类注释可视化

R脚本(prokka_cog_count_annotation_bar.R)代码：

#!/usr/bin/env Rscriptargs = commandArgs(T)route_file = unlist(strsplit(args[1], "/"))route = paste(route_file[1:(length(route_file)-1)], collapse="/")setwd(route)file_name = route_file[length(route_file)]library(ggplot2)data = read.table(file_name, header=T, sep="\t")data_sort = data[order(data[,3], data[,2], decreasing=F),]result = ggplot(data_sort, aes(x = func, y = Count, fill=level_1)) +    geom_bar(stat = "identity") +    coord_flip() +    #theme(axis.text.x=element_text(angle=45, hjust=1)) +    #angle：调整横轴标签倾斜角度    labs(x = "COG level 2", y = "Count", fill = "COG level 1") +    scale_x_discrete(limits=factor(data_sort[,1])) +    #scale_y_continuous(expand=c(0, 0)) +    theme(panel.grid=element_blank(), panel.background=element_rect(color='black', fill='transparent'))ggsave(result, filename=args[2], height = 7, width = 7)

对示例数据(联系我们获取)进行绘图，得到结果如下图，横轴是function计数，纵轴function代号(与level_2一一对应)，柱子的长度对应function计数，用level_1注释信息作为legend对应颜色。

图7

你可能还喜欢

1 技术贴 | 16S专题 | 简单介绍如何用自己的笔记本处理高通量16S数据

2 技术贴 | 宏基因组专题 | 组装工具盘点和比较

3 技术贴 | R语言菌群Alpha多样性分析和绘图

4 技术贴 | 宏转录组专题 | DDBJ数据库：宏转录组测序数据下载

5 技术贴 | R语言pheatmap聚类分析和热图

微生态科研学术群期待与您交流更多微生态科研问题

(联系微生态老师即可申请入群)

了解更多菌群知识，请关注“微生态”。