微生物组16S rRNA数据分析小结:从fastq测序数据到OTU table

笔记内容:
由二代测序产生的序列数据(fastq格式)到物种丰度的OTU table, 样本群落距离矩阵,物种多样性指数,序列的进化树及物种注释信息的分析结果,为常规分析流程。可以使用usearch, vsearch, qiime等分析软件实现,在必要的时候需要根据样本信息的具体情况编写脚本予以实现。以下为双端测序(pair-end)产生的fastq格式数据常规分析流程,包括:

  1. fastqc及multiqc
  2. (qiime1)multiple_join_paired_ends.py
  3. fastq_screen
  4. (qiime1)multiple_split_libraries_fastq.py
  5. cutadapt
  6. usearch 流程

需要先了解barcode, primer的概念,可以通过下图快速了解一下:barcode是每个样本对应的一段与样本一一对应的短序列,像超市货物条形码一样区别每个样本。Fwd和Rev primer(引物)是用于16S扩增与建库的引物序列,在测序技术中起到定位的作用。有些下机数据会已经将样本按照barcode分好,每个样本测得的序列片段都放在一个与样本对应的文件里。

input: fastq格式的下机数据(pair-end):
已经按照各样本的barcode分好,格式为每个样本两个序列文件,比方说sampleID_1.fastq.gz, sampleID_2.fastq.gz。1为正向测序一次得到的结果,2为反向测序一次得到的结果。这是pair-end双向测序的原理,所以每个样本得到两个序列文件。
(注:标注了(qiime1)表明需要在qiime环境中运行。)

  1. fastqcmultiqc
    可以先把raw data,即下机数据做一遍fastqc以了解其质量情况。在安装好fastqc, multiQC之后,在终端输入:
    $ fastqc -c dir/*.fastq.gz -o fastqout/
    即在fastqout/文件夹中得到所有样本的测序质量报告。可以将这些报告整合成一个大的report,为.html格式,双击在浏览器内打开即可。
    $ multiqc fastqout/
  2. (qiime1)multiple_join_paired_ends.py
    把正反两次测序结果join在一起,正常的运行结果应为每个样本一个文件夹,文件夹内为三个文件,一个是该样本join到一起的序列文件(fastqjoin.join.fastq),另外两个为1,2两个序列中无法join到一起的序列文件(fastqjoin.un1.fastq, fastqjoin.un2.fastq)。

multiple_join_paired_end.py为join_paired_ends.py的多个input文件用法,在qiime环境下的linux终端中使用。可以设置各种参数调整,比如常用--read1_indicator--read1_indicator用于根据文件名称定位哪个fastq文件为read1和read2。default设置为'_R1_''_R2_'。其他参数可详见官网。
介于运行结果为下图左边所示,每个sample对应一个文件夹,储存了上述三个文件,且文件名没有改变。这一步可以用python写个脚本用于文件名的调整,去掉没有join到一起去的两个文件,及统计每个样本join起来的reads数占总reads数的比例,以了解双端测序join的情况。

  1. fastq_screen
    fastq_screen是可以用于检查是否存在污染序列的软件,需要先install再在终端中使用。FastQ Screen Documentation非常详细。在终端中使用也不难,input为上一步multiple_join后,修改了文件名的.fastq文件。
    (需要将包含了污染序列的tag的文件去除,留下screen之后的文件,自己写脚本实现)

  2. (qiime1)multiple_split_libraries_fastq.py
    multiple_split_libraries的目的为把上一步得到的一系列,每个样本的fastq文件,整合成一个大fastq文件,且每行开头以“@样本ID”标注,如下所示。红色虚线框里是sample1的第一个序列片段,第一行为样本ID及第几条序列,以及一些其他的,后续可以去除的信息,接下来就是序列数据。加号单独占一行,接下来为一些质量信息,以上也为fastq文件的格式的一些说明。

如下图所示,样本ID整合到fastq文件中的结果:

  1. cutadapt
    使用cutadapt去掉序列头尾的引物(primer)
    测序公司会把头尾两段引物序列给你。cutadapt有两个用于指定这段引物序列的参数-a-g。你可以先在终端用grepless查看一下手头序列信息有没有去掉引物序列。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等字母即可。

  2. usearch流程

fastq_filter
这个命令是一个序列质量控制的过程,使用-fastq_maxee 来设置每个reads所引起的expected errors。比方说设置成0.5,则会删除expected errors大于0.5的reads.

bmp-Qiime2Uparse.pl
这是BMP的脚本,该脚本的目的为将上一步fastq_filter得到的fasta文件header重新格式化,以便于后续的分析。

fastx_uniques
fastx_uniques是一个查找重复序列的过程,将重复序列归拢在一起,从而只显示unique序列(里面同时包含了每条unique序列出现的次数信息)。
-fastout 则output为fasta文件,且以各unique序列丰度降序排列。
-sizeout则在output文件的header中加上其重复reads的数目。

cluster_otus
这一步将上述得到的unique进行聚类,按照97%相似度聚为各OTU。并输出每个OTU的代表序列,以及每个OTU包含的序列数目。
Pick OTU一般有三种策略,详见16s微生物组第十条,这里usearch用的是UPARSE算法,本质上应该是一种denove的策略。

这一步还去除了chimera(嵌合体序列)。 chimera是因为在扩增中出现了问题,导致两个来自不同DNA模板的序列被扩增为同一条序列。可以在一个chimera reference数据库中找到相关序列的chimera并去除。


以下为一个从fastq_filter到cluster_otus的过程中,样本与序列信息的结果小结:

otutab
这一步将根据前面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进行过滤。

otutab_trim
设置参数用于去除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

otutab_norm
指定参数-sample_size,将每个样本的reads数目抽平到指定的reads数目。比如 -sample_size 5000即把每个样本的reads总数目抽平到5000.这一步不涉及对样本或者OTU的删除或规整,只是归一化。

-cluster_agg
以OTU代表序列的fasta格式文件为Input,生成otu.tree(OTU的进化树)。

-alpha_div
以trim及norm后的OTU table为input, 生成alpha diversity的各种多样性指数。

-beta_div
以otu.tree及OTU table作为input, 生成beta diversity中的各种距离矩阵。注意想要生成Unifrac距离矩阵必须用指定-tree参数(tree file)。

-sintax
这一步用于根据各OTU的代表序列进行物种预测,并给每个OTU添加物种注释。以OTU代表序列的fasta文件及物种参考数据库为input。参考的物种数据库一般为一个包含详细物种信息的fasta格式的文件。usearch官网上推荐使用一个小而权威的参考数据库,可以参考RDP或SILVA LTP等16S数据库。
这一步得到的output为以下格式的文档:

OTU_ID tax_with_condifence strand tax
OTU1 d:Bacteria(1.0000),p:Firmicutes(0.9300),c:Clostridia(0.8400),o:Clostridiales(0.8400),f:Ruminococcaceae(0.3200),g:Clostridium_III(0.2400) - d:Bacteria,p:Firmicutes,c:Clostridia,o:Clostridiales

第二列中括号里为该物种注释的置信度,第三列为strand, 即正负链。一般为全是正链(+)或全是负链(-)。第四列为根据-sintax_cutoff参数设置界值,比如设置为0.8则只纳入第二列中置信度大于0.8的物种注释。

经过上述步骤,得到trim及标准化后的OTU table, otu.tree, alpha多样性指数文件,beta diversity各距离矩阵及物种注释信息sintax.txt。后续分析参考微生物组16S rRNA数据分析小结:从OTU table到marker物种

附:本流程中使用version10的usearch, qiime的版本参数如下(在qiime1环境中使用print_qiime_config.py查看版本信息)

推荐阅读更多精彩内容