推荐阅读

1.ggplot2绘制曼哈顿图示例2.phyloseq | 用 R 分析微生物组数据及可视化3.R语言PCA分析教程 | Principal Component Methods in R4.vegan包进行微生物群落主坐标分析(PCoA)及ggplot2作图5.R语言 | 相关性分析6.排序分析(Ordination Analysis)导读：由二代测序产生的序列数据(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.fastqc及multiqc
可以先把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的情况。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包含的序列数目。
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的OTUotutab_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为以下格式的文档：

第二列中括号里为该物种注释的置信度，第三列为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查看版本信息)
链接：https://www.jianshu.com/p/1928d2c87452
著作权归原作者所有。我就知道你在看?