Mothur2进阶_Mothur扩增子基因序列_数据预处理

本人在读研究生,方向环境微生物。之前在学习生物信息分析过程中在网络上四处奔走获取相关学习资料与解决问题,好生麻烦。于是,我就把与同学一起做的一些生物信息分析相关教程与经验总结搬运到这个CSDN这个大平台上来,希望能够与大家一起学习讨论。班门弄斧,大神见文多指教,抱拳抱拳抱拳抱拳!

本节重点讲解使用Mothur软件对扩增子序列进行拼接,切除与唯一化等预处理流程。

01数据的拼接

使用make.contigs完成每一个样品两端测序的拼接以及合并所有样本的数据。以stability.files作为输入文件,首先从fastq文件中提取序列和质量得分数据即生成单独的fasta和qual文件,创建反向读取的反向互补,然后连接序列形成contigs,最后在屏幕上会显示每个样本中的序列数。

命令注释:make.contigs用于读取一个正向fastq文件和一个反向fastq文件并输出一个新的fasta文件和report文件。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出文件:

stability.trim.contigs.fasta

stability.scrap.contigs.fasta

stability.contigs.report

文件注释
Name:序列名称;Length:长度;Overlap_length:重叠长度;Overlap_Start:重叠部分开始;Overlap_End:重叠部分结束;MisMatches:不匹配;Num_NsN:碱基数量;Expected_Error:预期错误率

输出文件:stability.contigs.groups

文件注释:第一列为序列名称,第二列为样本名称。

接下来使用summary.seqs命令查看序列,

命令注释:summary.seqs命令用于总结(概述)一个未排序或已排序的fasta格式文件的序列质量。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

屏幕输出结果显示:有152360个序列,大部分在248和253个碱基之间变化。数据集中最长的读取时间是502 bp。回想一下,每个读取应该是251 bp。说明这些片段显然没有很好地(或根本没有)组装起来。另外,序列中至少有2.5%不明确的碱基被识别。

输出文件:stability.trim.contigs.summary

文件注释 
Seqname:序列标识符; Start:起始位置;End:终止位置;Nbases:总碱基数,可以看做序列长度 ; ambigs:ambiguous bases 模糊碱基的数目;Polymer: homopolymer 碱基的最大长度;NumSeqs: 序列数,对于每条序列来说,其值总是为1

02数据的切除

运行screen.seqs删除所有碱基不明确和长度超过275 bp的序列,

命令注释:该命令用于保留符合用户自定义标准的序列,也可用于从names、group、contigsreport、alignreport、summary文件中挑选出不符合用户自定义标准的序列。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出文件:

stability.trim.contigs.good.fasta

stability.trim.contigs.bad.accnos

文件注释:第一列序列名称,第二列ambig。

输出文件:stability.contigs.good.groups

文件注释:第一列序列名称,第二列组名。

还有另一种方法可以使用summary.seqs的输出来运行,屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

该命令会使结果输出速度更快,因为summary.seqs输出文件已经为序列计算了模糊碱基的数量和序列长度。另外,Mothur记住了make.contigs使用了8个处理器,在当前的程序中也使用8个处理器。

运行get.current查看Mothur还记忆了什么,

命令注释:get.current命令用于显示数据处理过程中目前不同类型的最新版文件。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

从Mothur保存的当前文档可以看出Mothur记住了最新的fasta文件和group文件以及所拥有的处理器数量。

输出文件:current_files.summary

文件注释:由于文件名很长,对于许多命令使用“current”选项通常是最简单的。为了展示如何具体地使用Mothur,所以下面选择性地使用current选项。

运行命令:mothur > summary.seqs(fasta=stability.trim.contigs.good.fasta)

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出结果如上所示,测序错误率可能下降了一个数量级以上并且有128872个序列。

输出文件:stability.trim.contigs.good.summary

名词解释

Processors:电脑处理器线程

Gigabytes:十亿字节

Fastq文件:是一种包含质量值的序列文件,q为quality,一般用来存储原始测序数据,扩展名为fastq或者fq。目前Illumina测序,BGISEQ,Ion Torrent,pacbio,nanopore都以fastq格式存储测序数据,其中Illumina,BGISEQ一般是双末端测序,一般是一对命名为*_R1.fq.gz与*_R2.fq.gz的文件。

Fasta文件:无质量信息的序列,主要用于存储生物的序列文件,例如基因组,基因的核酸序列以及氨基酸等,是最常见的生物序列格式,一般扩展名为fa,fasta,fna等。该文件中,第一行是由大于号“>”开头的任意文字说明,用于序列标记(为了保证后续分析软件能够区分每条序列,单个序列的标识必须是唯一的,序列ID部分可以包含注释信息)。从第二行开始为序列本身,只允许使用既定的核苷酸或氨基酸编码符号。序列部分可以在一行也可以分成多行,序列中允许空格、换行、空行,直到下一个大于号,表示该序列的结束。

03数据的唯一化

数据中的许多序列是彼此重复的,所以使用unique.seqs命令来唯一化序列。

命令注释:unique.seqs命令只会输出fasta格式测序文件中的唯一化序列和一个表征它们与参考序列一致的文件。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出结果注释:若识别出两个序列具有相同序列,会被视为重复序列并合并。屏幕输出中有两列,第一列是特征序列的数量,第二列是剩余的唯一序列的数量。序列从128872变为16426。

输出文件:

stability.trim.contigs.good.names

stability.trim.contigs.good.unique.fasta

接下来运行count.seqs生成一个表,行是唯一序列的名称,列是组的名称。然后,用每个唯一序列在每个组中出现的次数填充表格。

命令注释:count.seqs又作make.table。该命令用于计算name文件中代表性序列的数量。如果提供了group文件,也可用于提供group中的组别数量。

屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出文件:stable.trim.contigs.good.count_tables

文件注释:第一列代表序列名称,第二列代表序列总数,其余每列可以看出每组代表序列个数。

通过使用count选项来使用这个文件,屏幕输出:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

输出文件:stability.trim.contigs.good.unique.summary

名词解释

重复序列:在基因组中,出现两次及以上的序列。

唯一化序列:将重复序列合并。

本文提供所有输出文件,百度网盘下载链接:

https://pan.baidu.com/s/11jMceGJYBFftLod6fqGEjw

提取码:1234

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

Mothur2进阶_Mothur扩增子基因序列_数据预处理

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

转载自原创文章:

Mothur2进阶_Mothur扩增子基因序列_数据预处理

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

上一篇:【Bioinfo Blog 005】【Python Code 001】——FASTA文件处理(未完)


下一篇:模拟一个简单计算器_阅读模拟器的简单介绍