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

微科盟宏基因组分箱数据分析结题报告

一概述

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

二项目流程

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^[31] 参数：ILLUMINACLIP:adapters_path:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50），和去宿主序列（Bowtie2^[32]参数：--very-sensitive），获取用于后续分析的有效序列（clean data）。测序数据预处理统计结果见表 3-1。序列质控步骤关键参数解释如下：

去除接头序列 (参数ILLUMINACLIP:adapters_path:2:30:10)；
扫描序列（4bp滑窗大小），如果平均质量分低于20（99%正确率），切除后续序列（参数SLIDINGWINDOW:4:20）；
去除最终长度小于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

Sample ID: 样品名;
InsertSize(bp): InsertSize是在建库切胶时选择的长度，合适的InsertSize能避免测序接头污染;
SeqStrategy: 测序策略（一般为双端，各150bp）;
RawReads(#)：测序Raw reads的数量;
Raw Base(GB): 以GB为单位的Raw reads数量, 测序原始数据的总碱基数，即为Raw reads数量乘以测序长度算得;
%GC：G/C碱基数占总碱基数量的百分比;
Clean Reads(#)：过滤(质控和去宿主序列)后获得的Clean reads 数量;
Cleaned(%)：过滤后剩余的序列数占Raw reads 的百分比;
Q20: 质量分高于20的碱基所占比例；
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的完整度与污染度评估]

分箱（Binning）分析步骤如下(除非特别说明, 否则各软件使用默认参数)：

组装：运用Megahit^[1] 对各样本去宿主后的clean reads 进行组装，得到样本contigs。
建索引与序列比对：运用minibwa^[3] 对每个样本的contigs分别建索引；每个样本的clean reads，需要比对自身样本的contigs，同时也需要比对组内其它样本的contigs(组内all-vs-all比对)；运用Samtools^[5]对比对结果排序并建立索引，生成BAM文件。
分箱：基于“建索引与序列比对”步骤的BAM文件计算各样本contigs的组内多样本测序深度(depth，分箱领域也叫coverage，即覆盖度)，挑选各样本长度大于1000bp的contigs(MetaBAT2除外, 它强制要求contig长度>1500bp)，运用MetaBAT2^[6] 对各样本的contigs分别进行分箱，获得原始MAGs。
MAGs提纯：运用RefineM^[10] 对MAGs进行进一步提纯，去除其中污染度较高的contigs。
质量评估与去冗余：合并所有样本的MAGs，运用CheckM^[11]的lineage_wf流程评估MAGs的完整度和污染度；运用dRep^[12]对MAGs进行去冗余（dRep参数： --completeness 90 --contamination 10 ），得到去冗余后的 MAGs 集合。
重组装：针对上一步去冗余后的 MAGs，根据前面“建索引与序列比对”步骤的比对结果BAM文件，以及各 MAG 的contigs集合信息，提取各MAG对应的特异reads(只提取MAG自身来源样本的reads，以避免破坏菌株的样本特异性)，运用MetaWRAP^[9]的reassemble_bins模块 (基于SPAdes^[13])对各MAG进行重组装。重组装完成后，再次运行 CheckM^[11] 评估这些最终 MAGs 的完整度和污染度。表3‑2 展示的即为重组装后最终 MAGs 集合的质量统计；与去冗余阶段的结果(checkm_raw.txt)相比，重组装常能使部分 MAGs 的完整度提高、污染度降低。

表3-2 最终MAGs统计

MAGID	Completeness	Contamination	N50	N90	Size	GC_percent
MAG.C2.12	89.31	2.784	6773	2350	3084828	0.520415
MAG.T11.11	95.53	3.078	36889	9301	2964067	0.383125
MAG.T11.14	87.69	5.379	14970	3445	5420119	0.662593
MAG.T11.1	95.95	1.010	277923	79559	3773525	0.498656
MAG.T11.24	93.53	1.870	62666	14464	6225202	0.626478
MAG.T11.25	90.09	1.485	23311	7024	4757249	0.496217
MAG.T11.56	94.95	1.293	107707	32654	3900513	0.681363
MAG.T11.57	96.42	0.000	19515	5452	4063474	0.387855
MAG.T11.59	92.52	1.165	18502	5492	3074489	0.428106
MAG.T11.72	94.44	5.102	32650	7917	6066297	0.663752
MAG.T11.77	91.57	1.662	20446	7422	2239752	0.656395
MAG.T12.26	91.70	0.434	41217	8279	3976735	0.624154
MAG.T12.46	99.52	0.000	130376	39923	3713881	0.454144
MAG.T12.69	94.63	0.000	198038	56006	5836656	0.514753
MAG.T13.20	96.05	0.000	46265	10517	3503545	0.487300
MAG.T13.29	92.89	2.494	89876	31751	4280487	0.670477
MAG.T13.57	99.03	1.020	112133	38719	6289436	0.611463
MAG.T13.59	98.52	0.862	240025	45685	3241034	0.372347
MAG.T13.60	95.75	1.674	21268	6307	2792083	0.402008
MAG.T13.71	99.70	1.092	185201	50675	4065898	0.421601
MAG.T13.83	91.94	2.169	12167	3448	3209872	0.366983
MAG.T13.85	91.05	3.374	16827	4635	3404888	0.425560
MAG.T13.90	93.49	1.666	37013	9074	3294058	0.727957
MAG.T21.30	94.49	0.591	63173	15313	3648595	0.644769
MAG.T21.39	90.35	5.646	26314	7460	3260498	0.692895
MAG.T21.70	95.85	2.272	37411	11363	4347203	0.594096
MAG.T22.36	96.88	3.402	63010	24210	2356144	0.635339
MAG.T22.48	93.10	0.507	79786	30043	2949438	0.553112
MAG.T22.49	98.87	0.187	177691	53214	4867090	0.463791
MAG.T23.20	94.28	1.733	41908	12484	4161403	0.697747
MAG.T23.24	91.42	1.149	6045	2538	6258264	0.541389
MAG.T23.60	96.02	1.136	322136	77184	3532390	0.642548
MAG.T23.71	90.01	0.675	11684	3344	2108079	0.446341

第一列：MAG编号；
第二列：MAG完整度；
第三列：MAG污染度；
第四列：Contig N50是指将一个MAG的Contig按照从大到小的顺序依次相加，当相加的长度达到Contig总长度的一半时，最后一个加上的Contig长度；
第五列：Contig N90是指将一个MAG的Contig按照从大到小的顺序依次相加，当相加的长度达到Contig总长度的90%时，最后一个加上的Contig长度；
第六列：一个MAG的所有Contig长度之和，即该基因组草图的总长度；
第七列：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都会在此处理（分箱）的过程中有一个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^[9] 的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分析^[34]中， 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软件^[14] 对每个MAG进行功能注释，Prokka是常用的基因组注释软件， Prokka先使用Prodigal^[15] 预测得到每个MAG中的基因，得到核苷酸序列(基因序列)和氨基酸序列（通过核苷酸序列翻译得到，也常称为蛋白序列），将序列与已有数据库信息比对，获得每个MAG中的rRNA、tRNA、tmRNA、CDS、直系同源蛋白簇（COG）、EC等注释信息。

我们接下来对Prokka得到的每个MAG中的氨基酸序列做进一步注释。使用KofamKOALA^[16] 的离线版本Kofamscan进行KEGG功能注释。使用emapper^[17] 注释得到多个数据库的注释信息(包括GO)。使用Diamond软件^[18] 比对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 Category	Level 1	Level 2
J	INFORMATION STORAGE AND PROCESSING	Translation, ribosomal structure and biogenesis
A	INFORMATION STORAGE AND PROCESSING	RNA processing and modification
K	INFORMATION STORAGE AND PROCESSING	Transcription
L	INFORMATION STORAGE AND PROCESSING	Replication, recombination and repair
B	INFORMATION STORAGE AND PROCESSING	Chromatin structure and dynamics
D	CELLULAR PROCESSES AND SIGNALING	Cell cycle control, cell division, chromosome partitioning
Y	CELLULAR PROCESSES AND SIGNALING	Nuclear structure
V	CELLULAR PROCESSES AND SIGNALING	Defense mechanisms
T	CELLULAR PROCESSES AND SIGNALING	Signal transduction mechanisms
M	CELLULAR PROCESSES AND SIGNALING	Cell wall/membrane/envelope biogenesis
N	CELLULAR PROCESSES AND SIGNALING	Cell motility
Z	CELLULAR PROCESSES AND SIGNALING	Cytoskeleton
W	CELLULAR PROCESSES AND SIGNALING	Extracellular structures
U	CELLULAR PROCESSES AND SIGNALING	Intracellular trafficking, secretion, and vesicular transport
O	CELLULAR PROCESSES AND SIGNALING	Posttranslational modification, protein turnover, chaperones
X	CELLULAR PROCESSES AND SIGNALING	Mobilome: prophages, transposons
C	METABOLISM	Energy production and conversion
G	METABOLISM	Carbohydrate transport and metabolism
E	METABOLISM	Amino acid transport and metabolism
F	METABOLISM	Nucleotide transport and metabolism
H	METABOLISM	Coenzyme transport and metabolism
I	METABOLISM	Lipid transport and metabolism
P	METABOLISM	Inorganic ion transport and metabolism
Q	METABOLISM	Secondary metabolites biosynthesis, transport and catabolism
R	POORLY CHARACTERIZED	General function prediction only
S	POORLY CHARACTERIZED	Function 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格式]
├── MAGs_taxonomy_detail.xls          [物种注释的详细情况统计]
├── MAGs_novel_species.xls            [根据GTDB-Tk分类结论推断的可能的新物种MAGs列表]
├── phylogenetic_tree_heatmap_ID.pdf [进化树可视化，标注MAG ID]
├── phylogenetic_tree_heatmap.pdf    [进化树可视化，标注MAG 所属的属水平分类]
└── phylogenetic_tree_table.xls      [与进化树中的热图对应的丰度表]

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

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

我们可以将多序列比对覆盖率 (msa_percent) > 80%, 同时GTDB-Tk分类结论为“taxonomic novelty determined using RED” 或“Genome not assigned to closest species as it falls outside its pre-defined ANI radius” 的MAGs推断为新物种(MAGs_novel_species.xls表)

图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用不同颜色区分。

四云流程个性化分析结果解读

4.1 个性化功能数据库注释

