16s扩增子测序分析全流程详解
扩增子测序是一种高靶向性方法（对特定长度的PCR产物或者捕获的片段进行测序），用于分析特定基因组区域中的基因变异。PCR产品（扩增子）的超深度测序可以有效地识别变异并对其进行特征分析。总体思路是靶向地捕获目标区域，然后进行新一代测序(NGS)，分析测序结果，得到相应信息。
根据研究目的不同扩增子测序可细分为16S rDNA测序、18S rDNA测序、ITS测序及目标区域扩增子测序等。下面小编对二代16S rDNA测序分析流程做一个详细介绍。
实验流程
提取样品总DNA后，根据保守区设计得到引物，在引物末端加上测序接头，进行PCR扩增并对其产物进行纯化、定量和均一化形成测序文库，建好的文库先进行文库质检，质检合格的文库用Illumina HiSeq 2500进行测序。高通量测序（如Illumina HiSeq等测序平台）得到的原始图像数据文件，经碱基识别（Base Calling）分析转化为原始测序序列（Sequenced Reads），结果以FASTQ（简称为fq）文件格式存储，其中包含测序序列（Reads）的序列信息以及其对应的测序质量信息。
生信分析流程
根据PE reads之间的Overlap关系，将Hiseq测序得到的双端序列数据拼接 (Merge)成一条序列Tags，同时对Reads的质量和Merge的效果进行质控过滤。主要有如下3个步骤：
1) PE reads拼接：使用 FLASH v1.2.7软件，通过overlap对每个样品的 reads 进行拼接，得到的拼接序列即原始 Tags 数据(Raw Tags)；
2) Tags过滤：使用 Trimmomatic v0.33软件，对拼接得到的 Raw Tags进行过滤，得到高质量的Tags 数据(Clean Tags)；
3) 去除嵌合体：使用 UCHIME v4.2软件，鉴定并去除嵌合体序列，得到最终有效数据 (Effective Tags)。
OTU即分类操作单元，是在系统发生学研究或群体遗传学研究中，为了便于进行分析，人为给某一个分类单元（品系，种，属，分组等）设置的同一标志。
使用QIIME（version 1.8.0）软件中的 UCLUST对Tags在97%的相似度水平下进行聚类、获得OTU，并基于Silva（细菌）和UNITE（真菌）分类学数据库对OTU进行分类学注释。
图1.各样品OTU个数分布图
将OTU的代表序列与微生物参考数据库进行比对可得到每个OTU对应的物种分类信息，进而在各水平（phylum，class，order，family，genus，species）统计各样品群落组成，利用QIIME软件生成不同分类水平上的物种丰度表，再利用R语言工具绘制成样品各分类学水平下的群落结构图。
图2.物种分布图 图3.物种丰度聚类热图
用QIIME软件挑选出属分类学水平上丰度最高的OTU的序列作为代表序列，进行多重序列比对并构建系统进化树，然后通过Python语言工具绘制图形。
图4.物种系统进化树
使用KRONA对物种注释结果进行可视化展示，展示结果中，圆圈从内到外依次代表不同的分类级别，扇形的大小代表不同OTU注释结果的相对比例。
图5.KRONA示意图
Alpha多样性（Alpha diversity）反映的是单个样品物种丰度（richness）及物种多样性（diversity），有多种衡量指标：Chao1、Ace、Shannon、Simpson。Chao1和Ace指数衡量物种丰度即物种数量的多少。Shannon和Simpson指数用于衡量物种多样性，受样品群落中物种丰度和物种均匀度（Community evenness）的影响。
使用Mothur （version v.1.30） 软件，对样品Alpha多样性指数进行评估。一般为比较样品间的多样性指数，分析时会将样品所含序列数进行标准化。
稀释性曲线(Rarefaction Curve)从样本中随机抽取一定数量的序列，统计这些序列所代表的物种数目，并以序列数与物种数来构建曲线，用于验证测序数据量是否足以反映样品中的物种多样性，并间接反映样品中物种的丰富程度。
图6.稀释性曲线图
利用 Mothur 软件和R语言工具依据各样品的测序量在不同测序深度时的Shannon指数（反映样品中微生物多样性的指数）绘制Shannon多样性指数稀释曲线，以此反映各样本在不同测序数量时的微生物多样性。Shannon index越大则OTU种类越多，物种越丰富，表明样品中已涵盖绝大多数的微生物物种信息。
图7.Shannon Index曲线
使用QIIME软件进行 Beta 多样性(Beta diversity)分析，比较不同样品在物种多样性方面存在的相似程度。Beta多样性分析主要采用binary jaccard、bray curtis、weighted unifrac(限细菌)、 unweighted unifrac (限细菌)等4种算法计算样品间的距离从而获得样本间的β值。基于这4种距离矩阵展开的Beta多样性分析，主要有以下几点。
1）PCA与PCoA分析：通过一系列的特征值和特征向量进行排序，选择主要的前几位特征值，采取降维的思想，找到距离矩阵中最主要的坐标，从而观察个体或群体间的差异。
图8.PCA分析图
2）NMDS分析：将对象间的相似性或相异性数据看成点间距离的单调函数，在保持原始数据次序关系的基础上，用新的相同次序的数据列替换原始数据进行度量型多维尺度分析，用于比对样本组之间的差异。
图9.NMDS分析图
3）UPGMA分析：基于各样品序列间的进化信息的差异来计算样品间的差异，可以反映样品在进化树中是否有显著的微生物群落差异。
4）组间物种差异热图：基于OTU-Table和上述4种距离矩阵进行物种聚类，从聚类中可以了解样品之间的相似性以及各分类水平上的群落构成相似性。
图10.UPGMA聚类树
组间差异显著性分析主要用于发现不同组间具有统计学差异的Biomarker。
使用PICRUSt软件通过比对16S测序数据获得的物种组成信息，推测样本中的功能基因组成，从而分析不同样本或分组之间在功能上的差异。
首先需要对生成的OTU-table进行标准化（不同的种属菌包含的16S拷贝数不相同）；然后，通过每个OTU对应的greengene id，即可获得OTU对应的KEGG和COG家族信息，从而计算该KEGG和COG的丰度并从KEGG数据库的信息中获得Pathway、EC信息、OTU丰度计算各功能类别的丰度。
图11.组间KEGG代谢途径差异分析图
图12.COG功能分类统计图