本人在读研究生，方向环境微生物。之前在学习生物信息分析过程中在网络上四处奔走获取相关学习资料与解决问题，好生麻烦。于是，我就把与同学一起做的一些生物信息分析相关教程与经验总结搬运到这个CSDN这个大平台上来，希望能够与大家一起学习讨论。班门弄斧，大神见文多指教，抱拳抱拳抱拳抱拳！
本文主要介绍使用Mothur处理改进的序列。
Mothur处理改进的序列
1 unique序列
许多序列往往是彼此重复的。将同一事物对比无数次在计算上是浪费的，所以使用unique.seqs命令来唯一化序列：
mothur > unique.seqs(fasta=stability.trim.contigs.good.fasta)
如果两个序列被识别为相同序列，则它们被视为重复序列并合并。在屏幕中输出两列，第一列是特征序列的数量，第二列是剩余的唯一序列的数量。
所以在运行unique.seqs之后，序列从128872变为16426，这将使接下来的工作更加轻松。简化名称和组文件也会使得工作简单化。如果查看这些文件的最新版本，可以看到总共有13MB。看起来似乎并不多，但是在完整的MiSeq运行下，长序列名称可能会累加起来使工作变得枯燥。因此，可以运行count.seqs生成一个表，其中行是唯一序列的名称，列是组的名称。然后，用每个唯一序列在每个组中出现的次数填充表格。
mothur>count.seqs(name=stability.trim.contigs.good.names,group=stability.contigs.good.groups)
这将生成一个名为 stable.trim.contigs.good.count_table的文件。
通过使用count选项来使用上述获取的table：
mothur>summary.seqs(fasta=stability.trim.contigs.good.unique.fasta,count=stability.trim.contigs.good.count_table)
2 数据库比对
现在需要将序列与参考序列比对。同样可以通过使用pcr.seqs命令针对感兴趣的区域定制一个数据库。要运行此命令，需要有参考数据库（silva.bacteria.fasta），并知道序列比对的起始和结束的位置。为了删除leading and trailing dots，将keepdots设置为false。也可以使用感兴趣的引物来运行此命令：
mothur > pcr.seqs(fasta=silva.bacteria.fasta, start=11894, end=25319, keepdots=F, processors=8)
使用rename.file命令将其重命名为更有用的名称：
mothur > rename.file(input=silva.bacteria.pcr.fasta, new=silva.v4.fasta)
看看此命令做了什么：
mothur > summary.seqs(fasta=silva.v4.fasta)
mothur>align.seqs(fasta=stability.trim.contigs.good.unique.fasta,reference=silva.v4.fasta)
这应该在几秒钟内完成，之后可以再次运行summary.seqs命令：
mothur>summary.seqs(fasta=stability.trim.contigs.good.unique.align,count=stability.trim.contigs.good.count_table)
这是什么意思？可以看到大部分序列从1968位置开始，在11550位置结束。一些序列从1250或1968位置开始，在10693或13400位置结束。这些偏离位置的情况可能是由于在路线的末端插入或删除。有时会看到序列在同一个位置开始和结束，这表明比对非常差，这通常是由于非特异性扩增所致。为了确保所有内容都overlaps在同一区域上，可以重新运行screen.seqs以获得在1968位置或其之前开始并在11550位置或之后结束的序列。将最大的均聚物长度设置为8，因为数据库中没有任何一个连续9个或更多相同的碱基(这在上面的第一步执行screen.seqs时完成)。
请注意，需要count表方便更新去除了一定序列的表，而且还使用summary文件，这样就不必再次计算出所有的开始和停止位置：
mothur>screen.seqs(fasta=stability.trim.contigs.good.unique.align,count=stability.trim.contigs.good.count_table,summary=stability.trim.contigs.good.unique.summary, start=1968, end=11550, maxhomop=8)
mothur > summary.seqs(fasta=current, count=current)
现在知道序列与相同的比对坐标重叠，但想确保它们仅与该区域重叠。接下来将过滤序列以消除两端的overhangs。由于已经完成了配对末端测序，因此这不是什么大问题，但是以防万一。此外，对齐中还有许多列只包含间隙字符(即‘-’)。这些都可以在不丢失任何信息的情况下取出。将用filter.seqs来做这一切，运行命令结束时，将获得以下信息：
mothur > filter.seqs(fasta=stability.trim.contigs.good.unique.good.align, vertical=T, trump=.)
mothur>unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta, count=stability.trim.contigs.good.good.count_table)
3 预聚类
这确定了3个重复序列，现在已将其与以前的唯一序列合并。接下来要做的是利用pre.cluster命令对序列进行预聚类，允许序列间最多有2个差异，使序列进一步去噪。这个命令将对序列进行分组拆分，然后按丰度进行选择，从最大丰度到最小丰度，并识别相互之间2 nt以内的序列。如果他们被确定，则他们就会合并。一般接受每100 bp序列有1个差异：
mothur>pre.cluster(fasta=stability.trim.contigs.good.unique.good.filter.unique.fasta, count=stability.trim.contigs.good.unique.good.filter.count_table, diffs=2)
4 去除嵌合体及“不良品”
mothur>chimera.vsearch(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.count_table, dereplicate=t)
使用计数文件运行chimera.vsearch将从计数文件中删除嵌合体序列。但是仍然需要从fasta文件中删除这些序列。可以使用remove.seqs来做到这一点：
mothur>remove.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta,accnos=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.accnos)
运行summary.seqs，会看到剩下的内容：
mothur > summary.seqs(fasta=current, count=current)
请注意，序列从128656减少到118171，减少了8.1％；这是可接受的嵌合体数量。作为最后的质量控制步骤，需要确认数据集中是否存在任何“不良品”。有时，选择一个引物组去扩增到达此点的其他物质，例如来自古细菌、叶绿体和线粒体的18S rRNA基因片段或16S rRNA。也希望去除一些随机的东西。现在您可能会说：“等等，我想要那些东西”。但是，使用的引物只能扩增细菌的成员，如果它们击中了真核生物或古细菌，那么这是一个错误。要意识到叶绿体和线粒体在微生物群落中没有任何功能作用。接下来继续使用带有classify.seqs命令的贝叶斯分类器对这些序列进行分类：
mothur>classify.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.count_table,reference=trainset9_032012.pds.fasta, taxonomy=trainset9_032012.pds.tax, cutoff=80)
现在，所有内容都已分类。
最后，要删除“不良序列”，使用remove.lineage命令执行此操作：
mothur>remove.lineage(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta,count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.vsearch.pick.count_table,taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pds.wang.taxonomy, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)
请注意，如果分类器无法将序列分类为一个域，则“unknown”仅作为分类弹出。另外，请记住，如果没有使用RDP参考分类法对序列进行分类，还需要自定义谱系的名称。例如，修改后的RDP版本称为线粒体“线粒体”。如果使用greengenes，它将它们称为“f_线粒体”。如果正在使用greengenes（或SILVA或其他），则需要适当更改这些名称。
如果运行summary.seqs，将看到现在有2469个唯一序列，总共118009个序列。这意味着大约有162个序列位于这些不同的组中。现在，要创建反映这些删除的更新的分类法摘要文件，使用summary.tax命令：
mothur > summary.tax(taxonomy=current, count=current)
这将创建一个pick.tax.summary文件，其中删除了“不良序列”。至此，已尽可能地整理了数据，接下来就是要看看错误率是多少。
这篇推文对你有帮助吗？喜欢这篇文章吗？喜欢就不要错过呀，关注本知乎号查看更多的环境微生物生信分析相关文章。亦可以用微信扫描下方二维码关注“环微分析”微信公众号，小编在里面载入了更加完善的学习资料供广大生信分析研究者爱好者参考学习，也希望读者们发现错误后予以指出，小编愿与诸君共同进步！！！
学习环境微生物分析，关注“环微分析”公众号，持续更新，开源免费，敬请关注！
转载自原创文章：
Mothur3_处理改进的序列_免费生信分析课持续更新，干货满满
最后，再次感谢你阅读本篇文章，真心希望对你有所帮助。感谢！