个性化功能数据库注释允许用户选择特定的功能数据库，对MAG进行功能注释。对于蛋白数据库，使用Diamond blastp软件将各MAG的蛋白序列比对到功能数据库；对于MGE核酸数据库，使用blastn软件将各MAG的基因核酸序列比对到数据库。统计比对到各功能分类的基因数目。

Diamond比对时，可通过evalue和identity两个参数控制比对的严格程度： evalue衡量的是在一个基因集（随机）中，存在比报告序列更匹配的其它序列的概率，越低越严格（默认上界0.001）； identity衡量的是查询序列和报告序列的相似程度，它=100×比对成功的长度/（查询序列长度和报告序列长度中的较小者），越大越严格（默认下界60）。

结果文件夹 01-Function_Databases 中包含所选数据库的注释结果，每个数据库一个子文件夹。每个子文件夹中，主结果文件为以数据库名命名的基因-MAG计数矩阵表（如 VFDB.txt），以及多个按功能分类层级汇总的汇总表（如 All.VFDB.Category.txt、All.VFDB.VF.Name.txt 等），还有包含详细注释信息的 detail 文件（如 VFDB.detail.txt）。所有表格均为制表符分隔的纯文本文件，数值表示对应MAG注释到该功能的基因数目，0.0表示该基因在该MAG中未检出。

各功能数据库介绍

1. VFDB（毒力因子数据库）

VFDB（Virulence Factor Database）是专门收录细菌毒力因子的数据库，涵盖病原菌中与致病性相关的基因。该数据库将毒力因子按功能分类（Category），如免疫调节（Immune modulation）、压力存活（Stress survival）、粘附（Adherence）、生物被膜（Biofilm）、效应递送系统（Effector delivery system）、外切酶（Exoenzyme）等。详情见：http://www.mgc.ac.cn/VFs/main.htm。

示例：All.VFDB.Category.txt（按毒力因子类别汇总各MAG注释到的基因数目）

Category	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
Adherence	3	8	10	4
Biofilm	0	0	1	0
Effector delivery system	1	0	5	1
Immune modulation	3	3	2	0
Stress survival	1	5	5	2

结果文件：

文件	说明
`VFDB.txt`	基因-MAG计数矩阵，行是VFDB基因ID（如VFG000036）
`VFDB.detail.txt`	含VF、Gene Name、VF Name、Gene Description、Category、Category Description、Species等详细注释
`All.VFDB.Category.txt`	按毒力因子类别汇总
`All.VFDB.VF.Name.txt`	按毒力因子名称汇总
`All.VFDB.Species.txt`	按来源物种汇总

2. CARD（抗性基因数据库）

CARD（Comprehensive Antibiotic Resistance Database）是综合性的抗生素抗性基因数据库，按抗性机制（Resistance Mechanism）分类：抗生素外排泵（antibiotic efflux）、抗生素靶标改变（antibiotic target alteration）、抗生素靶标保护（antibiotic target protection）、抗生素失活（antibiotic inactivation）等。详情见：https://card.mcmaster.ca/。

示例：All.CARD.Resistance.Mechanism.txt（按抗性机制汇总各MAG注释到的基因数目）

Resistance Mechanism	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
antibiotic efflux	3	0	10	1
antibiotic inactivation	0	0	2	0
antibiotic target alteration	0	0	1	0
antibiotic target protection	0	1	1	2
antibiotic target replacement	0	0	1	0

结果文件：

文件	说明
`CARD.txt`	基因-MAG计数矩阵，行是ARO编号（如ARO:3000165）
`CARD.detail.txt`	含Pathogen、ARO Name、AMR Gene Family、Drug Class、Resistance Mechanism、Short Name等详细注释
`All.CARD.AMR.Gene.Family.txt`	按抗性基因家族汇总
`All.CARD.Drug.Class.txt`	按药物类别汇总
`All.CARD.Resistance.Mechanism.txt`	按抗性机制汇总

3. ARDB（抗性基因数据库, VIP可用）

ARDB（Antibiotic Resistance Genes Database）是较早的抗生素抗性基因数据库，收录了各类抗生素抗性基因及其功能注释信息。详情见：http://ardb.cbcb.umd.edu/。

示例：All.ARDB.Level1.txt（按抗性基因家族汇总各MAG注释到的基因数目）

Gene Family	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1	Description
Aac2Ic	0	0	0	0	Aminoglycoside N-acetyltransferase; 6_n_netilmicin, dibekacin, gentamicin, netilmicin, tobramycin
Aac3IV	0	0	1	0	Aminoglycoside N-acetyltransferase; apramycin, dibekacin, gentamicin, netilmicin, sisomicin, tobramycin
BacA	1	1	1	1	Undecaprenyl pyrophosphate phosphatase; bacitracin
MexB	0	0	0	0	RND multidrug efflux pump; aminoglycoside, beta_lactam, glycylcycline, macrolide
tetA	0	0	1	0	MFS antibiotic efflux pump; tetracycline

结果文件：

文件	说明
`ARDB.txt`	基因-MAG计数矩阵，行是UniProt编号
`All.ARDB.Level1.txt`	按抗性基因家族（Level1）汇总

4. SARG（抗性基因数据库, VIP可用）

SARG（Structured Antibiotic Resistance Gene database）是结构化的抗生素抗性基因数据库，对抗性基因进行了多层次分类，包括Type（抗生素类型）、Subtype（亚型）、 Mechanism group（抗性机制组）、Mechanism subgroup（抗性机制亚组）等。详情见：https://smile.hku.hk/ARGs/Indexing/help。

示例：All.SARG.Type.txt（按抗生素类型汇总各MAG注释到的基因数目）

Type	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
aminoglycoside	0	0	1	0
bacitracin	0	0	0	0
beta_lactam	0	0	1	0
macrolide-lincosamide-streptogramin	0	0	0	0
multidrug	0	0	0	0

结果文件：

文件	说明
`SARG.txt`	基因-MAG计数矩阵，最后一列为抗生素类型
`SARG.detail.txt`	含Type、Subtype、HMM category、Mechanism group/subgroup/subgroup2、Pathogen等详细分类
`All.SARG.Type.txt`	按抗生素类型汇总
`All.SARG.Mechanism.group.txt`	按抗性机制汇总

5. BacMet（重金属抗性数据库）

BacMet（Antibacterial Biocide and Metal Resistance Genes Database）收录了细菌中与重金属和杀菌剂抗性相关的基因。每条记录包含基因名称、来源物种、抗性化合物等信息。详情见：http://bacmet.biomedicine.gu.se/。

示例：All.BacMet.Compound.txt（按抗性化合物汇总各MAG注释到的基因数目）

Compound	MAG.C2.12	MAG.T11.14	MAG.T11.1
Acriflavine	1	0	0
Antimony (Sb)	1	0	1
Arsenic (As)	0	2	0
Copper (Cu)	1	2	0
Zinc (Zn)	0	1	0

结果文件：

文件	说明
`BacMet.txt`	基因-MAG计数矩阵，最后一列为Organism; gene_name; Compound(s)
`BacMet.detail.txt`	含Accession、Gene Name、Organism、Compound、Location、NCBI Annotation等详细注释
`All.BacMet.Compound.txt`	按抗性化合物汇总
`All.BacMet.Gene.Name.txt`	按基因名称汇总

6. MGE（可移动遗传元件数据库）

MGE（Mobile Genetic Elements）数据库收录了各类可移动遗传元件，包括转座子（transposase）、插入序列（IS elements）、整合酶（integrase）、质粒（plasmid）等。可移动遗传元件是细菌基因组中能够进行水平转移的DNA片段，在抗生素抗性基因传播中起重要作用。详情见：https://github.com/KatariinaParnanen/MobileGeneticElementDatabase。

示例：All.MGE.Level1.txt（按MGE类型汇总各MAG注释到的基因数目）

MGE Type	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
IS1166	0	0	0	0
IS26	0	0	0	1
IS621	0	0	0	0
integrase	4	1	2	1
transposase	2	5	4	1

结果文件：

文件	说明
`MGE.txt`	基因-MAG计数矩阵，最后一列为gene_type; gene_name（如transposase; tnpA）
`All.MGE.Level1.txt`	按MGE Level1类型汇总
`All.MGE.Level2.txt`	按MGE Level2类型（如具体质粒类型IncFII、IncP等）汇总

7. PHI-base（病原-宿主互作数据库, VIP可用）

PHI-base（Pathogen-Host Interactions database）收录了经实验验证的病原菌与宿主相互作用相关的基因。该数据库提供详细的突变表型信息（Mutant Phenotype），包括致病性增强（increased virulence）、致病性减弱（reduced virulence）、无影响（unaffected pathogenicity）等。详情见：http://www.phi-base.org/。

示例：All.PHI.Disease.txt（按疾病类型汇总各MAG注释到的基因数目）

Disease	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
Fusarium ear blight	0	0	1	0
Fusarium head blight	0	1	1	1
Legionnaires' disease	1	0	1	0
salmonellosis	0	0	0	0
tuberculosis	0	0	0	0

结果文件：

文件	说明
`PHI.txt`	基因-MAG计数矩阵，最后一列为疾病名称
`PHI.detail.txt`	含Disease、PHI_MolConn_ID、Gene、Pathogen species、Host species、Gene Function、Mutant Phenotype、Exp. Technique等详细注释
`All.PHI.Disease.txt`	按疾病类型汇总
`All.PHI.Mutant.Phenotype.txt`	按突变表型汇总

8. QS（群体感应数据库, VIP可用）

QS（Quorum Sensing）数据库收录了与细菌群体感应相关的基因，主要包含酰基高丝氨酸内酯合成酶（Acyl-homoserine-lactone synthase）等群体感应信号分子合成相关蛋白。详情见：https://www.uniprot.org/keywords/KW-0673。

示例：All.QS.Gene.Name.txt（按基因名称汇总各MAG注释到的基因数目）

Gene Name	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
lasI	0	0	0	0
lasI_2	0	0	0	0
tral	0	0	0	0

结果文件：

文件	说明
`QS.txt`	基因-MAG计数矩阵
`QS.detail.txt`	含Gene Name、Protein Family、Organism、Molecular Function等详细注释
`All.QS.Gene.Name.txt`	按基因名称汇总

9. TCDB（转运蛋白分类数据库, VIP可用）

TCDB（Transporter Classification Database）是国际公认的转运蛋白分类系统，采用层次化分类体系（如1.A.1.11.15），从大类到小类逐级细分。第一级包括通道/孔蛋白（Channels/Pores）、电化学势驱动转运蛋白（Electrochemical Potential-driven Transporters）、初级主动转运蛋白（Primary Active Transporters）等。详情见：https://www.tcdb.org/。

