2020-09-15

Mothur的使用

1.数据的准备

使用mothur对分样之后的数据进行分析，每个样品的正向序列和反向序列各一个fastq文件，如mothur软件官网给出的测试数据：

2.设置数据输入输出文件夹

make.file(inputdir=MiSep_SOP,type=fastq, prefix=stability)；

此时，MiSep_SOP文件夹里必须有相对应的fastq文件，否则无法设置数据输入输出文件

3. #序列拼接（把正反向序列合并），减少测序和pcr错误#
make.contigs(file=stability.files, processors=8)

这一步是一个样品一个样品先后进行的，

全部样品拼接完成之后代码行会出现以下结果

显示每个样品的序列条数和全部样品的总序列数

4. 对序列拼接的结果进行总结

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

运行这条命令之后，会看到如下结果

显示序列大小及其分布情况；根据这个结果设置下一步的序列过滤参数

5. 对数据结果进行过滤

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

或者运行这一条命令

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

 

运行结束后我们来看看mothur知道些什么

get.current()

Current files

 saved by mothur:

fasta=stability.trim.contigs.good.fasta 

group=stability.contigs.good.groups

processors=8

 

6.挑选unique序列，因为总的数据量会非常大，挑选unique序列是一步去冗余的过程

unique.seqs(fasta=stability.trim.contigs.good.fasta)

 

7. 然后对得到的unique序列进行简单分析

 count.seqs(name=stability.trim.contigs.good.names, group=stability.contigs.good.groups)

这条领命运行结束之后得到结果文件 stability.trim.contigs.good.count_table

8. 对上一条命令得到的结果进行统计分析

summary.seqs(count=stability.trim.contigs.good.count_table)

9. 接下来，我们需要制作一个参考数据库

必须知道在比对从哪里开始、哪里结束。将keepdots设置为false，去除前后的点。

因为16S扩增子测序一般得到的序列都比较短也就300-400bp，所以我们必须要把从网上下载下来的接近全长的序列截到跟我们自己的序列位点相同的区域生成一个新的参比数据库，比对结果才有可能比较准确。

pcr.seqs(fasta=silva.bacteria.fasta, start=11894, end=25319, keepdots=F, processors=8)

这里的start和end都是根据测序区域来确定的，测试数据的V4区的序列；如果是V3V4或者其他区域，需要自行查找起始和结束位点。

运行结束后对制作的参比数据库进行重命名

 rename.file(input=silva.bacteria.pcr.fasta, new=silva.v4.fasta)

10. 对数据库进行简单统计

summary.seqs(fasta=silva.v4.fasta)

11. 进行序列比对

align.seqs(fasta=stability.trim.contigs.good.unique.fasta, reference=silva.v4.fasta)

12. 对得到的结果进行简单统计

summary.seqs(fasta=stability.trim.contigs.good.unique.align, count=stability.trim.contigs.good.count_table)

13. 根据summary的结果再一次对序列进行过滤

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)
summary.seqs(fasta=current, count=current)

14. 把不符合的序列过滤掉

 filter.seqs(fasta=stability.trim.contigs.good.unique.good.align, vertical=T, trump=.)

结束后会得到如下结果

15. 在修剪末端的时候，我们可能创造了一些冗余，所以运行unique.seqs

unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta, count=stability.trim.contigs.good.good.count_table)

16. 接下来进行预聚类，即去噪（允许出现2个错误）。这条命令将会对序列分组，按照丰度从高到低排列，识别出有2个差异的序列。如果序列被合并，我们倾向于允许100bp里有一个差异：

pre.cluster(fasta=stability.trim.contigs.good.unique.good.filter.unique.fasta, count=stability.trim.contigs.good.unique.good.filter.count_table, diffs=2)

17. 去嵌合

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)

18. 从fasta文件中取出嵌合体

 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)

19. 对上一条命令得到的结果进行简单统计

summary.seqs(fasta=current, count=current)

20. 用贝叶斯分类器检测数据集里是否还含有不符合条件的序列

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)