{示例企业}分箱数据分析结题报告

温馨提示:请使用火狐或者Chrome的网页浏览器来查看报告

{示例企业}宏基因组分箱数据分析结题报告


一 概述

宏基因组分箱(Binning)是将宏基因组测序得到的混合了 不同微生物的序列reads或序列组装得到的contigs或scaffolds按物种分开归类的过程。 这些分开归类的序列被称为宏基因组组装基因组(metagenome-assembled genomes,MAGs)。 传统的单物种全基因组序列都是经纯培养之后,再进行全基因组de novo测序才获得的, 但是环境中存在着大量的不可培养微生物(uncultured candidate bacterial species,未培养候选菌种), 而宏基因组分箱及相关技术不仅有助于获得不可培养微生物的基因组序列,还有以下诸多功能: (1)发现新物种,预测新物种基因,利用现有数据库分析新物种的功能; (2)扩充微生物基因组数据库,增加微生物多样性,提高宏基因组数据reads mapping率(检出率); (3)有助于宏基因组技术开发; (4)助力“感兴趣”微生物群结构和功能的研究; (5)为菌群的分类和功能描述提供了更多的解决方案。 早在2011年,science上的一篇文章就用了宏基因组Binning技术对来自牛瘤胃的样本进行了宏基因组测序研究。 该研究从268 Gbp的宏基因数据中成功Binning出了15个不能培养的微生物的全基因组序列 [17]

从那以后,宏基因组Binning技术开始被更多的人关注和重视, 也逐渐出现了很多宏基因组Binning相关的工具。 Metabat是近几年所有分箱工具中最受欢迎的工具(截止2020年, 引用达460+)。 2019年发表在PeerJ上的新版Metabat2更是在完成度、效率等多方面均优于Metabat和同类工具 。 仅在2019年当年Metabat2就已经被Cell、Nature Biotechnology、Genome Biology等多篇高水平期刊引用 。 因此{示例企业}从众多分箱工具中选择Metabat2进行分箱。 不仅如此,{示例企业}宏基因组分箱分析流程中使用的每一个软件都是根据效率、准确度等多个参数从众多同类软件中精选出来的, 例如序列组装工具Megahit,MAG质检工具Checkm,功能预测工具Prokka,蛋白比对工具Diamond等等。


二 项目流程

2.1 实验流程

图2-1 测序流程图

2.2 数据分析流程

本报告中,宏基因组分箱的数据分析流程如图2-2标准分析部分所示, 主要步骤为质控,去宿主,分箱,MAG去冗余和筛选,MAG丰度计算(定量), MAG功能注释,MAG物种注释,MAG圈图等。

图2-2 宏基因组分箱数据分析流程


三 结果文件夹解读