示例：All.TCDB.Level1.txt（按TCDB一级分类汇总各MAG注释到的基因数目）

Level1	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
1: Channels/Pores	3	1	3	0
2: Electrochemical Potential-driven Transporters	5	2	17	0
3: Primary Active Transporters	20	9	44	8
4: Group Translocators	0	0	0	0
5: Transmembrane Electron Carriers	1	0	4	0

结果文件：

文件	说明
`TCDB.txt`	基因-MAG计数矩阵，行是TCDB分类编号，最后一列为层次分类描述
`TCDB.detail.txt`	含Levels、Protein ID、Superfamily、Substrates、Gene Name、Organism等详细注释
`All.TCDB.Level1.txt`	按TCDB Level1分类汇总
`All.TCDB.Superfamily.txt`	按超家族汇总

10. RemeDB（污染物降解基因数据库, VIP可用）

RemeDB（Remediation Enzymes Database）收录了与污染物降解相关的酶类基因，按降解类型分为碳氢化合物降解酶（Hydrocarbon Degrading Enzymes）、塑料降解酶（Plastic Degrading Enzymes）、染料降解酶（Dye Degrading Enzymes）等。详情见：https://niot.res.in/Remedb/。

示例：All.RemeDB.Type.txt（按降解类型汇总各MAG注释到的基因数目）

Type	MAG.T11.11	MAG.T11.14	MAG.T11.1
Dye Degrading Enzymes	0	2	0
Hydrocarbon Degrading Enzymes	0	6	1
Plastic Degrading Enzymes	1	3	1

结果文件：

文件	说明
`RemeDB.txt`	基因-MAG计数矩阵，最后一列为酶名称
`RemeDB.detail.txt`	含Description、Organism、Gene Name、Type等详细注释
`All.RemeDB.Type.txt`	按降解类型汇总

11. CCycDB（碳元素循环数据库, VIP可用）

CCycDB（Carbon Cycling Database）是专注于碳元素循环的功能基因数据库，收录了参与碳固定（Carbon Fixation）、有机碳降解（Organic Degradation）、有机碳合成（Organic Biosynthesis）等碳循环过程的基因家族。详情见：https://ccycdb.github.io/。

示例：All.CCycDB.Gene.Family.Category.txt（按Category汇总各MAG注释到的基因数目）

Category	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
Carbon Fixation	42	38	85	36
Carbon Release	39	29	84	30
Organic Biosynthesis	160	177	353	155
Organic Degration	157	183	428	144
Organic Transformation	33	35	59	26

结果文件：

文件	说明
`All.CCycDB.Gene.Family.txt`	基因家族-MAG计数矩阵，最后一列为功能分类
`CCycDB_Gene_Family.detail.txt`	含Category、subCategoryI、subCategoryII等详细分类信息
`All.CCycDB.Gene.Family.Category.txt`	按Category汇总

12. NCycDB（氮元素循环数据库, VIP可用）

NCycDB（Nitrogen Cycling Database）是专注于氮元素循环的功能基因数据库，收录了参与固氮（Nitrogen fixation）、硝化（Nitrification）、反硝化（Denitrification）、硝酸盐异化还原（Dissimilatory nitrate reduction）、硝酸盐同化还原（Assimilatory nitrate reduction）等氮循环过程的基因家族。详情见：https://github.com/qichao1984/NCyc。

示例：All.NCycDB.Gene.Family.txt（按基因家族汇总各MAG注释到的基因数目）

Gene Family	MAG.T11.1	Pathway
NR	0	Assimilatory nitrate reduction
ansB	1	Organic degradation and synthesis
asnB	0	Organic degradation and synthesis
gdh_K00260	0	Organic degradation and synthesis
glnA	0	Organic degradation and synthesis

结果文件：

文件	说明
`All.NCycDB.Gene.Family.txt`	基因家族-MAG计数矩阵，最后一列为通路信息
`NCycDB_Gene_Family.detail.txt`	含Pathway、Description等详细注释，如nifH（固氮）、narG（硝酸盐还原）
`All.NCycDB.Gene.Family.Pathway.txt`	按通路汇总

13. PCycDB（磷元素循环数据库, VIP可用）

PCycDB（Phosphorus Cycling Database）是专注于磷元素循环的功能基因数据库，收录了参与有机磷酯水解（Organic phosphoester hydrolysis）、膦酸和次膦酸代谢（Phosphonate and phosphinate metabolism）、磷转运（Transporters）等磷循环过程的基因家族。详情见：https://github.com/ZengJiaxiong/Phosphorus-cycling-database。

示例：All.PCycDB.Gene.Family.Pathway.txt（按通路汇总各MAG注释到的基因数目）

Pathway	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1
Organic phosphoester hydrolysis	1	1	1	1
Phosphonate and phosphinate metabolism	2	1	3	1
Oxidative phosphorylation	0	4	1	2
Pentose phosphate pathway	5	2	6	2
Purine metabolism	3	2	5	1

结果文件：

文件	说明
`All.PCycDB.Gene.Family.txt`	基因家族-MAG计数矩阵，最后一列为功能描述
`All.PCycDB.Gene.Family.Pathway.txt`	按通路汇总

14. SCycDB（硫元素循环数据库, VIP可用）

SCycDB（Sulfur Cycling Database）是专注于硫元素循环的功能基因数据库，收录了参与有机硫转化（Organic sulphur transformation）、异化硫还原与氧化（Dissimilatory sulphur reduction and oxidation）、同化硫酸盐还原（Assimilatory sulphate reduction）、硫氧化（Sulphur oxidation）等硫循环过程的基因家族。详情见：https://github.com/qichao1984/SCycDB。

示例：All.SCycDB.Gene.Family.txt（按基因家族汇总各MAG注释到的基因数目）

Gene Family	MAG.C2.12	MAG.T11.11	MAG.T11.14	MAG.T11.1	Pathway
acuI	0	0	2	0	Organic sulphur transformation
aprA	1	0	0	0	Dissimilatory sulphur reduction and oxidation
betB	0	2	6	3	Organic sulphur transformation
cysA	0	0	1	0	Assimilatory sulphate reduction
dmdB	0	0	0	0	Organic sulphur transformation

结果文件：

文件	说明
`All.SCycDB.Gene.Family.txt`	基因家族-MAG计数矩阵，最后一列为功能描述
`SCycDB_Gene_Family.detail.txt`	含Pathway、Description等详细注释，如dsrA（异化硫还原）、soxB（硫氧化）
`All.SCycDB.Gene.Family.Pathway.txt`	按通路汇总

4.2 水平基因转移分析

水平基因转移（Horizontal Gene Transfer, HGT）是指生物体之间通过非生殖方式传递遗传物质的过程，是微生物进化的重要驱动力，尤其在抗生素抗性基因传播、代谢通路适应等方面发挥关键作用。本模块提供两种互补的HGT分析方法：MetaCHIP水平基因转移预测和基因潜在可移动性分析。

4.2.1 MetaCHIP水平基因转移预测

MetaCHIP（Community-level Horizontal gene transfer Identification Pipeline） ^[21]是一个用于识别微生物群落中水平基因转移事件的生物信息学工具（Song et al., 2019）。该工具通过结合最佳匹配和系统发育方法，识别不同分类群之间的基因转移事件，并能推断基因转移的方向（供体 → 受体）。

分析流程：首先将所有MAG的基因核酸序列进行BLASTN比对，寻找跨MAG的最佳匹配基因对；然后通过Prodigal预测基因并构建系统发育树，检测系统发育不一致性来判断HGT事件；最后使用Ranger-DTL推断基因转移方向。用户可指定在哪个分类水平（门、纲、目、科、属、种或MAG ID）间检测HGT事件，并通过覆盖度、比对长度、相似度等阈值控制检测的严格程度。

MetaCHIP结果文件

结果文件夹 OptionalResults/02-Gene_Transfer/MetaCHIP 中包含以下文件：

文件	说明
`detected_HGTs.tsv`	检测到的所有HGT事件列表，包含供体基因（Gene_1）、受体基因（Gene_2）、序列相似性（Identity）、是否末端匹配（end_match）、是否全长匹配（full_length_match）和基因转移方向（direction）
`cir_plot_matrix.csv`	不同分类群之间的HGT数量矩阵，行=受体，列=供体
`HGTs_among_provided_groups.svg`	HGT网络弦图，展示分类群之间的基因流动，连接带宽度与HGT数量成正比
`flanking_region_plots/`	每个HGT事件的侧翼区域可视化图，展示供体和受体基因在各自contig上的排列及序列相似性
`metachip_x_detected_HGTs_donor_genes.ffn`	供体基因的核酸序列（FASTA格式）
`metachip_x_detected_HGTs_recipient_genes.ffn`	受体基因的核酸序列（FASTA格式）
`metachip_x_detected_HGTs_donor_genes.faa`	供体基因的氨基酸序列（FASTA格式）
`metachip_x_detected_HGTs_recipient_genes.faa`	受体基因的氨基酸序列（FASTA格式）
`group_file.csv`	MAG分类群分组信息，格式为"分类群,MAG名称,分类学字符串"

示例：detected_HGTs.tsv

该文件记录了所有识别到的HGT事件，关键字段包括：

Gene_1	Gene_2	Identity	end_match	full_length_match	direction
bin.T11.18_00088	bin.T13.5_00793	73.058	no	no	bin.T11.18→bin.T13.5
bin.T11.42_02412	bin.T13.51_00639	74.662	no	no	bin.T11.42→bin.T13.51
bin.T11.42_03088	bin.T11.61_01127	77.363	no	no	bin.T11.42→bin.T11.61
bin.T11.47_00528	bin.T23.8_05307	75.817	no	no	bin.T23.8→bin.T11.47

其中，Identity表示两个基因的序列相似性百分比，高于70%通常表示近期发生的HGT事件； direction（箭头方向）指示基因从供体流向受体； flanking_region_plots/文件夹中为每个HGT事件生成了对应的侧翼区域图，命名格式为供体基因___受体基因.svg，图中浅蓝色/浅绿色框表示基因编码方向，蓝色大字体突出显示HGT基因，红色条显示BLASTN比对的序列相似性。

示例：HGT网络图

HGTs_among_provided_groups.svg 以弦图形式展示不同分类群之间的基因流动网络：

每个扇形区域代表一个分类群（如Actinobacteria、Proteobacteria等）
连接带宽度与HGT数量成正比，颜色表示供体分类群
可直观观察优势HGT方向和非对称转移模式

