Mothur—16s双端测序处理

老板不知道从哪弄了一堆肠道微生物的16s高通量数据丢给我,说了一通要求,然后让我分析。而我向来不善于拒绝,秉着“不会,硬着头皮也要弄出来”的原则,开始了回家前一段奇妙的mothur摸索之旅。

Question:

之前都是用qiime来处理数据,根据barcode来识别序列,给序列编号,同时去掉低质量的片段;而且不拼接。但这批来自苏州某测序公司的数据,既没有给我们barcode也没有给primer。所以有三个问题对我来说是新的:
1,去低质量的片段;
2,双端的拼接;
3,给每条序列编号。

言归正传,还是看看如何假借mothur之手,同时结合qiime来分析这批数据吧。

一.mothur的安装

bio-linux系统自带,这一块没有做,所以也就不提了,以后如果安装了,再来补一块的东西。
下载 http://www.mothur.org/wiki/Download_mothur

二.运行方式

三. 分析流程

分析大全:http://www.mothur.org/wiki/Analysis_examples
下面我们还是以MiSeqSOP (standard operating procedure)为例: http://www.mothur.org/wiki/MiSeq_SOP,以下内容都是这个的翻译,如果不嫌麻烦,还是尽量看英文的。
为了尊重Mothur工作组的劳动,如果用了他家的东西,感觉还不错的话,记得引上他们的文章 Kozich JJ,Westcott SL, Baxter NT, Highlander SK, Schloss PD. (2013): Development of adual-index sequencing strategy and curation pipeline for analyzing ampliconsequence data on the MiSeq Illumina sequencing platform. Applied andEnvironmental Microbiology. 79(17):5112-20。

 

四、实例

1,解压缩

2, trimmomatic过滤双端序列

3 paire end根据overlap合并

下面我们就来完成所有样品两端测序的拼接,通过make.contigs命令来实现,这个命令需要stability.files作为输入文件。他的原理在于:对一对序列进行比对后分析共有的那个区域,先看其中两条序列不匹配的地方,如果一条序列的某个碱基对应另一条序列的gap,这个碱基的质量必须高于25才算真实的;如果两条序列在哪个位置都一个碱基,我们需要其中一个条序列上该碱基质量大于6或者得分大于另一条上的;如果得分小于6,则认为该碱基为N。Processors为线程,我的电脑其实是12的,哈哈。

 查看合并后的结果

4. 对contigs的过滤

 

5. get.current()查看当前的有效对象:

#提升效率

6.unique序列

