16s分析之Qiime聚类OTU

首先Qiime官网命令都集中在这个网址，大家可以收藏，随用随查：

http://qiime.org/scripts/index.html

这次我们关注的是聚类OTU，与其相关的命令如下：

这里主要有三个封装的命令，也就是三种方式聚类otu：

de novo：无参考数据库，所有序列都进行聚类处理。

closed-reference：与参考数据库比对，留下比对上的序列，丢弃比对不上的序列。

open-reference：与参考数据库比对，比对不上的序列进行de novo聚类

因此，这三种聚类方式严格程度排序：

closed-reference > open-reference > de novo

本次我采用pick_open_reference_otus.py，详细参数参考（http://qiime.org/scripts/ pick_open_reference_otus.html）

准备数据，接之前去除嵌合体结果继续跑：

参数文件常调的参数我列出来，大家可以针对性去查看：

assign_taxonomy:reference_seqs_fp

assign_taxonomy:id_to_taxonomy_fp

assign_taxonomy:confidence #仅用于RDP、Mothur方法，

assign_taxonomy:rdp_max_memory #最大内存分配，对于java虚拟机使用RDP方法时。

pick_otus:max_rejects

pick_otus:stepwords

pick_otus:max_accepts

pick_otus:word_length

pick_otus:enable_rev_strand_match True/F#将序列反向并聚类out，会耗费双倍的时间；

align_seqs:template_fp

align_seqs:alignment_method

align_seqs:pairwise_alignment_method

需要注意的是Qiime使用的参考数据是gg_13_8，所以之后如果需要做功能预测时，要修改成gg_13_5

这次我编写参数文件（存为tab分隔的txt文件）：

pick_otus:enable_rev_strand_match True#是否反向匹配

pick_rep_set:rep_set_picking_method most_abundant#挑选最大丰度的序列为代表序列（默认不是这样）

pick_rep_set:log_fp pick_rep.log#输出log文件（默认不输出）

make_phylogeny:log_fp make_phylogeny.log#输出log文件（默认不输出）

#本次我采用open_reference方法聚类：

pick_open_reference_otus.py -o otus/ -i seqs_chimeras_filtered.fna-p uc_fast_params.txt

参数含义：

-o：输出文件夹

-i：输入文件（为fa格式注意：后缀名为fasta/fna/fa都是可以的，都一样）

-p：是参数文件

结果就像这样：

结果出来了，对于新手而言最终要的就是这些文件是什么？

final_otu_map.txt：就是一列otu编号，一列含有这个otu的样品，如下图：

大家发现这里面有的otu只有一个样品有，所以就有了下面的这个文件；

final_otu_map_mc2.txt：去除了singleton；

new_refseqs.fna：实际上就是参考数据库复制出来的文件

otu_table_mc2.biom：不含有物种注释的otu.table文件，当然是biom格式的；

otu_table_mc2_w_tax.biom：有物种注释的otu.table文件；

otu_table_mc2_w_tax_no_pynast_failures.biom：对齐后的otu.table文件，做进化树用的；

rep_set.fna：代表序列文件

rep_set.tre：进化树文件

uclust_assigned_taxonomy/rep_set_tax_assignments.txt：物种注释文件；

高能的文件来了：

log_20180106052656.txt：这个文件包含这次聚类全部过程和你的参数文件还有Qiime的配置信息，非常值得解读：

Qiime配置信息：

参数文件也在这里：

这是一条封装的命令，包含许多条命令，从聚类到注释到做进化树，因此当我们可以摘抄出来，分开跑：

下面分析open聚类的整个过程（这是之前我跑别的数据的参数文件，当时分析的，为了省事，我直接复制过来了）：

# Pick Reference OTUs command这步根据参考数据库聚类otu，与pick_close_reference_otus一样；使用的是qiime默认的数据库；13.8版本的GREENGENE数据库，可以更换数据库。

pick_otus.py -i seqs_chimeras_filtered.fna -o otus//step1_otus -r/usr/local/lib/python2.7/dist-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fasta-m uclust_ref --enable_rev_strand_match --suppress_new_clusters

参数解释

-o otus//step1_otus #输出文件夹

-m uclust_ref picking OTUs方法，sortmerna, mothur,trie, uclust_ref, usearch, usearch_ref, blast, usearch61, usearch61_ref,sumaclust, swarm, prefix_suffix, cdhit, uclust，有以上这几种，uclust是默认的，