图4-2-1 MetaCHIP水平基因转移网络弦图

说明：每个扇形区域代表一个分类群，弧长与该分类群涉及的HGT事件数量成正比。连接带表示分类群间的基因流动，带宽与HGT数量成正比，颜色对应供体分类群。可直观观察到优势基因流方向（如Verrucomicrobia向Actinobacteria的大量转移）和非对称转移模式。

cir_plot_matrix.csv 以矩阵形式展示各分类群间的HGT数量，可进一步用于统计分析和自定义可视化。该矩阵的行表示受体分类群，列表示供体分类群，数值为HGT事件数量。

示例：侧翼区域图（flanking_region_plots）

flanking_region_plots/ 文件夹中包含每个检测到的HGT事件的侧翼区域可视化图，文件命名格式为供体基因___受体基因.svg。以 bin.T11.18_00088___bin.T13.5_00793.svg 为例：

图4-2-2 MetaCHIP HGT侧翼区域图示例

说明：上方为供体contig（bin.T11.18），下方为受体contig（bin.T13.5）。浅蓝色框表示正向链编码的基因，浅绿色框表示反向链编码的基因。蓝色大字体突出显示预测为HGT的基因。红色条显示基于BLASTN结果的contig之间匹配区域的序列相似性。序列轨道左下方标注contig名称及HGT基因与contig端点的距离。图中部分基因名称为prokka软件预测结果。

本示例中共检测到77个HGT事件，对应的77个侧翼区域图均可在 flanking_region_plots/ 文件夹中查看。

4.2.2 基因潜在可移动性分析

本工具基于水平基因转移的分子机制，通过分析目标基因与移动遗传元件（Mobile Genetic Elements, MGEs）在基因组中的物理距离来评估基因的可移动性潜力。当目标基因（如抗性基因、毒力因子、重金属抗性基因等）与MGEs位于同一contig上且距离小于指定阈值（默认5 kb）时，认为该基因更有可能通过MGEs介导的机制发生水平基因转移。

分析流程：首先通过DIAMOND^[18]比对识别目标数据库基因和MGE基因在MAG上的位置；然后计算每个目标基因与同一contig上所有MGE基因的距离；保留距离小于阈值的基因-MGE对作为HGT候选；最后选取包含HGT候选最多的contig, 使用Matplotlib^[22]进行可视化图谱绘制。

支持的目标数据库包括：

CARD（综合抗生素抗性基因数据库）—— 按抗性机制分类，如抗生素外排泵、抗生素失活、抗生素靶标改变等
VFDB（毒力因子数据库）—— 按功能分类，如免疫调节、粘附、效应递送系统等
BacMet（重金属抗性数据库）—— 收录抗菌杀菌剂和重金属抗性基因
其它个性化功能数据

基因潜在可移动性分析结果文件

结果文件夹 OptionalResults/02-Gene_Transfer/Potential_HGT 中包含以下文件：

文件	说明
`{MAG_id}_hgt_results.tsv`	每个MAG的HGT候选详细信息，包含目标基因和MGE基因的ID、位置、链方向以及两者之间的距离（bp）
`{MAG_id}_hgt_{contig_id}.svg`	每个contig的HGT候选区域可视化图谱，红色箭头表示目标基因，蓝色箭头表示MGE，绿色弧线连接目标基因与其最近的MGE并标注距离

示例：TSV结果文件字段

字段	说明	示例
`contig`	基因所在的contig名称	contig_3105369
`target_gene_id`	目标基因ID（来自GFF文件）	EEINCMMK_05096
`target_db_id`	目标基因在数据库中的ID	GI1221342819
`target_start` / `target_end`	目标基因在contig上的位置	22 / 195
`target_strand`	目标基因链方向（1=正链，-1=负链）	1
`mge_gene_id`	MGE基因ID	EEINCMMK_05100
`mge_db_id`	MGE在数据库中的ID	193_IncP(6)_1__JF785550
`mge_start` / `mge_end`	MGE基因在contig上的位置	3368 / 4615
`mge_strand`	MGE基因链方向	1
`distance`	目标基因与MGE基因之间的距离（bp），距离为0表示基因嵌套在MGE内部	3173

示例：SVG图谱解读

每个SVG图谱展示一个contig上的HGT候选区域，以 bin.T23.41_hgt_contig_2909938.svg 为例：

图4-2-3 基因潜在可移动性分析contig图谱示例

说明：该图展示MAG bin.T23.41中contig_2909938上的HGT候选区域。红色箭头表示目标基因（如抗性基因、毒力因子等），箭头方向表示基因链方向（正向=正链，负向指向0=负链）。蓝色箭头表示移动遗传元件（MGEs），包括转座酶、整合酶、插入序列等。绿色弧线连接目标基因与其最近的MGE，并标注基因间距离（bp）。标题栏显示contig总长度和当前显示的基因组区域范围。

图谱元素说明：

红色箭头：目标基因（如抗性基因、毒力因子等），箭头方向表示基因链方向
蓝色箭头：移动遗传元件（MGEs），包括转座酶、整合酶、插入序列等
绿色弧线：连接目标基因与其最近的MGE，标注距离（bp）
距离为0：表示目标基因与MGE基因存在重叠，通常意味着基因嵌套在MGE内部，具有极高的水平转移潜力

本示例中，共38个MAG的110个contig生成了对应的SVG图谱，每个MAG最多选取5个包含HGT候选最多的contig进行可视化。

4.3 基因组比较分析

基因组比较分析模块提供三种互补的基因组间比较方法：Sibelia基因组共线性分析、 FastANI平均核苷酸一致性分析和POCP保守蛋白百分比分析。三种方法分别从基因组结构、核苷酸序列和蛋白质组成三个层面，全面评估基因组间的相似性和进化关系。

4.3.1 Sibelia基因组共线性分析

Sibelia（Ilia Minkin et al., 2013）^[23]是一种专门为密切相关的微生物基因组设计的共线性区块生成工具。共线性（Synteny）是指不同基因组中保守的基因或序列片段在染色体上保持相似的排列顺序和方向。通过识别基因组间的共线性区块（Synteny Blocks），可以推断基因组重排事件（如倒位、易位、插入/缺失等），评估基因组结构保守性，为进化关系研究提供重要线索。

分析流程：通过构建de Bruijn图来识别基因组间的保守区域，算法从较大的k-mer值开始，识别大的保守区块，然后逐渐减小k-mer值，识别更小的保守区域。用户可选择预设参数集（loose/fine/far）和最小共线性区块长度来控制分析的灵敏度和严格程度。分析结果以Circos圈图形式可视化展示，图中不同颜色区分不同MAG/基因组，连接线连接同源共线性区块。

Sibelia结果文件

结果文件夹 OptionalResults/03-Genome_Compare/Sibelia 中包含以下文件：

文件	说明
`blocks_coords.txt`	共线性区块坐标文件，包含序列信息和每个共线性区块在各基因组中的位置、方向和长度
`circos.svg`	Circos环形可视化图，展示基因组间共线性关系，适合论文发表
`coverage_report.txt`	覆盖率报告，统计不同Degree（出现次数）的共线性区块覆盖各基因组的比例
`d3_blocks_diagram.html`	交互式共线性区块图，支持缩放和平移查看详细信息
`genomes_permutations.txt`	基因组排列表示，将各基因组表示为共线性区块的排列顺序

示例：blocks_coords.txt

该文件分为两部分。第一部分列出输入序列的基本信息：

Seq_id	Size	Description
1	4740	contig_2559568
2	2141	contig_2559802
270	2594181	NZ_AP027055.1

第二部分详细描述每个共线性区块（Block），包含Seq_id（序列编号）、Strand（链方向，+/-）、Start/End（起止位置）、Length（区块长度）。

示例：coverage_report.txt

该文件统计不同Degree（共线性区块在输入基因组中出现的次数）的区块覆盖各基因组的比例：

Degree	Count	Total	含义
2	97	1.84%	出现2次的区块共97个，覆盖所有基因组总碱基数的1.84%
3	4	0.03%	出现3次的区块共4个
All	103	1.92%	全部103个共线性区块总覆盖率为1.92%

示例：Circos共线性环形图

circos.svg 以Circos圈图形式展示基因组间的共线性关系：

外圈表示各基因组的contig，不同颜色代表不同的MAG/基因组
内圈表示共线性区块所在的链：绿色为正向链上的区块，红色为反向链上的区块，灰色区域为非共线性区域
连接线连接同源共线性区块，相同颜色的连接线属于同一共线性区块

图4-3-1 Sibelia基因组共线性Circos图

说明：外圈不同颜色代表不同的MAG/基因组，内圈绿色表示正向链上的共线性区块，红色表示反向链上的共线性区块，灰色为非共线性区域。连接线连接同源共线性区块，相同颜色的连接线属于同一区块。本示例展示了两个MAG及一个参考基因组（NZ_AP027055.1）之间的共线性关系，共识别到103个共线性区块，总覆盖率为1.92%。

本示例中，共识别到103个共线性区块，覆盖所有基因组总碱基数的1.92%。低覆盖率可能表明基因组间存在较大的结构差异或序列分化。d3_blocks_diagram.html 提供了交互式可视化界面，可缩放和平移查看具体区块的详细信息。

4.3.2 FastANI平均核苷酸一致性分析

FastANI（Jain et al., 2018）^[24] 是一种快速计算全基因组平均核苷酸一致性（Average Nucleotide Identity, ANI）的比对自由算法。 ANI定义为两个微生物基因组之间共享直系同源区域的核苷酸一致性平均值，是微生物物种划分和基因组相似性评估的重要量化指标。与传统基于BLAST的ANI计算方法相比，FastANI在保持相同准确度的前提下，实现了2-3个数量级的速度提升。

分析流程：首先将输入基因组序列切割成可调长度的片段（默认3000 bp）；然后使用MashMap作为MinHash序列映射引擎，快速识别双向最佳匹配的直系同源区域；最后基于匹配片段计算平均核苷酸一致性。当仅比较两个基因组时，还会生成保守区域可视化图（conserved_regions.svg），以彩色连线展示两个基因组间的保守区域。

FastANI结果文件

结果文件夹 OptionalResults/03-Genome_Compare/FastANI 中包含以下文件：

文件	说明
`fastani.txt`	ANI计算结果表格，包含查询基因组、参考基因组、ANI值、共享片段数和查询基因组总片段数
`conserved_regions.svg`	保守区域可视化图，仅在两两比较时生成，展示两个基因组间的保守区域及一致性水平

