Nature Reviews：给医生的菌群分析指南(下)

本文转载自“态昌基因”，己获授权。

Nature Reviews：给医生的菌群分析指南(上），主要从实验阶段介绍了以下7方面，不清楚的请点击链接跳转原文阅读。

1. 样本选择

2. 样本的采集

3. 样本的保存与保护剂

4. DNA提取

5.我到底应该选择多样性测序还是宏基因组测序？

6.测序区域选择

7.测序平台的选择

今天我们开始菌群分析指南(下），从分析实战开始：

1. 多样性分析和宏基因组分析的基本流程
2. 多样性分析

3. 宏基因组分析

1多样性分析和宏基因组分析的基本流程

下面的图片是针对多样性分析和宏基因组分析的基本流程。

2多样性分析

QIIME、Mothur、USEARCH等是分析16S rRNA 多样性比较常用的软件。

我们的分析流程是这样的：

① 质控和过滤

Illumina平台的下机数据根据引物序列和barcode序列拆分，而后进行数据的质控和过滤（针对质量低、片段短、含有模糊碱基的序列）。

软件：USEARCH

② OTU聚类

为了研究样本的物种的组成多样性信息，对所有样本的有效序列，以97%的相似性将序列聚类成为OTUs (Operational Taxonomic Units )。OTU是物种多样性分析中的基本分类单元。在这里需要注意的是，“OTU”的概念与“种”的概念是不同的。OTU是在系统发生学或群体遗传学研究中，为了便于进行分析，人为给某一个分类单元（品系，属，种、分组等）设置的同一标志。有些时候，一个OTU可以对应一个种，但有些时候，一个OTU只能对应一个属，或者并不能划分到某个具体的分类单元。

软件：USEARCH

③ OTU物种注释

为了得到每个OTU对应的物种分类信息，采用RDP classifier贝叶斯算法对97%相似水平的OTU代表序列进行分类学分析，并分别在各个分类水平：domain（域），kingdom（界），phylum（门），class（纲），order（目），family（科），genus（属），species（种）统计各样本的群落组成。人类的16S细菌和古菌核糖体数据库一般使用silva数据库。大鼠、小鼠数据使用RDP在线数据库。

软件： Qiime

OTU注释之后，我们会对样本内的物种信息进行分析。最主要的分析内容和所用的软件如下：

④ α-多样性分析

多样性指数计算：mothur、R

稀释曲线：mothur、R

Rank-Abundance曲线：mothur、R

Specaccum物种累积曲线：R

⑤ 物种组成分析

相对丰度图： R

优势物种Heatmap图：R软件vegan包的vegdist和hclust进行距离计算和聚类分析；样本间距离算法：Bray-Curtis，聚类方法：complete；物种间距离算法：spearman，聚类算法：complete。

⑥ β-多样性分析

样本间距离计算： 样本间的物种丰度分布差异程度可通过统计学中的距离进行量化分析，使用统计算法Euclidean，Bray-Curtis，Unweighted_unifrac，weighted_unifrac，计算两两样本间距离，获得距离矩阵，可用于后续进一步的beta多样性分析和可视化统计分析。软件：Qiime计算beta多样性距离矩阵

PCA/PCoA/NMDS分析： Qiime计算beta多样性距离矩阵，R软件vegan软件包作PCoA分析和作图。

⑦ 差异分析

随机森林分析：Qiime（http://qiime.org/scripts/supervised_learning.html?highlight=random%20forest）

ROC曲线：R

Wilcoxon秩和检验分析：R

差异菌群Heatmap分析：R

Shannon多样性指数比较盒状图：R

LEfSe分析：可以去网站http://huttenhower.sph.harvard.edu/galaxy/root?tool_id=lefse_upload上根据分类学组成对样本按照不同的分组条件进行线型判别分析（LDA），找出对样本划分产生显著性差异影响的群落或物种。

ANOSIM相似性分析：R软件vegan包或QIIME

Adonis多因素方差分析：R软件vegan包或QIIME

⑧ 关联分析

RDA/CCA分析：R软件vegan包中rda或者cca分析和作图。

OTU共表达网络分析：Cytoscape或R

⑨ 16S代谢功能预测

PICRUSt详细预测过程可参考

http://picrust.github.io/picrust/tutorials/algorithm_description.html

啧啧，这些东西对于不懂代码的我一点也不难嘛（MMP），继续往下看宏基因组分析吧

3宏基因组分析

宏基因组不是针对单个目的基因进行分析，而是针对样本中的全部基因组信息进行分析。分析过程中涉及到序列的组装问题。

① 基因的组装

组装时，使用基于De-Brujin graph原理的拼接软件SOAPdenovo（http://soap.genomics.org.cn/，Version 1.06）对优化序列进行拼接组装，根据kmer间的重叠（overlap）关系，构建De-Brujin graphs，获得contigs；把reads mapping到contig上，根据成对reads之间的双端（pair-end）关系把contigs构建成scaffolds。拼接主要参数k-mer值设范围为39-47。在scaffolds内部gap处，将scaffolds打断成新的contigs，并对大于等于500bp的contigs进行统计，从中选择最优组装结果。

② 基因预测和功能分析

基因组装之后，可以使用MetaGene (http://metagene.cb.k.u-tokyo.ac.jp/)对拼接结果中的contig进行ORF预测。选择核酸长度大于等于100bp的基因，并将其翻译为氨基酸序列。

③ 基因序列聚类

来自相同环境的样品之间有很多微生物（或基因）是共有的，不同基因的丰度在样本之间的变化可以反映样本之间的共性和差异。因此可以通过构建一个非冗余基因集（non-redundant gene catalog），来描述该类环境所有基因的整体信息。将所有样品预测出来的基因序列，用CD-HIT软件（http://www.bioinformatics.org/cd-hit/）进行聚类（参数为：95% identity、90% coverage），每个类取最长的基因作为代表序列，构建非冗余基因集。

④ 基因丰度计算

使用SOAPaligner 软件(http://soap.genomics.org.cn/)，分别将每个样品的高质量reads与非冗余基因集进行比对（95% identity），统计基因在对应样品中的丰度信息。

⑤ 物种与功能注释

1) 物种分类学注释

使用BLASTP（BLAST Version 2.2.28+，http://blast.ncbi.nlm.nih.gov/Blast.cgi）将基因集与NR数据库进行比对（BLAST比对参数设置期望值e-value为1e-5），并通过NR库对应的分类学信息数据库获得物种注释，然后使用物种对应的基因丰度总和计算该物种的丰度，并在Domain（域）、Kingdom（界）、Phylum（门）、Class（纲）、Order（目）、Family（科）、Genus（属）、Species（种）各个分类学水平上统计物种在各个样品中的丰度，从而构建相应分类学水平上的丰度谱(abundance profile)。

2) COG功能注释

使用BLASTP（BLAST Version 2.2.28+，http://blast.ncbi.nlm.nih.gov/Blast.cgi）将基因集序列与eggNOG数据库进行比对（BLAST比对参数设置期望值e-value为1e-5），获得基因对应的COG （Clusters of orthologous groups of proteins，直系同源序列聚类），然后使用COG对应的基因丰度总和计算该COG的丰度。

3) KEGG功能注释

运用BLAST（BLAST Version 2.2.28+，http://blast.ncbi.nlm.nih.gov/Blast.cgi）将基因集序列与KEGG的基因数据库（GENES）进行比对，BLAST比对参数设置期望值e-value为1e-5。根据比对结果使用KOBAS 2.0（KEGG Orthology Based Annotation System，http://kobas.cbi.pku.edu.cn/home.do）进行功能注释。使用KO、Pathway、EC、Module对应的基因丰度总和计算该功能类别的丰度。

4) 碳水化合物活性酶注释

使用CAZy数据库的对应工具hmmscan（http://hmmer.janelia.org/search/hmmscan）将基因集与CAZy数据库（CAZy，http://www.cazy.org/）进行比对，比对参数设置期望值e-value为1e-5，获得基因对应的碳水化合物活性酶注释信息，然后使用碳水化合物活性酶对应的基因丰度总和计算该碳水化合物活性酶的丰度。

5) 抗生素抗性基因注释

运用BLAST（BLAST Version 2.2.28+，http://blast.ncbi.nlm.nih.gov/Blast.cgi）将基因集与ARDB数据库（ARDB，http://ardb.cbcb.umd.edu/）进行比对，比对参数设置期望值e-value为1e-5，获得基因对应的抗生素抗性功能注释信息，然后使用抗生素抗性功能对应的基因丰度总和计算该抗生素抗性功能的丰度。

其它的像是物种与功能组成分析、物种与功能差异分析、样本间比较分析等所用的方法与软件与多样性类似。

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：生信宝典 学术图表 高分文章 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外150+ PI，1300+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”