--enable_rev_strand_match反向比对，会使占用内存加倍

--suppress_new_clusters当序列与参考数据库不匹配时，不聚新的otu

# Generate full failures fasta file command从源文件中挑选出来第一次聚otu失败的序列，

filter_fasta.py -f seqs_chimeras_filtered.fna -s otus//step1_otus/seqs_chimeras_filtered_failures.txt-o otus//step1_otus/failures.fasta

参数解释：

-s，这个文件上一步生成的，里面就一列，样品名称后面跟序列编号，代表在参考数据库找不到这条序列

-n这个命令加上，是去除txt中序列，并输出剩余序列的

# Pick rep set command得到挑选好的代表序列文件，

pick_rep_set.py -i otus//step1_otus/seqs_chimeras_filtered_otus.txt-o otus//step1_otus/step1_rep_set.fna -f seqs_chimeras_filtered.fna

参数解释：

-i是mapping文件，上一步生成的，里面就是第一排序列编号，第二排拥有这条序列的样品名称

-o输出第一步得到的代表序列

-f从源文件中挑选出来的

备注：第一次聚otu完成都生成两个txt文件：

seqs_chimeras_filtered_failures.txt：就是那个样品中哪条序列没有参考比对上

seqs_chimeras_filtered_otus.txt：就是源聚otu文件中哪条序列参考比对上的，后面还跟有这条序列存在于哪个样品的哪条序列。就是第一步聚out成功的

抽取失败的文件使用API

# Subsample the failures fasta file usingAPI

python -c "import qiime;qiime.util.subsample_fasta('/home/qiime/Desktop/usearch/otus/step1_otus/failures.fasta','/home/qiime/Desktop/usearch/otus/step2_otus/subsampled_failures.fasta','0.001000')"Forcing --suppress_new_clusters as this is reference-based OTUpicking.

将原始序列第一步聚类失败的挑选出来，作为新的文件

# Pick de novo OTUs for new clusterscommand从第一步失败的子文件中再次聚otu

pick_otus.py -i otus//step2_otus//subsampled_failures.fasta -ootus//step2_otus/ -m uclust --denovo_otu_id_prefix New.ReferenceOTU --enable_rev_strand_match

-m uclust调用方法

--denovo_otu_id_prefix New.ReferenceOTU采用denovo方法，新的OTU命名前面加上New.ReferenceOTU

--enable_rev_strand_match反向比对，默认不进行此项，会使占用内存加倍

# Pick representative set for subsampledfailures command子序列合成的otu中，在subsampled_failures.fasta序列中挑选代表序列

pick_rep_set.py -i otus//step2_otus//subsampled_failures_otus.txt -ootus//step2_otus//step2_rep_set.fna -fotus//step2_otus//subsampled_failures.fasta

# Pick reference OTUs using de novo rep setcommand从参考数据库比对失败的序列文件之前就建好了，从中以第二步聚的otu得到的代表序列为参考序列在第一步失败文件中聚otu，再次得到otu

pick_otus.py -i otus//step1_otus/failures.fasta -o otus//step3_otus/-r otus//step2_otus//step2_rep_set.fna -m uclust_ref --enable_rev_strand_match--suppress_new_clusters

# Create fasta file of step3 failurescommand将第三步聚otu失败的序列根据其failures_failures.txt文件过滤出来，得到failures_failures.fasta，也就是再聚失败文件。

filter_fasta.py -f otus//step1_otus/failures.fasta -sotus//step3_otus//failures_failures.txt -o otus//step3_otus//failures_failures.fasta

# Pick de novo OTUs on step3 failurescommand将两次失败的序列进行--denovo_otu_id_prefix方法的聚otu

注意就是New.CleanUp.ReferenceOTU这次otu命名添加了CleanUp，这个新的名字，代表这是第三次聚otu的结果。

pick_otus.py -i otus//step3_otus//failures_failures.fasta -ootus//step4_otus/ -m uclust --denovo_otu_id_prefix New.CleanUp.ReferenceOTU --enable_rev_strand_match

# Merge OTU maps command合并otu

cat otus//step1_otus/seqs_chimeras_filtered_otus.txtotus//step3_otus//failures_otus.txtotus//step4_otus//failures_failures_otus.txt > otus//final_otu_map.txt