示例：fastani.txt

该文件以制表符分隔，记录ANI计算结果：

query	reference	ANI(%)	shared_fragment_count	total_query_fragments
bin.T11.55	bin.T13.2	75.8841	54	773

其中，ANI(%)为平均核苷酸一致性百分比，shared_fragment_count为双向共享片段数量，total_query_fragments为查询基因组总片段数。共享片段比例 = shared_fragment_count / total_query_fragments，表示查询基因组中与参考基因组有保守对应关系的区域占比。

ANI值解读

ANI ≥ 95%：通常认为属于同一物种（同一物种内的不同菌株）
ANI 90-95%：可能属于同一属的不同物种
ANI < 90%：很可能属于不同属

本示例中，bin.T11.55 与 bin.T13.2 的ANI值为75.88%，远低于95%的物种界定阈值，共享片段比例仅约6.99%，表明二者属于不同的微生物物种。

示例：保守区域可视化图

conserved_regions.svg 展示两个基因组之间的保守区域比对情况：

上方水平条表示参考基因组（reference）
下方水平条表示查询基因组（query）
彩色连线表示保守区域，颜色从浅色（低一致性）渐变到深色（高一致性）
右下角显示比例尺

图4-3-2 FastANI保守区域可视化图

说明：上方为参考基因组（bin.T13.2），下方为查询基因组（bin.T11.55）。彩色连线连接两个基因组间的保守区域，颜色深浅反映核苷酸一致性水平（74-100%）。连线数量少且颜色较浅，表明两个基因组间的同源性有限，与ANI值75.88%的低一致性结果一致。

4.3.3 POCP保守蛋白百分比分析

POCP（Percentage of Conserved Proteins，保守蛋白百分比）通过计算两个基因组之间满足双向最佳匹配且序列相似度≥50%、覆盖度≥50%的蛋白质对占总蛋白质数量的比例，来量化基因组间的蛋白质组成相似性。其核心作用是评估微生物基因组的保守性和亲缘关系，在微生物分类学中，POCP＞50%通常作为界定同一属的重要参考指标，同时也可用于比较不同基因组的功能模块保守程度。

分析流程：首先提取两个基因组的蛋白质序列（.faa），然后通过DIAMOND BLASTP ^[18]进行双向蛋白比对，筛选满足相似度阈值和覆盖度阈值的最佳匹配蛋白对，最后计算保守蛋白占总蛋白的百分比。用户可通过设置identity阈值和覆盖度阈值来控制筛选保守蛋白的严格程度。

POCP结果文件

结果文件夹 OptionalResults/03-Genome_Compare/POCP 中包含以下文件：

文件	说明
`pocp_results.tsv`	POCP计算结果，包含两个比较的基因组名称和POCP百分比值

示例：pocp_results.tsv

Genome1	Genome2	POCP(%)
MAG.T13.59	MAG.T13.60	51.14

POCP值表示两个基因组间保守蛋白占总蛋白的百分比。本示例中，MAG.T13.59 与 MAG.T13.60 的POCP值为51.14%，略高于50%的属水平界定阈值，表明二者可能属于同一属的不同物种。POCP值越高，说明两个基因组在蛋白质组成上越相似，亲缘关系越近。

4.4 基因组特殊序列分析

基因组特殊序列分析模块提供三种基因组序列层面的深度分析方法：16S rRNA物种分类、 Tandem重复序列预测和序列查找。三种方法分别从核糖体RNA序列、重复序列结构和用户自定义序列比对三个维度，全面挖掘基因组中具有生物学意义的特殊序列特征。

4.4.1 16S rRNA物种分类

16S rRNA基因是原核生物中高度保守的管家基因，包含9个可变区（V1-V9）和保守区交替排列的结构，被广泛用于细菌和古菌的物种鉴定与系统发育分析。本模块基于prokka（内部使用Barrnap）对MAGs 进行基因预测时找到的16S rRNA序列，将其提取后使用qiime2插件进行物种注释。

分析流程：首先从Prokka^[14]结果中识别并提取所有MAGs的16S rRNA基因序列；然后使用QIIME^[25]的sklearn（基于Naive Bayes分类器）或vsearch（基于序列相似性比对）方法，将16S序列与参考数据库（SILVA 138.2或Greengenes2）进行比对，获得物种分类信息。用户可通过设置sklearn分类可信度阈值或vsearch序列比对identity阈值来控制注释的严格程度。需要注意的是，由于并非所有MAGs都含有完整的16S rRNA基因，此方法可能无法对所有MAGs进行物种注释，结果仅供参考。

16S物种分类结果文件

结果文件夹 OptionalResults/04-Genome_Seq/Classify_16S 中包含以下文件：

文件	说明
`16S_sequences.fa`	从MAGs中提取的16S rRNA基因序列（FASTA格式），序列ID格式为MAG_ID\|gene_ID
`taxonomy.tsv`	物种注释结果表，包含序列ID、分类学信息和置信度（sklearn）或一致性（vsearch）

示例：taxonomy.tsv

该文件记录了所有16S rRNA序列的物种分类结果，关键字段包括：

Feature ID	Taxon	Confidence
MAG.T11.14\|MAG.T11.14_04766	k__Bacteria; p__Pseudomonadota; c__Alphaproteobacteria; o__Acetobacterales; f__Acetobacteraceae; g__Roseomonas	0.99999
MAG.T11.24\|MAG.T11.24_01912	k__Bacteria; p__Verrucomicrobiota; c__Verrucomicrobiia; o__Verrucomicrobiales; f__DEV007; g__Incertae Sedis; s__uncultured Verrucomicrobiales bacterium	0.98279
MAG.T11.56\|MAG.T11.56_02726	k__Bacteria; p__Actinomycetota; c__Actinobacteria; o__Micrococcales; f__Microbacteriaceae; g__Agromyces	0.99988
MAG.T23.71\|MAG.T23.71_00938	k__Bacteria; p__Chlamydiota; c__Chlamydiia; o__Chlamydiales; f__Parachlamydiaceae; g__Neochlamydia; s__uncultured bacterium	0.99911

其中，Feature ID格式为"MAG_ID|gene_ID"，可以追溯到具体MAG中的16S rRNA基因； Taxon字段以层级格式（k__界; p__门; c__纲; o__目; f__科; g__属; s__种）展示物种分类信息，注释到的层级越深，说明该序列与参考数据库的匹配度越高； Confidence字段表示sklearn方法给出的分类可信度，值越接近1表示分类结果越可靠。

4.4.2 Tandem重复序列预测

Tandem Repeats Finder（TRF）^[26]是一个经典的串联重复序列检测工具（Benson, 1999），通过k-mer匹配快速发现潜在的重复单元距离，并基于统计模型（预设匹配与插入缺失概率）自动过滤和验证，最终通过动态规划生成精确比对。串联重复序列（Tandem Repeats）是指基因组中相邻排列的重复DNA序列单元，广泛存在于原核和真核生物基因组中，在染色体结构稳定性、基因表达调控、遗传疾病发生、物种特异性特征及分子标记开发等方面具有重要生物学意义。

分析流程：用户可以指定要分析的MAG，并设置多个算法参数来控制检测的灵敏度和严格程度，包括匹配权重（Match）、错配惩罚（Mismatch）、插缺惩罚（Delta）、匹配概率（PM）、插缺概率（PI）、最低分数（Minscore）和最大周期（MaxPeriod）等。 TRF程序会以HTML格式输出每个contig的重复序列检测结果，包含重复序列的概要表格和详细比对视图。

Tandem重复序列预测结果文件

结果文件夹 OptionalResults/04-Genome_Seq/Tandem 中包含以下文件：

文件	说明
`{MAG}.fa.{Match}.{Mismatch}.{Delta}.{PM}.{PI}.{Minscore}.{MaxPeriod}.summary.html`	汇总索引文件（入口文件），列出所有检测到重复序列的contig，显示每个contig的重复序列数量，点击contig名称可跳转查看该contig的详细结果
`{MAG}.fa.{contig_id}.{Match}.{Mismatch}.{Delta}.{PM}.{PI}.{Minscore}.{MaxPeriod}.{max_tr}.html`	每个contig的重复序列概要表格，包含检测到的重复序列的索引位置、周期大小、拷贝数、匹配百分比、得分和碱基组成等信息
`{MAG}.fa.{contig_id}.{Match}.{Mismatch}.{Delta}.{PM}.{PI}.{Minscore}.{MaxPeriod}.{max_tr}.txt.html`	每个contig的重复序列详细比对视图，展示每个重复单元的序列比对结果

示例：TRF重复序列检测汇总结果

以MAG.C2.12为例，summary.html 作为入口文件汇总了该MAG中所有检测到串联重复序列的contig：

图4-4-1 Tandem重复序列预测汇总索引（点击此处打开新窗口查看）

说明：该汇总表格列出了MAG.C2.12中所有检测到串联重复序列的contig，仅显示含有重复序列的contig。表格包含三列：Sequence Index（contig在MAG中的编号）、Sequence Description（contig名称）、 Number of Repeats（该contig上检测到的重复序列数量）。点击contig名称可跳转到该contig的重复序列概要表格，在概要表格中再点击Indices列的链接可进一步查看重复单元的详细序列比对。

TRF结果解读要点：

Period Size（周期大小）：重复单元的长度，微卫星一般为1-6 bp，小卫星为10-60 bp
Copy Number（拷贝数）：重复单元串联重复的次数，拷贝数越高说明重复序列越长
Percent Matches（匹配百分比）：重复单元之间的序列一致性，值越高说明重复序列越保守
Score（得分）：基于匹配、错配和插入缺失加权计算的比对得分，得分越高说明重复序列越显著
Entropy（熵值）：基于碱基频率的信息熵，范围为0-2，0表示完全由单一碱基组成，2表示四种碱基均等分布

4.4.3 序列查找

序列查找功能允许用户输入任意核苷酸序列或氨基酸序列，与所有MAGs的contigs 构建的核酸数据库进行BLAST^[27]比对，快速定位目标序列在基因组中的存在情况和位置信息。该功能可广泛应用于：查找特定功能基因（如抗性基因、毒力因子）在MAGs中的分布；验证特定基因是否存在于目标MAG中；追踪感兴趣的基因序列在宏基因组数据集中的来源。

分析流程：首先将所有MAGs的contigs序列合并，构建本地BLAST核酸数据库；然后根据用户输入的序列类型，自动选择blastn（核苷酸比对）或tblastn（氨基酸序列比对核酸数据库）进行比对；最后根据用户设置的e-value和identity阈值过滤结果，输出详细的比对信息表。

