其他
微生物组16S rRNA数据分析小结:从fastq测序数据到OTU table
推荐阅读
fastqc及multiqc (qiime1)multiple_join_paired_ends.py fastq_screen (qiime1)multiple_split_libraries_fastq.py cutadapt usearch 流程
已经按照各样本的barcode分好,格式为每个样本两个序列文件,比方说sampleID_1.fastq.gz, sampleID_2.fastq.gz。1为正向测序一次得到的结果,2为反向测序一次得到的结果。这是pair-end双向测序的原理,所以每个样本得到两个序列文件。
(注:标注了(qiime1)表明需要在qiime环境中运行。)
可以先把raw data,即下机数据做一遍fastqc以了解其质量情况。在安装好fastqc, multiQC之后,在终端输入:
$ fastqc -c dir/*.fastq.gz -o fastqout/
即在fastqout/文件夹中得到所有样本的测序质量报告。可以将这些报告整合成一个大的report,为.html格式,双击在浏览器内打开即可。
$ multiqc fastqout/
把正反两次测序结果join在一起,正常的运行结果应为每个样本一个文件夹,文件夹内为三个文件,一个是该样本join到一起的序列文件(fastqjoin.join.fastq),另外两个为1,2两个序列中无法join到一起的序列文件(fastqjoin.un1.fastq, fastqjoin.un2.fastq)。
--read1_indicator
和--read1_indicator
用于根据文件名称定位哪个fastq文件为read1和read2。default设置为'_R1_'
和'_R2_'
。其他参数可详见官网。介于运行结果为下图左边所示,每个sample对应一个文件夹,储存了上述三个文件,且文件名没有改变。这一步可以用python写个脚本用于文件名的调整,去掉没有join到一起去的两个文件,及统计每个样本join起来的reads数占总reads数的比例,以了解双端测序join的情况。
fastq_screen是可以用于检查是否存在污染序列的软件,需要先install再在终端中使用。FastQ Screen Documentation非常详细。在终端中使用也不难,input为上一步multiple_join后,修改了文件名的.fastq文件。
(需要将包含了污染序列的tag的文件去除,留下screen之后的文件,自己写脚本实现)
multiple_split_libraries的目的为把上一步得到的一系列,每个样本的fastq文件,整合成一个大fastq文件,且每行开头以“@样本ID”标注,如下所示。红色虚线框里是sample1的第一个序列片段,第一行为样本ID及第几条序列,以及一些其他的,后续可以去除的信息,接下来就是序列数据。加号单独占一行,接下来为一些质量信息,以上也为fastq文件的格式的一些说明。
使用cutadapt去掉序列头尾的引物(primer)
测序公司会把头尾两段引物序列给你。cutadapt有两个用于指定这段引物序列的参数
-a
和-g
。你可以先在终端用grep
和less
查看一下手头序列信息有没有去掉引物序列。grep的简单用法可以参考答生信从业人员的N个问题(更新非常慢): 乱入一个grep怎么用grep -n 'ATGCATGC' seqs.fastq
把含有primer的行显示出来grep -c 'ATGCATGC' seqs.fastq
显示有多少行有primer如果显示出来大量序列在头尾均包含primer序列,那说明primer应该没被去掉。同理你可以在去掉之后这样看看,到底有没有去掉。在
grep -n
的时候你可以留意一下primer是出现在头部还是尾部。全部出现在头部则把这段序列指定在-g
之后,在尾部则指定在-a
之后。注意这里不能反了。这两个参数是用于决定把这段序列之前还是之后的序列切掉的。比方说把尾部的序列填充到了-g
, 则把尾部之前的序列都切掉,那output就没有东西了。其他cutadapt参数为质量控制的参数,比如设置一个长度阈值,在切掉primer后如果序列长度小于这个值则扔掉这个序列等。详细可见官网文件。
另外primer中可能使用了简并碱基,比如M,R,K,Y等这样不是ATGC四个密码子的值。在终端grep的时候需要使用正则把他们对应的ATGC写出来,但是在cutadapt中不用,直接填充M,R等字母即可。
这个命令是一个序列质量控制的过程,使用-fastq_maxee 来设置每个reads所引起的expected errors。比方说设置成0.5,则会删除expected errors大于0.5的reads.
这是BMP的脚本,该脚本的目的为将上一步fastq_filter得到的fasta文件header重新格式化,以便于后续的分析。
fastx_uniques是一个查找重复序列的过程,将重复序列归拢在一起,从而只显示unique序列(里面同时包含了每条unique序列出现的次数信息)。
-fastout 则output为fasta文件,且以各unique序列丰度降序排列。
-sizeout则在output文件的header中加上其重复reads的数目。
这一步将上述得到的unique进行聚类,按照97%相似度聚为各OTU。并输出每个OTU的代表序列,以及每个OTU包含的序列数目。
Pick OTU一般有三种策略,详见16s微生物组第十条,这里usearch用的是UPARSE算法,本质上应该是一种denove的策略。
这一步将根据前面qiime中得到的seqs.fastq,及OTU代表序列,将各个样本的序列与OTU代表序列相比对,生成OTU table。
-sample_delim 是一个指定样本ID和序列ID分隔符的参数,比方说qiime生成的seq.fastq是以“sample1_0”,“sample1_1”命名的,_ 之前的内容才是样本名称。则可以把sample_delim指定为 _。
-otutabout
这一步得到了最初的OTU table,即每个sample对应在每个OTU中有多少个reads。接下来可以查看一下样本中OTU的丰度如何,都分布在什么范围内。
biom summarize-table -i otu_raw.tab
随后可以用otutab_trim来对raw OTU table进行过滤。
设置参数用于去除OTU table中丰度太低的样本或者OTU。
-min_sample_size 5000
即去除reads数小于5000的样本
-min_otu_size 2
去除reads数小于2的OTU
-min_otu_freq 0.001
去除 单个OTUreads数/总OTUreads数目 小于0.001的OTU
指定参数-sample_size,将每个样本的reads数目抽平到指定的reads数目。比如 -sample_size 5000即把每个样本的reads总数目抽平到5000.这一步不涉及对样本或者OTU的删除或规整,只是归一化。
以OTU代表序列的fasta格式文件为Input,生成otu.tree(OTU的进化树)。
以trim及norm后的OTU table为input, 生成alpha diversity的各种多样性指数。
以otu.tree及OTU table作为input, 生成beta diversity中的各种距离矩阵。注意想要生成Unifrac距离矩阵必须用指定-tree参数(tree file)。
这一步用于根据各OTU的代表序列进行物种预测,并给每个OTU添加物种注释。以OTU代表序列的fasta文件及物种参考数据库为input。参考的物种数据库一般为一个包含详细物种信息的fasta格式的文件。usearch官网上推荐使用一个小而权威的参考数据库,可以参考RDP或SILVA LTP等16S数据库。
这一步得到的output为以下格式的文档:
经过上述步骤,得到trim及标准化后的OTU table, otu.tree, alpha多样性指数文件,beta diversity各距离矩阵及物种注释信息sintax.txt。后续分析参考微生物组16S rRNA数据分析小结:从OTU table到marker物种
链接:
著作权归原作者所有。