# Pick representative set for subsampledfailures command

pick_rep_set.py -i otus//step4_otus//failures_failures_otus.txt -ootus//step4_otus//step4_rep_set.fna -fotus//step3_otus//failures_failures.fasta

# Filter singletons from the otu map usingAPI去掉单个的序列，就是这条序列在所有样品中仅仅只有一个样品才有的序列

python -c "import qiime;qiime.filter.filter_otus_from_otu_map('/home/qiime/Desktop/usearch/otus/final_otu_map.txt','/home/qiime/Desktop/usearch/otus/final_otu_map_mc2.txt', '2')"

# Write non-singleton otus representativesequences from step1 to the final rep set file: otus//rep_set.fna

# Copy the full input refseqs file to thenew refseq file将完整的参考序列文件复制为一个新的文件new_refseqs.fna

cp/usr/local/lib/python2.7/dist-packages/qiime_default_reference/gg_13_8_otus/rep_set/97_otus.fastaotus//new_refseqs.fna

# Write non-singleton otus representativesequences from step 2 and step 4 to the final representative set and the newreference set (otus//rep_set.fna and otus//new_refseqs.fna respectively)

将第二步和第四步得到的无单序列的代表序列导入到rep_set.fna和完整的参考序列文件复制为一个新的文件new_refseqs.fna这个文件中。

这一步结果应该是rep_set.fna这个文件中包含第一，第二，第四步的代表序列# Make the otu table command

make_otu_table.py -i otus//final_otu_map_mc2.txt -o otus//otu_table_mc2.biom

基于去除单个序列的txt文件，制作out.table

# Assign taxonomy command注释，将物种信息加上

assign_taxonomy.py -o otus//uclust_assigned_taxonomy -iotus//rep_set.fna

# Add taxa to OTU table command这一步将注释加到otu_table中，并命名为otu_table_mc2_w_tax.biom

biom add-metadata -i otus//otu_table_mc2.biom--observation-metadata-fpotus//uclust_assigned_taxonomy/rep_set_tax_assignments.txt -ootus//otu_table_mc2_w_tax.biom --sc-separated taxonomy --observation-headerOTUID,taxonomy

下面的序列都是为了做树文件的前处理

# Align sequences command对齐序列

align_seqs.py -i otus//rep_set.fna -o otus//pynast_aligned_seqs

# Filter alignment command过滤代表序列文件，并存为rep_set_aligned_pfiltered.fasta

filter_alignment.py -o otus//pynast_aligned_seqs -iotus//pynast_aligned_seqs/rep_set_aligned.fasta

# Build phylogenetic tree command构建系统发育树

make_phylogeny.py -iotus//pynast_aligned_seqs/rep_set_aligned_pfiltered.fasta -o otus//rep_set.tre

这份参数文件有我不会的地方，比如API的使用？相信大家有能力吃掉它

#分门类统计相对丰度：

summarize_taxa_through_plots.py -o taxa_summary -i otus/otu_table_mc2_w_tax.biom-m map.txt

结果就像这样，当然这个子文件夹中还有图，直接打开浏览器就可以查看：

#其实我们可以用一个更简单的命令去跑，不出图，只有txt文件，之后做LEfSe分析也方便一些：

summarize_taxa.py -i otus/otu_table_mc2_w_tax.biom -osummarize_taxa_L6/ -m map.txt --delimiter ''

--delimiter ''#分隔符我用''也就是LEfSe分析要求的，默认是“;”

#采用如此命令总结，就是简单将序列统计信息输出到屏幕上简单看一下跑完的数据的序列统计

biom summarize-table -iotus/otu_table_mc2_w_tax.biom

#如果大家想输出统计文件：

biom summarize-table -i otu_table.biom -ootu_table.sum

#下面对otu_table进行抽平d值设置成最小值（序列数最小的样品的序列数），关于这个值如何设置，最小的，还是去掉最小序列的样品，用大部分样品有的序列数进行抽平，看上面统计结果自己安排，我用最小序列数抽平即可

single_rarefaction.py -i otus/otu_table_mc2_w_tax.biom -ootu_table_even6073.biom -d 6073

参数解释：

Øi, --input_path

Ø-o, --output_path

Ø-d, --depth抽样深度

到此，之后就是转换一些格式和做一些分析了，用Qiime作分析，用R出图是我默认的方式，谢谢大家查看