|
本文由Bayegy根据实践经验而整理，希望对大家有帮助。 原创微文，欢迎转发转载。 导 读 明明测了100000条序列，为什么我用dada2得到的OTU丰度只有8000？是人性的扭曲还是道德的沦丧？要回答这个问题，我们先从dada2的原理说起。 正 文 Qiime2 dada2用的是一个叫dada2的R包。Dada2处理16S测序数据分为以下几步： 这一步会根据大概的序列长度、最小质量分等标准将每条序列裁剪到指定长度，移除长度较短的序列，当然也可以根据固定长度裁剪掉序列5’端的引物，barcodes等。 这一步会对所有样本所有序列去重，也就是根据一定标准，去除重复序列，保留唯一的序列集合，并计算每个碱基的平均质量分，以及每一条唯一序列的丰度（同样的序列在各样本中分别有多少条）。 这一步是dada2的核心算法。由于测序错误的存在，一个OTU的序列可能测出多个序列集合（一个丰度较高的正确序列集、和一些丰度较低的错误的序列集合）。Dada2正是运用序列的丰度、质量分、序列之间的关系等信息，更正测序错误的碱基，推测真实的序列，达到去噪的目的。 嵌合体就是多条序列嵌合的序列，要鉴定嵌合体有两种方法。一种是对比嵌合体库，这种速度会相对快速。另一种是对比样本中丰度较高的序列，如果一条序列，本身丰度小，还和多条序列相似，就可以认为是嵌合体了。Dada2用的是后者。 dada2和很多其他算法的区别就在于这一步，因为它是运用在质控和去噪之后的。为什么要这样做？dada2官方解释是，拼接后的序列，在单正向、overlap、单反向三个区域的质量分不具可比性，根据这样的质量分来去噪可能带来较大误差。Qiime2 Dada2和其他算法比，拼接也是最严的，需要20bp的overlap，还不允许有错配，关键是参数还不能改，多么傲娇。  为什么要这样，dada2官方说去噪已经基本保证没有测序错误了，当然要不能有错配啦~ 看到这，PE250的小伙伴们眼泪流了下来(╥╯^╰╥) 很多小伙伴可能会疑问，为什么没有OTU聚类的步骤，那它的OTU丰度表和代表序列怎么来的。实际上第三步已经产生了OTU丰度表和代表序列（未拼接的）了，因为去噪后，假定的是所有序列都没有测序错误了，所以不需要按照97%的相似度来生成OTU，直接一条唯一序列一个OTU。当然，dada2是没有去除线粒体、叶绿体序列的步骤的，要去除线粒体、叶绿体的序列，需要另辟蹊径。 接下来就可以解释为什么100000条序列产生8000丰度的原因了。其实相信很多小伙伴已经猜出来了，症结在拼接。来来来，我们算笔账，V3+V4是从341F到806R（不同公司用的引物可能存在细微差别），一共465bp（包括引物）的长度，算上20bp的overlap，我们至少需要正向、反向共485bp的序列（包括引物）才能顺利拼接。如果我们用PE250策略测V3-V4区，也就是正向、反向各测250bp，真实测下来平均大概也就能测245bp左右，减去双端barcodes（大概6bp左右），245+245-6*2=478。不好意思，完美错过。即使--p-trunc-len-f/--p-trunc-len-r参数选择保留整条序列（以0作为参数值），不进行裁剪，还是连最低标准485都达不到，你也就只能期待PE250偶尔爆发一下测个250bp以上，还能勉强拼接一下。 这个问题目前在qiime2 dada2基本是无解的。如果你用的是PE250测序策略测V3-V4区，想用qiime2来分析数据。其一，你可以选择先根据一个较低的标准（如10bp的overlap，允许一定的错配碱基数）顺利完成拼接，再把拼接好的序列当成单端序列用dada2分析，当然，dada2官方是不推荐这么做的，因为提前拼接对去噪不利。其二，你也可以选择只分析正向序列，或者只分析反向序列，这样得到的丰度一定很高，OTU的精度也不会低太多，但这样就丢弃了大量测序数据了。其三，你可以选择抛弃傲娇的dada2，用一个较低的标准拼接，再用其他方法聚类生产OTU，比如qiime2 Deblur方法。 如果你用PE250测V3-V4区，推荐官网的这个流程：https://docs./2018.6/tutorials/read-joining/。注意修改一下qiime vsearch join-pairs的--p-minovlen（控制overlap的碱基数）和—p-maxdiffs（控制错配碱基数）参数就行。 如果你用了PE300策略测V3-V4区，或者用PE250测V4，那满足dada2还是绰绰有余的。最后附上dada2双端测序分析的一个例子： 介绍几个常用的参数，剩下的可以用—help大法解决： --p-trim-left-f /--p-trim-left-r用来去除序列5’端的引物、barcodes、和低质量序列， 三者长度之和加起来作为参数即可。 --p-trunc-len-f/--p-trunc-len-r，用最终保留序列的长度作为参数，用来去除序列后端的低质量区间，尽可能去除，注意不要去过了，保留长度太低，会导致不能完成拼接；也不要保留太长，因为达不到长度标准的序列会被丢弃，长度标准太高，得到的丰度也会降低。 --p-n-threads如果用0做为参数，能调动计算机的所有核来运行程序，加快运算速度，很好用的。 --verbose用来显示分析进度，能让你漫长的等待不是那么无聊。 单端用的命令是qiime dada2 denoise-single，大同小异，官网”moving pictures”的例子就是单端https://docs./2018.6/tutorials/moving-pictures/。 |
|
|
