Mothur2_减少测序和PCR错误

本文主要介绍通过 Mothur减少测序和PCR错误。

数据选择

由于原始数据集（3.9GB）很大，作者提供了362对fastq文件中的21对。例如，F3D0_S188_L001_R1_001.fastq和F3D0_S188_L001_R2_001.fastq。这两个文件对应于第0天（即断奶当天）的雌性3，第一个和所有带有R1的对应于序列1，而第二个和所有具有R2的对应于第二个或反向序列。这些序列为250 bp，在16S rRNA基因的V4区overlap，该区域长约253 bp。查看已下载的MiSeq_SOP文件夹中的文件，可以看到文件中包含22个fastq文件，代表来自雌3和“1”模拟群落的10个时间点。还包括HMP_MOCK.v35.fasta文件，文件中是以fasta格式排序的模拟群落中使用的序列。最后，有一个名为stability.files的文件，该文件的第一行如下所示：

可以使用make.file命令创建此文件。

mothur > make.file(inputdir=路径, type=fastq, prefix=stability)

第一列是每个样本的名字，第二列是正向测序序列的名字，第三列是反向测序序列的名字。最后，包含一个批处理文件。

减少测序和PCR错误

该部分要做的第一件事是完成每一个样品两端测序的拼接，然后合并所有样本的数据。这是使用make.contigs命令完成的，这个命令需要stability.files作为输入文件。此命令将从fastq文件中提取序列和质量得分数据，创建反向读取的反向互补，然后连接序列形成contigs。有一个非常简单的算法可以做到这一点。首先比对序列，接下来查看比对结果，找出两个读数不一致的位置。如果一个序列具有碱基，而另一个序列具有缺口，则该碱基的质量得分必须超过25，才能视为真实序列。如果两个序列在该位置均具有碱基，则要求其中一个碱基的质量得分比另一个要好6分或更多。如果少于6分，则将共同的碱基（consensus base）设置为N。

mothur > make.contigs(file=stability.files, processors=8)

首先处理fastq文件生成单独的fasta和qual文件。然后通查每组文件生成contigs。这在我的计算机上花费了大约31秒。显然，完整的数据集将花费更长的时间。最后显示每个样本中的序列数：

典型的序列名称将看起来像“M00967：43：000000000-A3JHG：1：1101：18327：1699”。除了异常长之外，这些序列名称还包含“：”，如果尝试利用这些序列创建系统发育树，这些序列名称将引起很多麻烦。因此，将所有“：”字符都转换为“_”字符。此命令还将生成以后需要的几个文件：stability.trim.contigs.fasta和stable.contigs.groups。这些包含每个序列的序列数据和组标识。stability.contigs.report每次读取时，该文件都含有有关contig程序集的信息。接下来使用summary.seqs命令查看这些序列是什么样的：

mothur > summary.seqs(fasta=stability.trim.contigs.fasta)

结果显示，有152360个序列，大部分在248和253个碱基之间变化。有趣的是，数据集中最长的读取时间是502 bp。应对此保持警惕：回想一下，每个读取应该是251 bp。这些片段显然没有很好地（或根本没有）组装起来。

请注意，序列中至少有2.5%不明确的碱基被识别。下一步运行screen.seqs时可以处理这些问题：

mothur > screen.seqs(fasta=stability.trim.contigs.fasta, group=stability.contigs.groups, maxambig=0, maxlength=275)

该命令将删除所有碱基不明确和长度超过275 bp的序列。还有另一种方法可以使用summary.seqs来运行它：

mothur > screen.seqs(fasta=stability.trim.contigs.fasta, group=stability.contigs.groups, summary=stability.trim.contigs.summary, maxambig=0, maxlength=275)

这可能会更快，因为summary.seqs输出文件已经为您的序列计算了模糊碱基的数量和序列长度。另外，Mothur也可以记住我们在make.contigs中使用了8个处理器，因此它将在当前的程序中使用它。要查看Mothur还记忆了什么，可以运行以下命令：

mothur > get.current()

这意味着Mothur会记住最新的fasta文件和group文件以及所拥有的处理器数量。这样就可以运行：

mothur > summary.seqs(fasta=stability.trim.contigs.good.fasta)
mothur > summary.seqs(fasta=current)
mothur > summary.seqs()

获得相同的输出。测序错误率可能下降了一个数量级以上并且有128872个序列。

这篇推文对你有帮助吗？喜欢这篇文章吗？喜欢就不要错过呀，关注本知乎号查看更多的环境微生物生信分析相关文章。亦可以用微信扫描下方二维码关注“环微分析”微信公众号，小编在里面载入了更加完善的学习资料供广大生信分析研究者爱好者参考学习，也希望读者们发现错误后予以指出，小编愿与诸君共同进步！！！

学习环境微生物分析，关注“环微分析”公众号，持续更新，开源免费，敬请关注！

转载自原创文章：

Mothur2_减少测序和PCR错误_免费生信分析课持续更新，敬请关注​

最后，再次感谢你阅读本篇文章，真心希望对你有所帮助。感谢！