00-QCStats序列质控和去宿主序列

    |--00-QCStats
    | |--1-QC_report_Rawfastq/*.html [原始数据fastqc质检结果]
    | |--2-QC_report_Filtered/*.html [序列质控和去宿主序列之后的fastqc质检结果]
    | |--reads_summary.txt [数据产出质量情况一览表]

本项目采用Illumina测序平台对测序样本进行双端测序。 基于FASTQ格式的测序文件是一种存储序列信息的特定文件,推荐用Notepad++等文本编辑器或者在电脑终端中打开。 FASTQ文件每四行对应一条测序Read:第一行以符号“@”起始,对应于序列ID和相应的描述信息; 第二行为实际测得的碱基序列;第三行以符号“+”起始; 而第四行的字符串则记录了第二行序列中每个碱基所对应的测序质量。

采用 Illumina测序平台测序获得的原始数据(Raw Data)存在一定比例低质量数据, 为了保证后续分析的结果准确可靠,需要对原始的测序数据进行预处理,包括质控 (Trimmomatic[16] 参数:ILLUMINACLIP:adapters_path:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50), 和去宿主序列(Bowtie2[2]参数:--very-sensitive),获取用于后续分析的有效序列(clean data)。 测序数据预处理统计结果见表 3-1。序列质控步骤关键参数解释如下:

  1. 去除接头序列 (参数ILLUMINACLIP:adapters_path:2:30:10);
  2. 扫描序列(4bp滑窗大小),如果平均质量分低于20(99%正确率),切除后续序列(参数SLIDINGWINDOW:4:20);
  3. 去除最终长度小于50bp的序列(参数MINLEN:50)。
表3-1 数据产出质量情况一览表

Sample ID InsertSize(bp) SeqStrategy RawReads(#) Raw Base(GB) %GC Raw Q20(%) Raw Q30(%) Clean Reads(#) Cleaned(%) Clean Q20(%) Clean Q30(%)
C1 350 (150:150) 23613034 7.08 61 97.70 93.78 22416119 94.93 98.80 95.62
C2 350 (150:150) 27584954 8.28 61 97.10 92.62 26060963 94.48 98.55 94.94
C3 350 (150:150) 23777022 7.13 61 98.45 95.71 22889304 96.27 99.27 97.08
C4 350 (150:150) 27007788 8.10 61 97.63 93.65 25636029 94.92 98.77 95.52
T11 350 (150:150) 29809022 8.94 62 97.62 93.48 28346006 95.09 98.63 95.16
T12 350 (150:150) 24702558 7.41 61 97.72 93.78 23570823 95.42 98.72 95.43
T13 350 (150:150) 28274974 8.48 61 97.63 93.54 26703496 94.44 98.69 95.32
T21 350 (150:150) 25016746 7.51 60 97.74 93.74 23831131 95.26 98.70 95.35
T22 350 (150:150) 24980036 7.49 60 97.55 93.21 23555813 94.30 98.57 94.96
T23 350 (150:150) 24026221 7.21 62 97.64 93.58 22758509 94.72 98.69 95.34

  1. Sample ID: 样品名;
  2. InsertSize(bp): InsertSize是在建库切胶时选择的长度,合适的InsertSize能避免测序接头污染;
  3. SeqStrategy: 测序策略(一般为双端,各150bp);
  4. RawReads(#):测序Raw reads的数量;
  5. Raw Base(GB): 以GB为单位的Raw reads数量, 测序原始数据的总碱基数,即为Raw reads数量乘以测序长度算得;
  6. %GC:G/C碱基数占总碱基数量的百分比;
  7. Clean Reads(#):过滤(质控和去宿主序列)后获得的Clean reads 数量;
  8. Cleaned(%):过滤后剩余的序列数占Raw reads 的百分比;
  9. Q20: 质量分高于20的碱基所占比例;
  10. Q30,质量分高于30的碱基所占比例。

01-MAG分箱与MAGs质控

文件说明:
|--01-MAG
| |--1-MAG_raw/*.fa [所有样本分箱得到的MAGs]
| |--2-MAG_pick/*.fa [筛选和去冗余后的MAGs]
| |--MAGs_picked_summary.txt [统计所有筛选得到的MAGs:完成度、污染度、N50、N90、Size、GC含量]
| |--checkm_raw.txt [所有MAGs的污染度、完整度]

对于每个样本,运用Megahit[1]软件, 将样本去宿主基因后的clean reads进行组装(megahit默认组装参数), 得到contigs;运用Bowtie2[2] 进行建索引和序列比对获得clean reads在contigs中的比对信息; 用Samtools软件[3]对比对结果进行格式转换和排序; 用Metabat2[4] 自带的jgi_summarize_bam_contig_depths程序计算contigs深度; 挑选长度大于1500bp的contigs(Metabat2参数“-m 1500”, 分箱的contig最短1500bp, 软件算法要求, 无法再降), 结合上一步得到的contigs深度,用Metabat2对contigs进行分箱; 用RefineM [5] 对得到的MAGs进行提纯,去除MAGs中污染度较高的contigs。

合并所有样本中得到的MAGs,用CheckM[6] 中的lineage_wf流程评估MAGs的完成度和污染度, 并用dRep[7] 软件对MAGs进行去冗余(dRep 参数: --P_ani 0.9 --S_ani 0.95 --completeness 50 --contamination 50) (完成度、污染度阈值, 可在创建流程前联系技术支持修改, 创建流程后, 无法修改)。

表3-2 筛选得到的MAGs统计

MAGID Completeness Contamination N50 N90 Size GC_percent
MAG.C3.5 91.61 3.46 16222.0 5435.0 2330867.0 0.487967
MAG.T11.11 96.46 1.01 136750.0 61426.0 3819437.0 0.499357
MAG.T11.1 95.14 3.60 23709.0 8074.0 3869115.0 0.681603
MAG.T11.53 91.38 2.59 17579.0 5549.0 2851598.0 0.383465
MAG.T12.14 90.04 1.74 13545.0 4350.0 3882823.0 0.624192
MAG.T12.17 95.81 0.00 97302.0 20803.0 6036892.0 0.515289
MAG.T12.29 95.60 1.44 25269.0 8491.0 2707859.0 0.696592
MAG.T12.45 90.17 4.96 12372.0 4035.0 3176878.0 0.728136
MAG.T13.10 91.72 1.01 20692.0 6182.0 2282946.0 0.446710
MAG.T13.13 93.52 9.08 12099.0 3571.0 3550691.0 0.690979
MAG.T13.17 99.51 0.65 103918.0 25255.0 3253597.0 0.372338
MAG.T13.41 99.71 1.09 55759.0 17241.0 4043162.0 0.421992
MAG.T13.54 94.17 2.22 17187.0 4623.0 3401258.0 0.489488
MAG.T13.5 97.68 3.74 35657.0 10168.0 6222019.0 0.611574
MAG.T13.7 91.54 5.49 9910.0 3889.0 3087226.0 0.429199
MAG.T21.37 94.21 4.55 20916.0 6551.0 3683299.0 0.645074
MAG.T21.41 92.67 4.09 17877.0 4667.0 4582529.0 0.594384
MAG.T21.45 99.49 1.52 105090.0 32404.0 4627794.0 0.369790
MAG.T21.66 91.99 2.91 48087.0 15377.0 4138707.0 0.671726
MAG.T22.14 91.62 1.17 42437.0 14457.0 1289198.0 0.528735
MAG.T22.15 96.71 1.93 30859.0 9781.0 3398348.0 0.456205
MAG.T22.30 91.29 4.86 15261.0 5611.0 3017272.0 0.616883
MAG.T22.39 93.30 2.84 47933.0 14119.0 2380561.0 0.594701
MAG.T22.3 100.00 0.19 103466.0 28233.0 4931295.0 0.463812
MAG.T22.40 96.02 2.76 29263.0 9011.0 2817253.0 0.555783
MAG.T22.42 96.97 1.52 146404.0 36830.0 1591561.0 0.628842
MAG.T22.46 91.14 6.08 19126.0 5693.0 4338382.0 0.651421
MAG.T23.25 94.02 4.55 13586.0 4862.0 4184703.0 0.697441
MAG.T23.3 96.02 1.14 97496.0 27779.0 3556789.0 0.642735
MAG.T23.45 95.53 1.59 80961.0 24507.0 2302337.0 0.634462

  1. 第一列:MAG编号;
  2. 第二列:MAG完成度;
  3. 第三列:MAG污染度;
  4. 第四列:Contig N50是指将一个MAG的Contig按照从大到小的顺序依次相加, 当相加的长度达到Contig总长度的一半时,最后一个加上的Contig长度;
  5. 第五列:Contig N90是指将一个MAG的Contig按照从大到小的顺序依次相加, 当相加的长度达到Contig总长度的90%时,最后一个加上的Contig长度;
  6. 第六列:一个MAG的所有Contig长度之和,即该基因组草图的总长度;
  7. 第七列:GC含量是指一个MAG中GC碱基占总碱基的比例。

02-MAG_Plot 分箱效果散点图

文件说明:
|--02-MAG_Plot
| |--MAGs_contig_summary.txt [散点图输入文件,含每个MAG的contig的GC含量和contig的深度信息]
| |--MAGs_contig_summary.svg [绘图结果,利用GC和depth数据把每个MAG的contig 画到图中]
		

利用可视化的方法直观展示MAG的由来和所含信息。 宏基因组分箱的原理是根据序列(contig/scaffold) 四核苷酸频率和序列丰度变化模式将序列分成一个个MAGs。 每个MAG的每个contig都会在此处理(metabat2分箱)的过程中有一个depth深度数据。 另外我们利用自己的python脚本计算每个MAG的contig GC含量。 有了contig GC含量和depth数据即可进行MAG可视化,绘制每个MAG中每个contig的散点图(图3-1)。

图3-1 Binned contigs可视化

横坐标是contig的GC含量;纵坐标是contig depth; 一个点代表一个contig,相同颜色的contig来自同一个MAG。

03-MAG_Abundance 丰度分析

03-MAG_Abundance
├── 1-Barplots [丰度柱形图]
├── 2-Heatmaps [丰度热图]
├── 3-Circos   [丰度圈图]
├── 4-SignificanceAnalysis  [丰度组间显著性差异比较]
├── 5-CorrelationAnalysis   [丰度相关性分析]
└── MAGs_abundance_table.xls [MAGs在各样本的丰度汇总表]
		

使用Metawrap[8] 的quant_bins模块(Salmon算法)计算每个MAG的丰度。 原理是将样本clean reads比对到MAG中的contigs, 计算单个样本某个MAG中所有contigs的总TPM。

我们用R语言pheatmap绘图函数绘制样品-MAG丰度热图。 热图可以以色块颜色深浅的方式表达MAG丰度的大小,还可以进行MAG-MAG聚类和样品-样品聚类,相似的丰度模式会被聚到一起。 然后,统计每个MAG的丰度总数,用ggplot geom_bar绘制柱形图并排序,展示整批数据中所有MAG的丰度情况。 如果客户提供的样品分组>=2且每组样品在三个以上,那么还可以用lefse进行组间差异分析寻找与分组有关的MAG。

1-Barplots文件夹中包含按照分组顺序排列的每个样本的MAG丰度柱形图(图3-2),以及计算分组均值后画的柱形图。

图3-2 丰度柱形图

横坐标(Sample Name)是样品名,纵坐标(Sequence Number Percent)表示MAG的丰度占总丰度的比率,柱状图自上而下的颜色顺序对应于右侧的图例颜色顺序。图例中最多显示最优势的20个MAGs,余下的相对丰度较低的MAG被归类为Other在图中展示。

2-Heatmaps文件夹中包含按照分组顺序排列的每个样本的MAG丰度热图(图3-3),以及计算分组均值后画的热图。

图3-3 丰度聚类热图

说明:纵轴为样品名称信息,同时也包括了分组信息。横轴为MAG ID。图中上方的聚类树为MAG在各样本中丰度分布的相似度聚类,左侧的聚类树为样品聚类树,中间的热图是MAGs的相对丰度(log10(TPM))热图,颜色与相对丰度的关系见图上方的刻度尺。

3-Circos 文件夹中包含按照分组顺序排列的每个样本的MAG丰度Circos图(图3-4),用于展示每个样本各个MAG(丰度前10)的丰度比例,以及每个MAG在各个样本中的丰度比例。

图3-4 丰度Circos图

说明:左半圈为丰度最高的十个MAG,每个MAG内,不同颜色代表不同样本来源的丰度比例;右边半圈为样本,每个样本内不同颜色代表不同MAG的丰度比例。

4-SignificanceAnalysis文件夹中, DunnTest文件夹包含了DunnTest组间多重比较结果表格,第一列为MAG ID, KW_pvalue是Kruskal-Wallis算出总的p值,DunnTest_comparison是比较的分组对, DunnTest_Z是DunnTest检验统计量,DunnTest_PValueAdjusted是Bofferoni校正错误发现率后的DunnTestp值; LEfSe文件夹中包含了LEfSe分析[18]中, LDA分值大于2的MAGs柱形图(图3-5), LEfSe寻找每一个分组的特征MAGs(默认为LDA>2的MAG), 也就是相对于其他分组,在这个组中丰度较高的MAG。

图3-5 LEfSe分析LDA柱形图

说明:每一横向柱形体代表一个MAG,柱形体的长度对应LDA值,LDA值越高则分组间该MAG的丰度差异越大。柱形的颜色对应该MAG是那个分组的特征MAG(在对应分组中的丰度相对较高)。

只有提供了环境因子,5-CorrelationAnalysis中才会包含内容。其中,CorrelationHeatmap文件夹中包含了相关性热图(图3-6);RDA文件夹中包含了RDA分析结果图(图3-7)。CCA/RDA的分析主要依赖R语言VEGAN包,以及用ggplot2进行可视化。CCA/RDA(DCA判断用哪一种分析)分析是基于对应分析发展的一种排序方法,将对应分析与多元回归分析相结合,每一步计算均与环境因子进行回归,又称多元直接梯度分析。RDA是基于线性模型,CCA是基于单峰模型。本报告先进行DCA分析,看最大轴的值是否大于4,如果大于4.0,就选CCA,否则选RDA。该分析主要用来反映菌群与环境因子之间的关系,可以检测环境因子、样品、菌群(抗性基因,KEGG功能)三者之间的关系或者两两之间的关系,可得到影响样品分布的重要环境驱动因子。该分析给出的所有p值都是反映解释变量(连续的数值变量,或者分类变量)对MAGs群落丰度结构的解释程度是否显著(简单的说就是解释变量对MAGs 丰度是否有影响,影响是否显著),所有p值都是用R语言VEGAN包里的置换检验得出的(permutation test),*_features_location_plot图中(图3-7)的p值反映了所有连续的数值变量(环境因子)对MAGs丰度差异的解释程度(总的p值),表格*_RDA.envfit中的p值反映了每个环境因子对MAGs丰度差异的解释程度;*_RDA_sample_location_plot图中的p值反映了分组对MAGs丰度差异的解释程度,p<0.05,解释方差显著,图中样本点之间的距离近似于样本之间MAGs丰度结构差异程度,样本点投影到环境因子对应的值近似于该样本真实的环境因子值。

图3-6 MAGs与环境因子之间的相互关系热图

说明:X轴上为环境因子,Y轴为MAGs。计算获得R值(秩相关)和校正错误发现率的P值。R值在图中以不同颜色展示,右侧图例是不同R值的颜色区间。* 0.01≤ P <0.05,** 0.001≤P < 0.01,*** P < 0.001。

图3-7 MAGs CCA/RDA排序图

说明:在CCA/RDA MAGs排序图内,环境因子用箭头表示, 箭头连线的长度代表某个环境因子与MAGs丰度分布间相关程度的大小(解释方差的大小), 箭头越长,说明相关性越大,反之越小。箭头连线和排序轴的夹角代表某个环境因子与排序轴的相关性大小, 夹角越小,相关性越高;反之越低。 环境因子之间的夹角为锐角时表示两个环境因子之间呈正相关关系,钝角时呈负相关关系。 每个点代表一个MAG,点越大,对应MAG丰度越高, 灰色点代表丰度较低的MAGs,未在图中注释MAGs名称,将MAGs投影到各个环境因子, 对应的值即为该MAGs倾向于存在的环境(喜欢的环境)。或者说,将MAGs点与原点连线,MAGs间, MAGs与环境因子间,环境因子间的夹角的余弦值近似于相关系数,至于有多近似, 就要看RDA1/CCA1和RDA2/CCA2两个坐标轴的解释方差百分比有多大,越大越近似。

04-MAG_Function MAGs功能注释情况统计

04-MAG_Function [功能注释结果文件夹]
├── CAZyme [CAZy碳水化合物数据库注释结果]
├── COG    [eggNOG数据库COG注释结果]
├── GO     [GO数据库注释结果]
├── KEGG   [KEGG数据库注释结果]
└── *
├── *_detail_MAG           [包含每个MAG 注释到数据库基因家族的种类和CDS个数统计图]
|	└── *_annotations.txt [每个MAG的原始注释结果文件]
├── *_detail_MAG_L1        [包含每个MAG 注释到数据库基因家族level1水平的种类和CDS个数统计图]
├── *_detail_gene          [包含每个注释到的数据库基因家族来源于各个MAG 的CDS个数统计图]
├── *_summary              [注释情况汇总统计图,由于图片大小限制,只会展示CDS个数排名靠前的基因家族和MAGs]
├── *_GeneCount.xls        [MAGs注释到数据库基因家族的CDS个数统计表]
└── *_GeneCount.L1.xls     [MAGs注释到数据库基因家族level1水平的CDS个数统计表]
		

使用Prokka软件[9] 对每个MAG进行功能注释,Prokka是常用的基因组注释软件, Prokka先使用Prodigal[10] 预测得到每个MAG中的基因, 得到核苷酸序列(基因序列)和氨基酸序列 (通过核苷酸序列翻译得到,也常称为蛋白序列),将序列与已有数据库信息比对, 获得每个MAG中的rRNA、tRNA、tmRNA、CDS、直系同源蛋白簇(COG)、EC等注释信息。

我们接下来对Prokka得到的每个MAG中的氨基酸序列做进一步注释。 使用KofamKOALA[11] 的离线版本Kofamscan进行KEGG功能注释。 使用emapper[12] 注释得到多个数据库的注释信息(包括GO)。 使用Diamond软件[13] 比对CAZy数据库, 对每个MAG进行碳水化合物酶注释, 获取每个MAG中的碳水化合物酶信息。

有了每个MAG中的基因组注释信息,我们就可以统计每个基因组(MAG)上各个基因家族的数目 (同一个基因家族在该基因组上注释到的次数,可理解为拷贝数)。

04-MAG_Function文件夹包含了各个功能数据库的注释情况统计;接下来我们将对每个功能数据库展开说明:

04-MAG_Function/CAZyme CAZy碳水化合物数据库注释

碳水化合物酶(CAZy)数据库是关于能够合成或者分解复杂碳水化合物和糖复合物的酶类的数据库。CAZy数据库基于蛋白质结构域中的氨基酸序列相似性将碳水化合物活性酶类归入不同蛋白质家族。CAZy数据库提供了酶分子序列的家族信息,物种来源,基因序列,蛋白序列,三维结构,EC分类,相关数据库链接。此数据库可以将酶分子的序列、结构、催化机制关联起来。CAZy有三个level,第一个level是六大功能类,即GH GT AA CE PL CBM;第二个level是CAZy family;第三个是有EC编号的具体酶信息,如EC2.4.1.129。

利用蛋白比对工具diamond将prokka预测得到的所有功能区(CDS)与CAZy数据库进行比对可获取基因组CAZyme family注释信息。

碳水化合物活性酶数据库中CAZymes功能分类如下:

名称 缩写
糖苷水解酶类 GHs
糖苷转移酶类 GTs
多糖裂解酶类 PLs
糖水化合物脂酶类 CEs
碳水化合物结合模块 CBMs
辅助模块酶类 AAs

部分结果统计图如下:

图3-8 Level1 CAZymes数量统计箱图

说明:横坐标是 CAZymes的level1大类,纵坐标是每一类CAZyme在各个MAG中的CDS数量, 每个箱子的高度表示CAZy数量的变化幅度,箱子的首和尾是分位数,中间的横线是中位数。

图3-9 Level1 CAZymes数量统计热图

说明:pheatmap热图可以展示每个MAG的CAZyme CDS计数,不同MAG之间的相似度(横向)和不同CAZyme类别之间的相似度(纵向)。图中每一个色块表示一个MAG中的一类CAZyme CDS计数,颜色越蓝表示个数越低,越红表示个数越高。

图3-10 Level1 CAZymes数量统计堆叠图

说明:一种颜色表示一类CAZyme;色柱的高度表示该CAZy的CDS计数

图3-11 Level1 CAZymes数量统计饼图

说明:该饼图展示每个MAG中的每一个CAZyme类别CDS计数与该MAG中总数量的比例。

图3-12 某个CAZyme 在每个MAG中的CDS数量统计饼图

说明:图片篇幅限制,只展示CDS数量排名前10的MAGs

04-MAG_Function/COG COGs注释

COG,即Clusters of Orthologous Groups of proteins(同源蛋白簇)。COG是由NCBI创建并维护的蛋白数据库,根据细菌、藻类和真核生物完整基因组的编码蛋白系统进化关系分类构建而成。COG分为两类,一类是原核生物的(一般称COG),另一类是真核生物(一般称KOG)。通过比对可以将某个蛋白序列注释到某一个COG中,每一簇COG由直系同源序列构成,从而可以推测该序列的功能。

图3-13 展示了COG level2水平其中一个MAG的CDS计数统计柱形图,其它MAGs请查看结果文件夹。

图3-13 COG数据库分类注释统计图

说明:该图描述的是预测基因(CDS)注释到COG不同分类水平数量情况;纵坐标是COG level 2分类(字母表示,共26种分类);横坐标是注释到的CDS计数;颜色是COG level 1分类(四大类,分别是:细胞过程和信号传递、信息储存和加工、代谢和尚未明确);柱子越长,属于该分类的预测基因越多。

COG level 2分类的26个字母的意义如下:

COG Function CategoryLevel 1Level 2
JINFORMATION STORAGE AND PROCESSINGTranslation, ribosomal structure and biogenesis
AINFORMATION STORAGE AND PROCESSINGRNA processing and modification
KINFORMATION STORAGE AND PROCESSINGTranscription
LINFORMATION STORAGE AND PROCESSINGReplication, recombination and repair
BINFORMATION STORAGE AND PROCESSINGChromatin structure and dynamics
DCELLULAR PROCESSES AND SIGNALINGCell cycle control, cell division, chromosome partitioning
YCELLULAR PROCESSES AND SIGNALINGNuclear structure
VCELLULAR PROCESSES AND SIGNALINGDefense mechanisms
TCELLULAR PROCESSES AND SIGNALINGSignal transduction mechanisms
MCELLULAR PROCESSES AND SIGNALINGCell wall/membrane/envelope biogenesis
NCELLULAR PROCESSES AND SIGNALINGCell motility
ZCELLULAR PROCESSES AND SIGNALINGCytoskeleton
WCELLULAR PROCESSES AND SIGNALINGExtracellular structures
UCELLULAR PROCESSES AND SIGNALINGIntracellular trafficking, secretion, and vesicular transport
OCELLULAR PROCESSES AND SIGNALINGPosttranslational modification, protein turnover, chaperones
XCELLULAR PROCESSES AND SIGNALINGMobilome: prophages, transposons
CMETABOLISMEnergy production and conversion
GMETABOLISMCarbohydrate transport and metabolism
EMETABOLISMAmino acid transport and metabolism
FMETABOLISMNucleotide transport and metabolism
HMETABOLISMCoenzyme transport and metabolism
IMETABOLISMLipid transport and metabolism
PMETABOLISMInorganic ion transport and metabolism
QMETABOLISMSecondary metabolites biosynthesis, transport and catabolism
RPOORLY CHARACTERIZEDGeneral function prediction only
SPOORLY CHARACTERIZEDFunction unknown

04-MAG_Function/GO GO注释

GO数据库是基因本体联合会(Gene Onotology Consortium)所建立的数据库,旨在建立一个适用于各种物种的,对基因和蛋白质功能进行限定和描述的,并能随着研究不断深入而更新的语言词汇标准。GO是多种生物本体语言中的一种,是OBO(Open BiomedicalOntologies)组织中的一员,GO提供了一系列的语义(terms)用于描绘基因、基因产物的特点,这些语义通过三个概念维度展开:细胞学组件(Cellular Component)用于描述某个节点的亚细胞结构、位置和大分子复合物,如外部封装结构(external encapsulating structure)等;分子功能(molecular function),用于描述基因以及基因产物的功能,比如蛋白质结合转录因子活性(protein binding transcription factor activity);生物学途径(biological process)指的是分子功能的有序组合以实现更复杂的生物功能,例如树突状细胞的抗原处理和提呈(dendritic cell antigen processing and presentation)。

使用eggnog-mapper功能注释软件和eggNOG数据库进行GO注释。该方法除了获得GO的注释信息还能获得KEGG、COG等信息。获得GO结果后统计注释到每个GO中的基因数量,然后进一步归类获得GO的category分类信息。最后通过R语言对该数据进行柱形图绘制。

图3-14 展示了其中一个MAG预测基因所属GO的计数,其它MAGs请查看结果文件夹。

图3-14 GO注释统计图

该图说明的是其中一个MAG预测基因所属GO的情况;横坐标是数量;纵坐标是GO;颜色是GO分类;柱子越高,该GO中包含的预测基因(CDS)越多。

04-MAG_Function/KEGG KEGG注释

KEGG数据库于 1995 年由 Kanehisa Laboratories 推出 0.1 版,目前发展为一个综合性数据库,其中最核心的为 KEGG PATHWAY 和 KEGG ORTHOLOGY 数据库。在 KEGG ORTHOLOGY 数据库中,将行使相同功能的基因聚在一起,称为 Ortholog Groups (KO entries),每个 KO 包含多个基因信息,并在一至多个 pathway 中发挥作用。而在 KEGG PATHWAY 数据库中,将生物代谢通路划分为 6 类,分别为:细胞过程(Cellular Processes)、环境信息处理(Environmental Information Processing)、遗传信息处理(Genetic Information Processing)、人类疾病(Human Diseases)、新陈代谢(Metabolism)、生物体系统(Organismal Systems),其中每类又被系统分类为二、三、四层。第二层目前包括有 57个种子 pathway;第三层即为其代谢通路图;第四层为每个代谢通路图的具体注释信息。

使用KEGG注释软件KofamScan和KOfam数据库对prokka预测基因的蛋白序列进行KEGG注释。其中一个MAG的KEGG pathway level2水平的CDS数量如图3-15所示,其它MAGs请查看结果文件夹。

图3-15 KEGG pathway level2水平CDS数量可视化

横坐标是KEGG level2水平CDS数量;纵坐标是KEGG level 2名称;颜色用于区分KEGG level 1的类型。

KEGG数据库有比较详细的通路图,我们将每个MAG注释到的KO(KEGG Ortholog Groups)在通路图中高亮显示。其中一个MAG的其中一个通路如图3-16所示,其它结果请看结果文件夹。

图3-16 MAG 注释到的基因通路图

高亮显示的基因是在该MAG中注释到的基因。对应的KO ID可以在相应的网页中查看(打开map*.html,鼠标悬浮在色块上即可)

KEGG通路的符号图例如下:

05-MAG_Taxonomy MAGs物种注释和进化分析

05-MAG_Taxonomy/                     [MAGs物种注释和进化分析结果文件夹]
├── MAGs_abundance_taxonomy.xls       [MAGs丰度表,同时在最后一列加上了MAG的物种分类注释]
├── MAGs_tree.tre                     [MAGs 进化树,nwk格式]
├── phylogenetic_tree_heatmap_ID.pdf [进化树可视化,标注MAG ID]
├── phylogenetic_tree_heatmap.pdf    [进化树可视化,标注MAG 所属的属水平分类]
└── phylogenetic_tree_table.xls      [与进化树中的热图对应的丰度表]
		

PhyloPhlAn是发表在Nature Communications上的一个用于分析微生物之间进化关系的软件, PhyloPhlAn3在PhyloPhlAn基础上作了很多改进,再次发表于Nature Communications [14]。 我们使用PhyloPhlAn3软件,基于Prokka得到的MAG中的氨基酸序列, 通过分析400多种通用微生物蛋白, 计算MAGs之间的进化树。

使用GTDB-Tk软件[15](默认参数), 将MAGs(基因组)与GTDB-Tk软件配套的 GTDB数据库比较(Release 220), 获取每个MAG的物种分类信息。

图3-17是MAG之间的进化关系树状图,图中的热图展示的是MAG丰度信息,树状图由R语言ggplot2和ggtree包绘制

图3-17 MAGs 之间的进化树

树状图中距离越近的MAG,进化关系也越亲近;用不同颜色标注出了进化树中的不同微生物门水平分类;右侧热图的颜色与MAG在各样本中的丰度对应。

06-MAG_Circos MAGs基因组Circos图

Circos是由加拿大生物信息学科学家Martin Krzywinski利用Perl语言开发的一款可以用于描述关系型数据和多维数据的基因组圈图可视化软件。2009年Circos发表在Genome research。Circos不仅能将一个物种的整个基因组呈现在一张图片中,还可以给基因组添加丰富的注释信息,如功能注释信息,差异统计信息等等。

利用每个MAG(基于contig)的Prokka蛋白预测信息,功能区注释信息(含正负链、CDS、RNA类型等信息),CAZy数据库注释结果(GH、GT、CE、PL、CBM、AA),以及GC content、GC skew的统计结果绘制Circos圈图。图3-18是其中一个MAG的基因组Circos图,其它MAG请查看结果文件夹。

图3-18 MAG基因组圈图

由于图片篇幅限制,图中只展示了长度最长的20条contigs。从外到内,第一圈表示属于该MAG的contigs,长度用刻度表示,单位为bp;第二圈用不同颜色区分contigs上的正链和负链;第三圈用不同颜色的三角标出tRNA, rRNA, tmRNA以及CDS(编码蛋白)的编码区;第四圈标注CDS注释到CAZy数据库的碳水化合物活性酶分类,如果没有注释到任何基因,该圈可能为空;第五圈标注contigs分段(每1kb)GC含量,其中用不同颜色区分大于和小于所有contigs分段GC含量总平均值;第六圈标注contigs 分段(每1kb)GC skew值,大于和小于0用不同颜色区分。


四 常见问题解答

1. 为什么大多数MAGs注释不到种水平?
    Binning是一个比较新的分析技术,与之配套的物种分类数据库还在不断发展和完善;分箱得到的MAGs,往往是存在很大变异的微生物,在数据库中还没有收录;分箱得到的MAGs,多少存在无法通过分析手段解决的序列污染,会使一个MAG注释到多个物种,从而导致最终注释结果无法精确到种水平。
2. 宏基因组分析得到的优势菌为何无法通过Binning得到其基因组?
    我们通常会觉得样本中的优势细菌更容易通过分箱得到其基因组,这种认识是不全面的。Binning算法是通过计算contigs的测序深度(同一个样本中,属于同一物种的contig, 测序深度应该近似,不同物种之间,丰度存在一定差异),结合contigs的碱基组成(属于同一个物种的contig,GC含量会比较近似,不同物种之间GC含量存在一定差异)区分不同物种的contigs,达到分箱的效果。这种算法的局限性也很明显,假如两个物种在样本中的丰度很近似,碱基组成也很近似,那分箱技术将无法区分这两个物种,从而无法分箱,或者得到一个污染严重的MAG。简而言之,不是所有的微生物的基因组都能通过分箱得到,只有那些变异比较大(通常丰度也比较小)的物种,才能够较好地和其它物种区分开,通过分箱技术得到其基因组。

五 分析所用软件的版本

软件 版本
Megahit 1.2.9
Bowtie2 2.3.5.1
Samtools 1.7
Metabat2 2.15
RefineM 0.1.2
CheckM 1.1.3
dRep 2.6.2
Metawrap 1.3.2
Prokka 1.14.6
Prodigal 2.6.3
KofamKOALA 1.2.0
emapper 2.0.1
diamond 0.9.14
PhyloPhlAn3 3.0.60
GTDB-Tk 2.4.0
Trimmomatic 0.39
LEfSe 1.0.8

六 参考文献

  • [1] D. Li, R. Luo, C.-M. Liu, C.-M. Leung, H.-F. Ting, K. Sadakane, H. Yamashita, and T.-W. Lam, “MEGAHIT v1.0: A fast and scalable metagenome assembler driven by advanced methodologies and community practices,” Methods, vol. 102, pp. 3–1 1, Jun. 2016.
  • [2] Langmead, B. , & Salzberg, S. L. . (2012). Fast gapped-read alignment with bowtie 2. Nature Methods, 9(4), 357-359.
  • [3] Li, H. , Handsaker, B. , Wysoker, A. , Fennell, T. , Ruan, J. , & Homer, N. , et al. (2009). The sequence alignment/map format and samtools. Bioinformatics, 25(16), 2078-2079.
  • [4] Kang, D. D. , Li, F. , Kirton, E. , Thomas, A. , & Wang, Z. . (2019). Metabat 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ, 7(7), e7359.
  • [5] Parks DH et al. 2017. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat Microbiol 2: 1533-1542.
  • [6] Parks, D. H. , Imelfort, M. , Skennerton, C. T. , Hugenholtz, P. , & Tyson, G. W. . (2015). Checkm: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Research, 25(7), 1043-1055.
  • [7] Olm, M. R. , Brown, C. T. , Brooks, B. , & Banfield, J. F. . (2017). Drep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. Isme Journal.
  • [8] Gherman, V, Uritskiy, Jocelyne, DiRuggiero, & James, et al. (2018). Metawrap-a flexible pipeline for genome-resolved metagenomic data analysis. Microbiome.
  • [9] Seemann T, "Prokka: Rapid Prokaryotic Genome Annotation",Bioinformatics, 2014 Jul 15;30(14):2068-9.
  • [10] Hyatt D., Chen, GL., et. Al., (2010) Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics, 11:119.
  • [11] Aramaki, T. , Blanc-Mathieu, R. , Endo, H. , Ohkubo, K. , & Ogata, H. . (2019). Kofamkoala: kegg ortholog assignment based on profile hmm and adaptive score threshold. Bioinformatics, 36(7).
  • [12] Jaime, H. C. , Kristoffer, F. , Pedro, C. L. , Damian, S. , Juhl, J. L. , & Christian, V. M. , et al. (2016). Fast genome-wide functional annotation through orthology assignment by eggnog-mapper. Molecular Biology & Evolution(8), 2115.
  • [13] Buchfink B, Xie C, Huson DH. (2015) Fast and sensitive protein alignment using DIAMOND. Nat Methods 12:59-60.
  • [14] Asnicar, F. , Thomas, A. M. , Beghini, F. , Mengoni, C. , & Segata, N. . (2020). Precise phylogenetic analysis of microbial isolates and genomes from metagenomes using phylophlan 3.0. Nature Communications, 11(1).
  • [15] Chaumeil PA, et al. 2022. GTDB-Tk v2: memory friendly classification with the Genome Taxonomy Database. Bioinformatics, btac672.
  • [16] Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics, btu170.
  • [17] Matthias Hess 1, Alexander Sczyrba, Rob Egan, et al., (2011). Metagenomic discovery of biomass-degrading genes and genomes from cow rumen. Science, 331(6016), p.463-467.
  • [18] Segata, N., et. al., (2011). Metagenomic biomarker discovery and explanation. Genome Biol.12,R60 .