序列查找结果文件

结果文件夹 OptionalResults/04-Genome_Seq/Seq_Search 中包含以下文件：

文件	说明
`input.fa`	用户输入的查询序列（FASTA格式）
`blast_result.xls`	BLAST比对结果表，包含比对上的序列ID、相似度、比对长度、错配数、gap数、比对位置、e-value和得分等信息

示例：blast_result.xls

该文件为标准BLAST tabular输出格式（outfmt 6），各列含义如下：

qseqid	sseqid	pident	length	mismatch	gapopen	qstart	qend	sstart	send	evalue	bitscore
input_seq	MAG.T13.29\|MAG.T13.29_contig_13	61.423	998	380	2	1	994	48931	51921	0.0	1068
input_seq	MAG.T21.70\|MAG.T21.70_contig_8	71.930	57	16	0	938	994	10	180	3.92e-17	89.4
input_seq	MAG.T23.20\|MAG.T23.20_contig_23	60.413	533	206	2	354	882	1597	2	0.0	605
input_seq	MAG.T11.14\|MAG.T11.14_contig_157	61.187	219	85	0	776	994	276	932	1.61e-65	244

BLAST结果字段说明：

qseqid / sseqid：查询序列和目标序列ID，sseqid格式为"MAG_ID|contig_ID"，可追溯到具体MAG的contig
pident：序列相似性百分比，即比对区域内完全匹配的碱基/氨基酸比例
length：比对长度（bp），即查询序列与目标序列匹配上的区域长度
mismatch / gapopen：错配碱基数/氨基酸数与gap（插入缺失）数量
qstart / qend / sstart / send：比对在查询序列和数据库序列上的起始和终止位置
evalue：期望值（E-value），表示在随机情况下获得该比对或更好比对结果的期望次数，值越小越可靠，0.0表示极显著
bitscore：比特得分，基于比对统计模型计算，得分越高表示比对越可靠，可用于不同查询之间的比较

从示例结果可以看出，输入的氨基酸序列（tblastn模式）在MAG.T13.29的contig_13上获得了最长的比对区域（998 aa），e-value为0.0且bitscore最高（1068），表明该MAG最可能含有该目标基因。在不同MAG的contig上均有比对结果，说明该基因序列在多个MAG中可能存在同源序列。

4.5 基因组尺度代谢模型 (VIP 可用)

基因组尺度代谢模型（Genome-scale Metabolic Model, GEM）是通过整合基因组注释信息和生化反应数据库，构建的包含生物体所有代谢反应和代谢物的数学模型。基于GEM可进行代谢网络分析、通量平衡分析和物种互作分析，有助于理解微生物的代谢能力、营养偏好以及群落中物种间的代谢互作关系。本模块提供四种GEM相关分析：单物种基因组尺度代谢模型、群落基因组尺度代谢模型、 通量平衡分析和物种代谢互作分析。

4.5.1 单物种基因组尺度代谢模型

使用CarveMe软件（1.6.5版本）^[28]，采用自上而下的重建方法，以精心整理的通用模型为模板，依据基因组中预测的蛋白质（酶）序列，推断单个MAG的基因组尺度代谢模型，并绘制代谢网络图。用户可选择通用细菌模板（bacteria），或针对革兰氏阳性菌（grampos）、革兰氏阴性菌（gramneg）的特化模板。

单物种GEM结果文件

结果文件夹 OptionalResults/05-GEM/01-GEM_Single 中包含以下文件：

文件	说明
`{MAG}.xml`	SBML格式的代谢模型文件，包含区室列表、代谢物列表、反应列表和约束条件
`{MAG}_network.svg`	静态代谢网络图，以SVG格式呈现，同心圆布局展示不同区室的代谢物和反应关系
`{MAG}_network_interactive.html`	交互式代谢网络图，使用Plotly库实现，支持缩放、平移和悬停查看详细信息
`{MAG}_nodes.xls`	代谢物节点信息表，包含代谢物ID、名称、所在区室和连接度
`{MAG}_edges.xls`	代谢反应边信息表，包含反应物、产物、反应ID、反应名称和是否可逆

示例：代谢网络图

以MAG.T21.45为例，代谢网络图展示了该MAG的代谢反应网络：

图4-5-1 单物种代谢网络图（MAG.T21.45）

说明：图中每个节点代表一个代谢物，节点颜色区分不同区室（C_c: 胞质，C_p: 周质空间，C_e: 细胞外环境）。节点大小反映连接度（degree），连接度越高表示该代谢物参与的代谢反应越多。边代表代谢反应，箭头方向表示反应方向。图中采用同心圆布局，内圈为胞质（C_c），外圈为周质空间（C_p）和细胞外环境（C_e），仅标注连接度最高的代谢物名称。

示例：nodes.xls 和 edges.xls 字段说明

节点表字段	说明	示例
id	代谢物唯一标识符，可能来自BiGG数据库	M_10fthf_c
name	代谢物名称	10-Formyltetrahydrofolate
compartment	所在区室（C_c/C_p/C_e）	C_c
degree	连接度，即参与的代谢反应数量	12

边表字段	说明	示例
source	反应物代谢物ID	M_10fthf_c
target	产物代谢物ID	M_fprica_c
reaction	反应ID	R_AICART
reaction_name	反应名称	Phosphoribosylaminoimidazolecarboxamide formyltransferase
reversible	是否可逆反应	True

4.5.2 群落基因组尺度代谢模型

群落基因组尺度代谢模型（Community GEM）是在单物种GEM的基础上，使用CarveMe软件的 merge_community命令，合并多个MAG的代谢模型，创建共享的细胞外环境区室（ext），构建群落尺度的代谢网络，为研究种间代谢物交换和交叉喂养提供模型基础。

群落GEM结果文件

结果文件夹 OptionalResults/05-GEM/02-GEM_Community 中包含以下文件：

文件	说明
`community.xml`	合并后的群落SBML模型文件，各MAG的区室和代谢物均带有MAG后缀标识
`community_network.svg`	群落代谢网络图，采用网格+同心圆混合布局，展示群落中所有MAG的代谢网络
`community_nodes.xls`	群落代谢物节点信息表，MAG后缀标识物种归属
`community_edges.xls`	群落代谢反应边信息表，跨MAG的边表示潜在的种间代谢交互

示例：群落代谢网络图

以下示例合并了2个MAG（MAG.T21.45和MAG.T11.53），构建的群落模型包含2,942个代谢物和13,499条反应边：

图4-5-2 群落代谢网络图（2个MAG）

说明：每个MAG占据网格中的一个位置，MAG内部以同心圆布局展示（内圈为胞质C_c，外圈为周质空间C_p）。群落共享的细胞外环境（ext）代谢物节点分布在所有MAG外围的矩形框上。节点颜色区分不同区室，节点大小反映连接度。跨MAG的边表示潜在的种间代谢交互关系。

群落模型中，每个MAG的区室和代谢物均带有MAG后缀标识（如 C_c_MAG_T21_45），群落共享的 ext 区室是模拟物种间交叉喂养的关键——若某代谢物被一个MAG分泌到ext区室，而被另一个MAG从ext区室摄取，则可能表示两个物种间存在营养互作关系。 community_biomass 是群落总生物量代谢物，各MAG的生物量反应产物均汇入此节点，在本示例中其连接度高达106，是群落代谢网络的核心枢纽。

4.5.3 通量平衡分析

通量平衡分析（Flux Balance Analysis, FBA）是一种基于约束的建模方法，使用COBRApy软件（0.30.0版本）^[29]，以最大化细胞生长为目标函数，在给定的环境约束下，通过线性规划预测基因组尺度代谢模型的最优代谢通量分布。支持同时对多个MAG进行FBA分析，并对多MAG结果进行合并比较，包括碳源和微量元素利用热图。

FBA结果文件

结果文件夹 OptionalResults/05-GEM/03-GEM_FBA 中包含以下文件：

文件	说明
`{MAG}/nonzero_flux.xls`	单个MAG的非零通量反应列表，包含反应ID、名称、通量值和上下限约束
`{MAG}/flux_results.xls`	单个MAG的全部反应通量（含零通量反应）
`{MAG}/{MAG}_flux_filtered_nodes.xls`	通量过滤后的网络节点信息
`{MAG}/{MAG}_flux_filtered_edges.xls`	通量过滤后的网络边信息，包含通量值列
`{MAG}/{MAG}_flux_filtered_network.svg`	通量过滤后的代谢网络图，红色边=正向通量，绿色边=负向通量
`{MAG}/visualization.html`	Escher通路图，基于E. coli核心代谢模板展示通量分布
`merged_flux.xls`	多MAG通量合并宽表，每行一个反应，每列一个MAG
`carbon_sources_flux.xls`	各MAG在碳源交换反应上的通量值
`trace_elements_flux.xls`	各MAG在微量元素交换反应上的通量值
`carbon_sources_utilization_heatmap.svg`	碳源标准化消耗量热图
`trace_elements_utilization_heatmap.svg`	微量元素标准化消耗量热图

示例：单个MAG通量网络图

以MAG.C3.5为例，通量过滤后的代谢网络图展示了该MAG在最优通量分布下的主要代谢反应：

图4-5-3 单MAG通量网络图（MAG.C3.5）

说明：节点表示代谢物，颜色区分不同区室。边表示代谢反应，红色边表示正向通量（产物生成），绿色边表示负向通量（底物消耗）。图中仅显示非零通量的反应，标注连接度最高的代谢物名称。通量大小反映反应在最优生长条件下的活跃程度。

示例：碳源利用热图

碳源利用热图展示了各MAG对不同碳源的相对利用量，基于标准化消耗数据（正值置零→取绝对值→加1→ln对数→除以矩阵最大值）生成：

图4-5-4 碳源标准化消耗量热图

说明：行为不同碳源，列为不同MAG。颜色深浅表示标准化处理后的消耗强度，越深（越红）表示该MAG对该碳源的相对消耗量越大。通过热图可直观比较不同MAG的碳源利用偏好，偏好差异反映了生态位分化，可指导设计选择性培养基。

示例：微量元素利用热图

图4-5-5 微量元素标准化消耗量热图

说明：与碳源热图类似，展示各MAG对微量元素（如铁、钙、锰等）的相对需求差异，对培养基优化具有指导意义。

nonzero_flux.xls 解读要点

该文件是FBA分析的核心结果，需要重点关注：