7.pcr.seqs 优化数据库(http://www.mothur.org/wiki/Pcr.seqs)

8. align.seqs序列比对(http://www.mothur.org/wiki/Align.seqs)

下面这一步是干什么,没明白?(http://www.mothur.org/wiki/Filter.seqs)
mothur > filter.seqs(fasta=stability.trim.contigs.good.unique.good.align, vertical=T, trump=.)

因为上一步已经去掉了首尾,这个时候需要重新在unique一下
mothur > unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta, count=stability.trim.contigs.good.good.count_table)

9.pre.cluster(http://www.mothur.org/wiki/Pre.cluster)

10.去嵌合体

11.remove.lineage(http://www.mothur.org/wiki/Remove.lineage)

12 评估错误(这一步我没弄,需要用他们的数据集)

13  聚OTUs

这个命令可以生成a .list, .rabund, .sabund files三个文件 Nearest neighbor: Each of the sequences within an OTU are at most X% distant from the most similar sequence in the OTU. Furthest neighbor: All of the sequences within an OTU are at most X% distant from all of the other sequences within the OTU. Average neighbor: This method is a middle ground between the other two algorithms. 默认的是Average neighbor 在聚类之前,先是把文件切割成几部分,有3种方法可以用来切割文件:通过距离文件、分类地位、分类地位加上fasta文件。splitmethod对应的3种方法:distance、classify、fasta。默认的为distance。

为什么上面的有的距离丢失了呢
Missing distances Perhaps the second most commonly asked question is why there isn’t a line for distance 0.XX. If you notice the previous example the distances jump from 0.003 to 0.006. Where are 0.004 and 0.005? mothur only outputs data if the clustering has been updated for a distance. So if you don’t have data at your favorite distance, that means that nothing changed between the previous distance and the next one. Therefore if you want OTU data for a distance of 0.005 in this case, you would use the data from 0.003.原来如果矩阵相同,那么label就合并了啊,所以看到的第一列这个label是断开的。 Next we want to know how many sequences are in each OTU from each group and we can do this using the make.shared command.

上面的几个策略都试了一下,生成的数据过于庞大(10T),没办法,就只有采用qiime来接着分析了。qiime在聚OTU的时候选取cd-hit算法能够减少硬盘使用量

ABNR1_1_1101_26846_9215 G-G-G-A-G-GC-A-GC-A-G-T-G-G-G-G-A-A–TA-TTGC-A-C-AA-T-G-G–GG–GA-A-A-C-CC-TG-A-TGCAGC-GACGCCGCGT-G-G-A-G–GA-A-G-A-A–G-G-TC———-TT-CG—————GA-TTG-T-A–AA-C-T-CC—–TG-TT-G-T–T-GGG—-G–A-A–G—A————————T———–AA——————————-T-G-A-C-G—–G-T-A-C-CC—A-A-C-A-A-G—G–AA–GT-G–AC-G–GC-TAA-C–T-A-C-G-T-G-CCA-G-C–A-GC-CG-C—-GG–TA-AA–AC-GT-AG-GTC–ACA-A-G-C-GT—T-GT-C-CGG-AA-TT-A-C-T–GG-GT-GT–A—AA-GG-GA-GC—–G-CA-G-G-C-G–G–G-AA-G-A-C-AAG-TT-G-G-A-A-G–TG–A-AA-TC-T-A-TG-G-G–CT-C-AA-C-C-C-A-T-A-A-A-C-T-G-C-T-T-T-C–AA-A-ACT-G-T–TT–T-T-C-T-T-GA-GT-A-G-TG—-CA-G-A–G-G-T-A–GG-C—-GG-A-ATT-C-C-C-G-GTGT-A-GCGGT-G-G-AA-TGCGTAGATA-TC–G-GG-A-G-G-A-ACA-CC-AG-T-GGC-GAA-G–G-C–G-G–C-C-T-A—CTG-G–GC-A-C-C-A-A-CT-GA-CG-CT-G-A-GG-CT-CG-AAA-G-T-G-TG-GG-T-AG-C-AAA-CA–GG-AT-TA-G-ATA–C-C-C-T-G-GTA-G

下面就是进化分析,alpha,beta等分析。

总之,我们也可以通过mothur一键式来完成所有任务
$ ./mothur “#make.contigs(file=stability.files, processors=8); screen.seqs(fasta=current, maxambig=0, maxlength=275); unique.seqs(); count.seqs(name=current, group=current); align.seqs(fasta=current, reference=silva.v4.fasta); screen.seqs(fasta=current, count=current, start=1968, end=11550, maxhomop=8); filter.seqs(fasta=current, vertical=T, trump=.); pre.cluster(fasta=current, count=current, diffs=2); unique.seqs(fasta=current, count=current); pre.cluster(fasta=current, count=current, diffs=2); chimera.uchime(fasta=current, count=current, dereplicate=t); remove.seqs(fasta=current, accnos=current); classify.seqs(fasta=current, count=current, reference=trainset9_032012.pds.fasta, taxonomy=trainset9_032012.pds.tax, cutoff=80); remove.lineage(fasta=current, count=current, taxonomy=current, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota); remove.groups(count=current, fasta=current, taxonomy=current, groups=Mock); cluster.split(fasta=current, count=current, taxonomy=current, splitmethod=classify, taxlevel=4, cutoff=0.15); make.shared(list=current, count=current, label=0.03); classify.otu(list=current, count=current, taxonomy=current, label=0.03); phylotype(taxonomy=current); make.shared(list=current, count=current, label=1); classify.otu(list=current, count=current, taxonomy=current, label=1);”

参考资料:

一个有趣的博客

《Mothur—16s双端测序处理》有14个想法

  1. 博主,你好!非常感谢你的经验分享!我是生物信息学零开始的菜鸟,最近在捣鼓16srDNA分析工具QIIME,在做beta_diversity分析执行beta_diversity_through_plots.py命令运行时提示RuntimeWarning错误提示/usr/localb/python2.7/dist-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:107: RuntimeWarning: The result contains negative eigenvalues. Please compare their magnitude with the magnitude of some of the largest positive eigenvalues. If the negative ones are smaller, it’s probably safe to ignore them, but if they are large in magnitude, the results won’t be useful. See the Notes section for more details. The smallest eigenvalue is -0.0189198266496 and the largest is 0.227502762093.
    不知道博主之前有没遇到过这种问题?该如何解决?谢谢!

    1. 具体怎么解决的,我不清楚,但是建议你可以去qiime官网(需要翻墙)上咨询,回复很快的。

  2. 我的输入:qiime@qiime-190-virtual-box:~/mywork/second$beta_diversity_through_plots.py -i ./usearch_qf_results_80/otu_table.biom -obdiv_even100/ -m ./samples_Map.txt -e 8349 -t./phylogeny/rep_set_align/rep_set.tre
    报错提示如下:
    /usr/local/lib/python2.7/dist-packages/skbio/stats/ordination/_principal_coordinate_analysis.py:107:RuntimeWarning: The result contains negative eigenvalues. Please compare theirmagnitude with the magnitude of some of the largest positive eigenvalues. Ifthe negative ones are smaller, it’s probably safe to ignore them, but if theyare large in magnitude, the results won’t be useful. See the Notes section formore details. The smallest eigenvalue is -0.0187522128398 and the largest is0.226088912333.
    RuntimeWarning

  3. 楼主 我能问一下 没有barcode 和 Primer 最后怎么把序列归到对应的样本上吗?
    我菜鸟,啥都不懂的…还望不吝赐教呀!
    我最近做的和你的类似,需要用sra数据库上别人发表的数据进一步分析,但是没有barcode等信息,不知道如何进一步处理了

  4. 另外,我按照那个链接下载了silva 119文件,解压之后找不到fasta文件呀,求教博主~

  5. 楼主,你好,文章最后一段:(总之,我们也可以通过mothur一键式来完成所有任务),为什么我每次把所有的命令放入,提交后就出现问题,具体是怎么操作的,能不能详细地分解一下?谢谢

  6. 博主你好,我是一名研究生,最近有一批16s的数据想用QIIME处理,自己装他的教程装了好久,不知博主是否有流程可以让我学习一下?十分感谢。

  7. 博主你好,我有个问题想请教您一下:就是我们在做alpha多样性之前,需要对otu表根据最低测序深度进行随机重抽样,即拉平处理,我知道qiime中的parallel_multiple_rarefactions.py可以实现otu矩阵稀疏化,并会输出好几个稀疏化的文件,问了公司的人 说用R取均值则可得到最终的拉平后的otu矩阵,请问这个取“均值”是什么意思啊,应该怎么理解,是好几个文件中重复出现的么?q期待您的解答。

    1. 随机抽样处理的本质是要做normalization,保证两两数据之间可比较性。你可以网上查一下normalization的方法,方法很多,他们说的应该是其中的一种方式。

  8. 博主,您好!生信菜鸟,想用的您流程跑16s 双端测序的结果,但是,mothur运行的第一步batch files就没有运行成功……invalid command……希望得到您的帮助。谢谢!

发表评论

电子邮件地址不会被公开。 必填项已用*标注