其他
微生物组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物种
链接：
著作权归原作者所有。