技术贴 | 16S专题 | 初学者如何深入解读16S rDNA扩增子测序数据，从而选择自己的分析步骤（满满干货～）
导 读
网络上有很多16S rDNA扩增子测序数据的详细分析流程。但是很多初学者在拿到测序公司给的测序数据时，仍然不知道从何下手。究其原因，我们从测序公司拿到的数据是五花八门的，网上的分析流程虽然详细，但由于拿到的测序数据可能已经做了其中的几个步骤，而初学者很难判断自己拿到的测序数据进行到哪一步，也不知道从哪里开始。这里就教大家如何根据拿到的测序数据（fastq文件），来决定自己的分析步骤。
从16S rDNA扩增子测序原始数据到生成物种组成丰度表，是很多初学者倍感困惑的步骤。要经历（正反序列拆分）、去除barcodes、割库、（双端序列拼接）、去除引物、质控、序列聚类成OTU、挑选OTU代表序列、得到OTU丰度（频数）表（qiime2里叫feature table）、去除嵌合体、代表序列对比参考数据库得到每个OTU的微生物物种分类信息的步骤。除少数步骤外，这些步骤的总体先后顺序不是一成不变的，每一步的实现方法也不是固定的，在处理顺序和方法选择上有较大的选择空间，要根据具体数据进行选择。这里先对每个步骤做一个简要概述（后续会持续更新具体操作方法）：
正反序列拆分是把正向（forward）序列和反向（reverse）序列放到不同文件中去。原始的测序数据可能是比较混乱的，可能所有样本，所有正反序列都放在了一个文件里，如果测序公司给的数据是这样的，就需要先进行正反序列拆分，再割库。
测序前，会给每条序列加barcodes，用于区分不同的样本，一般根据序列barcodes来割库，也就是区分每条序列属于哪个样本。在去除barcodes前，一定要确保已经割库，或者将barcodes单独存放在一个fastq文件里（如extract_barcodes.py的输出文件）。在qiime1里叫split libraries，用split_libraries_fastq.py或split_libraries.py实现；qiime2里可以用dumux（割库被定义为demultiplex）插件实现，处理双端序列相应的命令是qiime demux emp-paired，它要求所有正向序列放一个文件（去除barcodes），所有反向序列放一个文件（去除barcodes），barcodes放一个文件，它们刚好是extract_barcodes.py的输出文件；如果公司给的测序数据已把不同样本的序列放到不同文件中，同一个样本正反序列也区分开了，那也可以用qiime2的qiime tools import命令，自己制作Fastq manifest文件作为输入参数，来导入文件，直接实现割库（这种方法不需要，也没有处理barcodes）。
序列拼接是将双端测序的双端序列（paired-ends），根据正反序列的重叠情况（overlap，一般判断标准是至少有20个碱基的重合）拼接成完整的序列。单端测序数据（single-ends）没有这一步。
引物是PCR扩增时为了定位目标序列加上去的，去除引物顾名思义就是把引物从测序序列中删除，以防干扰后续的OTU聚类和嵌合体去除。
质控有很多种方法，但都是大多是基于质量分的，是将测序质量较低的序列或者区间删除，也会删除测序长度较短的序列，不同分析流程（分析软件、脚本等）有自己独特的质控方式。
OTU聚类，就是根据序列相似度指定那些序列是同一种微生物（最小分类水平或者叫做操作分类水平）的序列，属于同一种微生物的所有序列就叫一个OTU（这时我们并不知道这个OTU界门纲目科属种的分类信息），同时会选出每个OTU的代表序列。这一步的方法有很多，传统的有open-reference、closed-reference等，这些方法聚类时，需要提前去重（去除重复序列，减少运算时间），它们是根据一个比较宽松的相似度标准（95%或97%）来聚类的。最近，领域内各位大神又发明了deblur、dada2、unoise等算法来聚类OTU，这些算法实现了去噪和100%相似度聚类（严格地说已经不叫聚类了，直接一步到位）生成OTU（重命名为feature，zotu等，但意思一样），用feature代替OTU是趋势。qiime2里重点推荐了deblur和dada2算法，而冷淡了closed-reference和open-reference的方法（2018.4版本的qiime2也可以实现旧版qiime1的这两个传统OTU聚类方法了）。
得到OTU丰度（频数）表。有了代表序列以后，对于每个样本，将这个样本的每条序列和代表序列比对，就可以得到该样本每条序列属于哪个OTU，从而得到了每个样本的OTU组成丰度（频数）表。当然，要是一条序列和其中多个代表序列相似，那就要考虑是不是嵌合体了。
去除嵌合体，嵌合体是由于扩增错误造成两条序列的部分嵌合成一条，不去除的话，研究者会以为自己发现了新物种。去除嵌合体也有很多方法，可以根据一个参考数据库的序列来判断一条序列是否是嵌合体，也可以根据自己数据分析得到的feature（otu）代表序列，同时结合feature丰度来判断一条序列是不是嵌合体代表序列。前者可以在OTU（feature）聚类前进行，后者只能在OTU聚类后进行。
代表序列对比参考数据库得到分类信息。参考数据库里，有参考序列，每条参考序列有对应的微生物分类信息。将自己分析得到的代表序列和参考序列比对，达到一定匹配率则认为是同一种微生物，从而得到自己OTU(Features)的分类信息。这一步在数据库上也可以有很多选择。常用的数据库有Silva和GreenGene（GG），GG数据库已经很久没有更新了，Silva数据库更新比较快，数据较详细，很受16S数据分析研究者们的欢迎。
在分析自己拿到的16S rDNA扩增子测序数据（fastq文件）之前，你需要判断你的测序数据已经进行了哪些步骤，从而决定自己还需要进行哪些步骤，或者选择合适的方法。某些软件、脚本或命令会同时包括其中几个步骤，所以要注意仔细看说明书，选择自己的分析流程（不同分析命令、软件、程序的组合）时，不要漏，也不要重复其中的步骤，否则会影响分析结果。如用qiime2的dada2分析插件聚类OTU时，需要先割库，去除barcodes和引物，但是千万不要将双端序列拼接，因为这一分析插件已经涵盖了拼接、质控、features聚类、挑选代表序列、去除嵌合体、生成feature丰度表的步骤（也可以用于粗略去除引物和barcodes，看你怎么选择参数，后续会具体讲解）。
如何判断自己数据进行到哪一步，同时获取分析需要的信息呢？这就需要详细解读一下自己拿到的fastq文件。在这之前，公司还需要给你一个map文件，对于双端测序，map文件里需要包括样本ID，LinkPrimerSequence（包括正向引物和反向引物），BarcodeSequence（包括正向barcodes，和反向barcodes），在qiime2里，这个文件也叫sample metadata，可以放入样本的分组信息，文件一般用制表符分隔，类似下面的格式：
SampleID BarcodeSequence LinkerPrimerSequence Group
Sample1 GCCAAT,CGTACG CCTAYGGGRBGCASCAG,GGACTACNNGGGTATCTAAT T
Sample2 ACAGTG,CGTACG CCTAYGGGRBGCASCAG,GGACTACNNGGGTATCTAAT T
Sample3 TGACCA,CGTACG CCTAYGGGRBGCASCAG,GGACTACNNGGGTATCTAAT C
Sample4 ATCACG,GTTTCG CCTAYGGGRBGCASCAG,GGACTACNNGGGTATCTAAT C
有的测序公司给的双端测序数据map文件给的barcodes和引物都只有正向的，这样的话，你需要找测序公司补充表格，实在没法补充，告诉你大概多长也行（去除引物和barcodes的时候，粗略去除序列左边这个长度的碱基也能达到去除引物和barcodes的效果）。
对于fastq文件，上文已经提到一些点，你需要注意你拿到的fastq文件（不包括样本map文件）是一个（所有样本、所有正向或反向序列混合在一起），两个（正向序列一个文件，反向序列一个文件，当然，我也见过虽然是两个，但是正向和反向还是混合在一起的），三个（正向序列一个文件，反向序列一个文件，barcodes一个文件），还是许多个（正向，反向，样本都分开了）。不同的情况割库方式有所区别，一个的话，要先拆分正反向序列，在提取barcodes（如extract_barcodes.py，同时去除了原序列barcodes），再割库（或拼接后割库，根据你聚类用的程序的要求来决定）；如果是两个，也要判断一下正反向是否分开放，如果分开，就提取barcodes，再割库，如果不分开，那就要先拆分；如果是多个，正反向、样本都分开放了，可以考虑qiime2，制作manifest文件来直接导入实现割库。
接下来要打开fastq文件（如果是压缩包，先解压），fastq格式的文件可以用一般文本编辑器打开（记事本，notepad++等）。
图1 fastq文件
如图1所示，在图中给的fastq文件例子里，每条序列以@HISEQ（不用@是因为在质量符号行里也会出现@符号，在用python写脚本处理fastq文件时，要特别注意这一点）开头，到下一个@HISEQ（不包括）截止（不包括），我，所以图中有5条序列。每条序列对应四行：
1）第一行从@开始，HISEQ表示测序平台是Illumina HISEQ，如果是M开头的一般是Illumina Miseq测序平台。接着是该序列在测序过程中的一些基本信息，这些信息对我们的分析过程大多没有价值，经常在分析过程中把这一行换成别的信息（如样本和序列编号），感兴趣的可以参阅http://boyun.sh.cn/bio/?p=1901。值得一提的是图中红框圈出来的数字（1或2）表示的是双端测序的正向序列或反向序列，一般双端测序正向序列文件中这个数字是1，反向序列文件中这个数字是2。我也见过一个文件中既有1也有2的（不常见），这表示这个测序数据的正向序列和反向序列是混合在一起的，这就要求分析的时候，先将正向序列和反向序列放到不同的文件中（步骤1，拆分），这个过程可以用一个简单的python脚本可以实现（感兴趣的可以在评论区留言，我根据情况决定是否在下一期介绍）。
2）第二行是最重要的测序序列，序列内是允许换行的（这样一条序列对应的就不是4行了，这里是特例，也可以用python脚本将换行的序列改成不换行的，以迎合一些程序的要求）：一般情况下，正向序列以该样本对应的正向barcodes开头，紧接着是正向引物，然后是正向序列主体；反向序列以反向barcodes开始，紧接着反向引物，然后是反向序列主体。如果是单向测序（single-ends），每个样本对应一个唯一的barcode；如果是双端测序，每个样本对应一个唯一的正向和反向barcodes组合（两个样本可能有相同的正向barcodes，但是正向和反向的barcodes组合一定是唯一的）。在一批分析数据中，所有样本使用的引物应该是相同的，否则就很难放在一起分析了。我们在拿到测序数据时，一定打开看一看序列的情况，如果fastq文件中的序列barcodes和引物都还在，这种情况就要求在分析中去除barcodes和引物，如果不去除，可能会干扰后续嵌合体的鉴定和OTU聚类（qiime2官网给的分析例子中，并没有特别指出去除barcodes和引物的步骤，这不意味我们自己的分析不需要去除）；当然我也见过一些测序公司给的fastq文件中，已经擅自去除了barcodes，这种情况下就要求这些fastq文件是进行过割库操作的（不同样本放到不同文件中），如果公司给的fastq文件中的序列已经去除barcodes，同时又没有进行割库操作（不知道序列属于哪个样本），也没有给对应的barcodes.fastq文件，那这样的测序公司，趁早解约吧；我还见过一个测序公司给的测序数据，正向序列已经去除barcodes，没有去除引物，反向序列barcodes和引物都去除了，这种情况下，我们分析时就只需要去除正向序列的引物就行。举个例子，我把图1的5条序列对齐，如图2所示，这五条序列是我一个项目中得到的测序数据的其中一个样本正向的5条序列的开始部分，公司给我的map文件告诉我这个样本（上文map文件里的Sample4）的正向barcodes是ATCACG，正向引物是CCTAYGGGRBGCASCAG（为什么会有Y，R，S这些字母，他们是混合碱基符号，表示这个位置可能是多个碱基的其中一种，如Y表示C或T，R表示A或G），我经过比较这几条序列和map文件给的barcodes和引物，发现这些序列barcodes和引物都还在，其他序列也一样，但是已经进行过割库操作（不同样本的序列放在不同文件中），我们把这样的文件导入qiime2后，就直接实现了割库，但需要单独去除barcodes和引物。
3） 第三行是“+”，后面可以跟一些注释信息，不重要，对分析过程没有影响。
4） 第四行与第二行等长，对应于第二行每个碱基的测序质量。这些五花八门的符号是美国规定的最开始的初代计算机能够使用的符号。这里稍微科普一下，很久以前，一群美国人发明了计算机，发现用八个晶体管的不同状态（包括00000000，一共可以有256种）就足以区分存储他们文化里的任意符号，每个晶体管有开和关两种状态，于是8个晶体管就相当于一个8位的2进制数字（对应十进制里的0-255），也就是一个字节，他们只编到了127（十进制）号，命名为ASCII编码表。你可以简单理解为一个符号对应于一个32-127（十进制）的数字。为什么没有0-32呢，因为0表示空字符（真正意义的空，啥也没有，不是空格，也不是制表符），1-32是一些控制字符，控制打印机的，没法显示的。那什么是测序质量呢，目前最流行的是phred质量分，计算公式Q = -10log10p，其中p是根据测序时荧光信号强度计算出的错误概率，错误概率一般是小于等于1的，所以这个值应该是大于零的，与测序错误概率负相关，一般情况下二代测序数据质量分是小于60的。到这里大家应该发现问题了，如果我们想用一个符号代表一个质量分（毕竟一个3位十进制数字还是太长，即使十六进制也得2位，二进制就想都别想了），质量分是可能在0-32范围内的，但是一个能显示的符号对应的编码值不可能是0-32的。所以为了能用一个符号代表一个质量分，我们就得把原始质量分加上33（phred33质量体系）或者64（phred64质量分体系），再来对应编码值。相应的，我们解读fsstq质量行时，对于phred33体系，就把编码值减去33得到质量分，对于phred64体系，就把符号编码值减去64得到质量分。我们在分析过程中必须要清楚自己fastq文件使用的是哪个质量体系，不然在质控的时候就乱套了。比如qiime2导入fastq时，如果用phred64参数导入phred33的数据，运行是时是就会报错的。那么怎么判断你的序列质量用的是哪个体系呢？最可靠的方法是问测序公司，看他们的测序平台使用的是哪个质量体系。要是公司不知道呢，那我们也可以根据fastq文件做一个简单判断，如图3所示是ASCII编码表，phred64体系应该是从@（对应0质量分）开始的，如果质量行符号里（可以往后一点找，序列后面测序质量会偏低）出现了编码表@前面的字符（数字、<、>、=、？、：、；等），那无疑一定是phred33体系，我这5条序列最后一条就发现了?（图1），所以一定是phred33体系。如果所有序列后面，完全找不到这些符号，那可以推测是phred64了。
图2 对齐的5条序列
图3 ASCII编码表（不包括控制符）
感谢大家的阅读，本文仅供交流和参考，如果发现我介绍的有不对的地方，请一定纠正。
本文由董小橙、江舜尧编辑。