Growth反应通量：FBA优化的主要目标，值越高表示生长条件越有利
EX_开头的交换反应：反映细胞与环境之间的代谢物交换，负通量表示细胞摄取该化合物（底物消耗），正通量表示细胞分泌该化合物（产物生成）。通量大小反映该化合物对细胞生长的重要程度
交叉喂养线索：若某MAG对某代谢物为正通量（分泌），而另一MAG为负通量（摄取），可能暗示两者间存在营养互作关系

4.5.4 物种代谢互作分析

SMETANA（Species METabolic iNterAction analyzer）^[30]是一种基于约束的建模方法，基于微生物的基因组尺度代谢模型，计算群落的SMETANA评分，量化种间的交叉喂养和资源竞争，预测物种互作潜力（Zelezniak et al., 2015, PNAS）。分析有全局模式和详细模式两种：

全局模式（Global）：计算MRO（Metabolic Resource Overlap，代谢资源重叠，0-1，越高竞争越激烈）和MIP（Metabolic Interaction Potential，代谢互作潜力），快速评估整个群落的互作特性
详细模式（Detailed）：计算SCS（物种耦合分数）、MUS（代谢物摄取分数）、MPS（代谢物生产分数）和SMETANA综合分数（SCS × MUS × MPS），分数越高表示交叉喂养互作越可靠

SMETANA结果文件

结果文件夹 OptionalResults/05-GEM/04-GEM_SMETANA 中根据分析模式生成以下文件之一：

文件	说明
`global.tsv`	全局模式结果，包含群落名称、培养基类型、物种数量、MIP和MRO值
`detailed.tsv`	详细模式结果，包含接收物种、供体物种、化合物、SCS、MUS、MPS和SMETANA分数

示例：全局模式结果（global.tsv）

本示例中，3个MAG的群落全局分析结果如下：

community	medium	size	mip	mro
all	complete	3	n/a	0.667

该结果表明，在完全培养基条件下，3个物种组成的群落中MRO值为0.667，表明物种间存在较高的资源竞争； MIP值为n/a，表明在完全培养基条件下物种间不存在明显的交叉喂养互作。

示例：详细模式结果（detailed.tsv）

本示例中，详细模式结果展示了物种间具体的代谢互作关系：

receiver	donor	compound	scs	mus	mps	smetana
MAG.T11.11	MAG.T11.1	M_fe2_e (Iron)	1.0	1.0	1	1.0
MAG.T11.11	MAG.T11.1	M_ile__L_e (L-Isoleucine)	1.0	1.0	1	1.0
MAG.T11.11	MAG.T11.1	M_val__L_e (L-Valine)	1.0	1.0	1	1.0
MAG.T11.24	MAG.T11.1	M_glyc_e (Glycerol)	1.0	0.88	1	0.88
MAG.T11.24	MAG.T11.1	M_nh4_e (Ammonium)	1.0	0.83	1	0.83
MAG.T11.24	MAG.T11.1	M_fe2_e (Iron)	1.0	1.0	1	1.0
MAG.T11.24	MAG.T11.1	M_alaala_e (Ala-Ala)	1.0	0.01	1	0.01

详细模式结果解读：MAG.T11.1作为供体物种，向MAG.T11.11和MAG.T11.24提供多种代谢物。其中，铁离子（M_fe2_e）的SMETANA分数在两个受体物种中均为1.0，表明MAG.T11.11和MAG.T11.24 对MAG.T11.1提供的铁离子存在高度可靠的交叉喂养依赖关系。 MAG.T11.1还为MAG.T11.11提供了异亮氨酸和缬氨酸（SMETANA=1.0），为MAG.T11.24提供了甘油和铵离子（SMETANA分别为0.88和0.83）。 SMETANA分数接近1.0表示高度可靠的交叉喂养互作，分数越低表示互作关系越弱。

4.6 其它分析

此处介绍三种云流程可用分析工具：回归分析、环形进化树和动态物种丰度桑基图，分别用于探究变量间的线性关系、MAGs进化关系与功能基因数的综合可视化，以及MAG物种分类层级间丰度流向的交互式展示。

4.6.1 回归分析

回归分析（Linear Regression）用于计算两个连续变量之间的线性回归关系，评估自变量对因变量的解释程度。分析基于R语言lm函数，对每个分组分别拟合线性回归模型 Y = βX + α，输出回归系数、决定系数R²以及显著性P值，并通过散点图加拟合直线的方式进行可视化。

用户可指定自变量表格和因变量表格（如MAGs物种丰度表、分组表中的环境因子等），选择特定的指标列（如某个MAG或某个环境因子），并可选择按分组分别进行回归分析，以观察不同分组中变量关系的差异。当生物学重复>=3且自变量和因变量组内方差不为0时，才能进行有效的回归分析。

回归分析结果文件

结果文件夹 OptionalResults/06-Other/LR 中包含以下文件：

文件	说明
`LM.svg`	回归分析散点图，包含每个分组的散点、拟合直线及回归方程

示例：回归分析结果图

LM.svg 以散点图形式展示自变量与因变量之间的线性关系：

横坐标（X轴）为自变量（如某MAG的丰度），纵坐标（Y轴）为因变量（如环境因子）
不同颜色的散点代表不同分组，每条直线对应一个分组的线性回归拟合线
图例中标注了每个分组的回归方程（Y = βX + α）、决定系数R²和P值
R²越接近1，说明自变量对因变量的解释程度越高；P值越小，回归关系越显著

图4-6-1 回归分析散点图

说明：图中展示了两个分组（C组和T组）的回归分析结果。横坐标为MAG.T11.14的丰度，纵坐标为环境因子EnvFac5。C组（橙色）回归方程为Y=-0.05742*X+6.99795，R²=0.1959，P=0.5574，表明自变量对因变量的解释程度较低且不显著；T组（蓝色）回归方程为Y=-0.03199*X+14.19568，R²=0.79997，P=0.10559，表明自变量可解释约80%的因变量变异，但样本量较小导致P值未达到显著水平。散点代表各样本的实际观测值，直线代表回归拟合结果。

4.6.2 环形进化树

环形进化树（Circular Phylogenetic Tree）是一个综合性的可视化工具，由内到外层层展示MAG的进化关系、 MAGs分组丰度比例饼图、各数据库基因数、功能类别基因数和物种分类信息。该图将进化树以环形布局呈现，可同时展示多个维度的MAG信息，便于直观比较不同MAG之间的进化关系、丰度分布和功能特征。

分析流程：首先根据MAG丰度表筛选丰度排名靠前的MAG（默认前30个），然后基于MAG进化树构建环形树图，由内向外依次叠加：（1）进化树及MAG ID标签；（2）各MAG在不同分组中的丰度比例饼图；（3）各数据库基因计数热图（如MGE、VFDB、PHI、TCDB等）；（4）外圈单独展示的功能类别基因数热图（如CARD耐药基因drug class）；（5）物种分类（门水平）瓦片。图形使用R语言的ggtree、ggtreeExtra等包绘制。

环形进化树结果文件

结果文件夹 OptionalResults/06-Other/Circul_Tree 中包含以下文件：

文件	说明
`tree_heatmap.svg`	环形进化树复合图，包含进化树、丰度饼图、基因数热图、分类瓦片
`sum_counts.xls`	各MAG中各数据库基因计数汇总表，行为MAG ID，列为数据库名称

示例：环形进化树结果图

tree_heatmap.svg 以环形复合图的形式展示多维度信息：

最内圈（进化树）：展示MAG之间的系统发育关系，树的分支越近，进化关系越亲近；MAG ID标签标注在树的外围
第二圈（丰度饼图）：每个MAG对应一个饼图，展示该MAG丰度在不同分组中的占比，颜色对应分组颜色
第三圈（基因数热图）：以多个环形热图分别展示各数据库（如MGE、VFDB、PHI、TCDB等）在不同MAG中的基因计数
第四圈（功能类别热图）：展示外圈单独选定的功能类别（如SARG耐药基因的Type）的基因数，按基因数从高到低选取Top N个类别展示
最外圈（分类瓦片）：用不同颜色标注每个MAG的门水平物种分类（如Chlamydiota等），颜色对应左侧图例中的分类标签

图4-6-2 环形进化树复合图

说明：该图以环形布局展示了MAG进化树及多维度注释信息。从内到外分别为：进化树和MAG ID标签（如MAG.T12.69、MAG.T23.24等），分组丰度比例饼图，MGE/VFDB/PHI/TCDB等数据库基因计数热图（颜色从蓝到红表示基因数从低到高），SARG耐药基因的Type计数热图，以及门水平物种分类瓦片（如Chlamydiota）。左侧图例显示了各数据库热图的颜色刻度。piecol图例中的颜色对应分组方案中的分组颜色。

4.6.3 动态物种丰度桑基图

动态物种丰度桑基图（Dynamic Species Abundance Sankey Diagram）是一个交互式可视化工具，以桑基图（Sankey Diagram）的形式展示不同分类层级之间的物种对应关系以及丰度在不同样本或分组中的分布。桑基图通过流动带的宽度来表示丰度的大小，直观地展示物种丰度从高分类层级（如门）到低分类层级（如种）的流向和变化。

分析流程：读取物种丰度表（最后一列为分号分隔的分类层级字符串，如"k__Bacteria;p__Firmicutes;c__Bacilli;..."），解析为界（kingdom）、门（phylum）、纲（class）、目（order）、科（family）、属（genus）、种（species）等层级。可添加分组信息表，按分组汇总丰度后展示不同分组间的物种丰度分布差异。使用R语言plotly包生成交互式HTML网页，支持节点拖拽、缩放、标签显示调整等功能。

动态物种丰度桑基图结果文件

结果文件夹 OptionalResults/06-Other/Sankey 中包含以下文件：

文件	说明
`sankey.html`	交互式桑基图网页，可在浏览器中打开，支持拖拽节点、缩放、下载SVG/PNG等操作
`sankey_readme.html`	桑基图使用说明文档，介绍交互操作方式和结果解读

示例：动态物种丰度桑基图

sankey.html 是一个交互式网页，打开后可通过以下方式探索数据：

节点（Node）：每个矩形块代表一个分类层级中的物种/分类单元，块的高度表示该分类单元的丰度大小
流动带（Link）：连接不同层级节点的带子，宽度表示丰度流向的大小，颜色来源于源节点的颜色
交互操作：可拖拽节点调整布局；点击节点显示来源个数和目标个数；点击流动带显示源和目标信息；支持缩放和导出SVG/PNG图片
控制面板：可调整节点标签位置（左/中/右）、标签颜色和字体大小、标题文字和颜色、图片尺寸等
支持分组展示：若添加了分组信息表，可展示不同分组中物种丰度在不同分类层级间的分布差异

