微生物组16S rRNA数据分析小结：从fastq测序数据到OT

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

1. (qiime1)multiple_join_paired_end.py
2. fastq_screen
3. (qiime1)multiple_split_libraries_fastq.py
4. cutadapt
5. 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. (qiime1)multiple_join_paired_end.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的情况。

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

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

   

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

4. cutadapt
   使用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等字母即可。

5. 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包含的序列数目。
这一步还去除了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物种