图4-6-3 动态物种丰度桑基图（点击此处打开新窗口查看）

说明：该图展示了从界（kingdom）到种（species）各分类层级间的物种丰度流动关系。最左侧为界层级（Bacteria），向右依次为门（如Nitrospirota、Pseudomonadota、Bacteroidota、Planctomycetota、Actinomycetota、Verrucomicrobiota等）、纲、目、科、属、种。流动带的宽度表示丰度值大小，颜色由源节点的颜色决定。用户可以拖拽节点调整布局，点击节点或流动带查看详细信息，使用控制面板调整显示效果。图中展示了不同分类层级之间的对应关系及丰度分布，可用于观察物种在各级分类中的分布情况和物种组成的层次结构。

五常见问题解答

1. 为什么大多数MAGs注释不到种水平？

Binning是一个比较新的分析技术，与之配套的物种分类数据库还在不断发展和完善；分箱得到的MAGs，往往是存在很大变异的微生物，在数据库中还没有收录；分箱得到的MAGs，多少存在无法通过分析手段解决的序列污染，会使一个MAG注释到多个物种，从而导致最终注释结果无法精确到种水平。

2. 宏基因组分析得到的优势菌为何无法通过Binning得到其基因组？

我们通常会觉得样本中的优势细菌更容易通过分箱得到其基因组，这种认识是不全面的。Binning算法是通过计算contigs的测序深度(同一个样本中，属于同一物种的contig, 测序深度应该近似，不同物种之间，丰度存在一定差异)，结合contigs的碱基组成（属于同一个物种的contig，GC含量会比较近似，不同物种之间GC含量存在一定差异）区分不同物种的contigs，达到分箱的效果。这种算法的局限性也很明显，假如两个物种在样本中的丰度很近似，碱基组成也很近似，那分箱技术将无法区分这两个物种，从而无法分箱，或者得到一个污染严重的MAG。简而言之，不是所有的微生物的基因组都能通过分箱得到，只有那些变异比较大（通常丰度也比较小）的物种，才能够较好地和其它物种区分开，通过分箱技术得到其基因组。

六分析所用软件的版本

软件	版本
Megahit	1.2.9
MMseqs2	12b7931f517415bc69adf30d62a31701907694b0
minibwa	0.3
BWA	0.7.17
Samtools	1.7
MetaBAT2	2.15
MaxBin2	2.2.7
CONCOCT	1.1.0
MetaWRAP	1.3.2
RefineM	0.1.2
CheckM	1.1.3
dRep	2.6.2
SPAdes	3.13.0
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
MetaCHIP	1.10.13
Matplotlib	3.7.5
Sibelia	3.0.7
FastANI	1.33
QIIME	2022.2
Tandem	4.09
blast	2.12.0+
CarveMe	1.6.5
COBRApy	0.30.0
SMETANA	1.2.1
Trimmomatic	0.39
Bowtie2	2.3.5.1
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] Steinegger M and Soeding J. Clustering huge protein sequence sets in linear time. Nature Communications, doi: 10.1038/s41467-018-04964-5 (2018).

[3] Li, H., & Homer, N. (2026). Fast genomic read alignment with minibwa. arXiv. https://arxiv.org/abs/2606.15357

[4] H. Li and R. Durbin, “Fast and accurate short read alignment with Burrows-Wheeler transform,” Bioinformatics, vol. 25, no. 14, pp. 1754–1760, Jul. 2009.

[5] 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.

[6] 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.

[7] Yu-Wei Wu, Blake A. Simmons, Steven W. Singer, MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets, Bioinformatics, Volume 32, Issue 4, February 2016, Pages 605–607, https://doi.org/10.1093/bioinformatics/btv638

[8] Alneberg, J., Bjarnason, B., de Bruijn, I. et al. Binning metagenomic contigs by coverage and composition. Nat Methods 11, 1144–1146 (2014). https://doi.org/10.1038/nmeth.3103

[9] Gherman, V, Uritskiy, Jocelyne, DiRuggiero, & James, et al. (2018). MetaWRAP-a flexible pipeline for genome-resolved metagenomic data analysis. Microbiome.

[10] Parks DH et al. 2017. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat Microbiol 2: 1533-1542.

[11] 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.

[12] 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.

[13] Bankevich, Anton , et al. "SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. " #i{Journal of Computational Biology} 19.5(2012):455-477.

[14] Seemann T, "Prokka: Rapid Prokaryotic Genome Annotation",Bioinformatics, 2014 Jul 15;30(14):2068-9.

[15] Hyatt D., Chen, GL., et. Al., (2010) Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics, 11:119.

[16] 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).

[17] 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.

[18] Buchfink B, Xie C, Huson DH. (2015) Fast and sensitive protein alignment using DIAMOND. Nat Methods 12:59-60.

[19] 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).

[20] Chaumeil PA, et al. 2022. GTDB-Tk v2: memory friendly classification with the Genome Taxonomy Database. Bioinformatics, btac672.

[21] Song, W.Z., Wemheuer, B., Zhang, S., Steensen, K., & Thomas, T. (2019). MetaCHIP: community-level horizontal gene transfer identification through the combination of best-match and phylogenetic approaches. Microbiome, 7(1), 36. https://doi.org/10.1186/s40168-019-0649-y

[22] J. D. Hunter, "Matplotlib: A 2D Graphics Environment", Computing in Science & Engineering, vol. 9, no. 3, pp. 90-95, 2007.

[23] Minkin, I., Patel, A., Kolmogorov, M., Vyahhi, N., & Pham, S. (2013). Sibelia: a scalable and comprehensive synteny block generation tool for closely related microbial genomes. In Algorithms in Bioinformatics (pp. 215-229). Springer Berlin Heidelberg.

[24] Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T., & Aluru, S. (2018). High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nature Communications, 9(1), 5114.

[25] Caporaso J. G., et al., QIIME allows analysis of high-throughput community sequencing data. Nat Methods, 7(5):335-6 (2010).

[26] G. Benson, Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acid Research(1999). Vol. 27, No. 2, pp. 573-580.

[27] Altschul, S. F. , Gish, W. , Miller, W. , Myers, E. W. , & Lipman, D. J. . (1990). Basic local alignment search tool.

[28] D. Machado et al, "Fast automated reconstruction of genome-scale metabolic models for microbial species and communities", Nucleic Acids Research, gky537, 2018. doi: https://doi.org/10.1093/nar/gky537

[29] Ebrahim, A., Lerman, J.A., Palsson, B.O. et al. COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Syst Biol 7, 74 (2013). https://doi.org/10.1186/1752-0509-7-74

[30] A. Zelezniak, S. Andrejev, O. Ponomarova, D.R. Mende, P. Bork, & K.R. Patil, Metabolic dependencies drive species co-occurrence in diverse microbial communities, Proc. Natl. Acad. Sci. U.S.A. 112 (20) 6449-6454, https://doi.org/10.1073/pnas.1421834112 (2015).

[31] Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics, btu170.

[32] Langmead, B. , & Salzberg, S. L. . (2012). Fast gapped-read alignment with bowtie 2. Nature Methods, 9(4), 357-359.

[33] 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.

[34] Segata, N., et. al., (2011). Metagenomic biomarker discovery and explanation. Genome Biol.12,R60 .

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

微科盟宏基因组分箱数据分析结题报告

一 概述

二 项目流程

2.1 实验流程

图2-1 测序流程图

2.2 数据分析流程

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

三 一键化分析结果

表3-1 数据产出质量情况一览表

表3-2 最终MAGs统计

图3-1 Binned contigs(重组装前)可视化

图3-2 丰度柱形图

图3-3 丰度聚类热图

图3-4 丰度Circos图

图3-5 LEfSe分析LDA柱形图

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

图3-7 MAGs CCA/RDA排序图

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

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

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

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

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

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

图3-14 GO注释统计图

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

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

图3-17 MAGs 之间的进化树

图3-18 MAG基因组圈图

四 云流程个性化分析结果解读

各功能数据库介绍

1. VFDB（毒力因子数据库）

2. CARD（抗性基因数据库）

3. ARDB（抗性基因数据库, VIP可用）

4. SARG（抗性基因数据库, VIP可用）

5. BacMet（重金属抗性数据库）

6. MGE（可移动遗传元件数据库）

7. PHI-base（病原-宿主互作数据库, VIP可用）

8. QS（群体感应数据库, VIP可用）

9. TCDB（转运蛋白分类数据库, VIP可用）

10. RemeDB（污染物降解基因数据库, VIP可用）

11. CCycDB（碳元素循环数据库, VIP可用）

12. NCycDB（氮元素循环数据库, VIP可用）

13. PCycDB（磷元素循环数据库, VIP可用）

14. SCycDB（硫元素循环数据库, VIP可用）

4.2.1 MetaCHIP水平基因转移预测

MetaCHIP结果文件

示例：detected_HGTs.tsv

示例：HGT网络图

图4-2-1 MetaCHIP水平基因转移网络弦图

示例：侧翼区域图（flanking_region_plots）

图4-2-2 MetaCHIP HGT侧翼区域图示例

4.2.2 基因潜在可移动性分析

基因潜在可移动性分析结果文件

示例：TSV结果文件字段

示例：SVG图谱解读

图4-2-3 基因潜在可移动性分析contig图谱示例

4.3.1 Sibelia基因组共线性分析

Sibelia结果文件

示例：blocks_coords.txt

示例：coverage_report.txt

示例：Circos共线性环形图

图4-3-1 Sibelia基因组共线性Circos图

4.3.2 FastANI平均核苷酸一致性分析

FastANI结果文件

示例：fastani.txt

ANI值解读

示例：保守区域可视化图

图4-3-2 FastANI保守区域可视化图

4.3.3 POCP保守蛋白百分比分析

POCP结果文件

示例：pocp_results.tsv

4.4.1 16S rRNA物种分类

16S物种分类结果文件

示例：taxonomy.tsv

4.4.2 Tandem重复序列预测

Tandem重复序列预测结果文件

示例：TRF重复序列检测汇总结果

图4-4-1 Tandem重复序列预测汇总索引（点击此处打开新窗口查看）

4.4.3 序列查找

一概述

二项目流程

三一键化分析结果

四云流程个性化分析结果解读

五常见问题解答

六分析所用软件的版本

七参考文献