Mothur配合qiime生成otutable

qiime设置傻瓜化，安装起来欲仙欲死，操作起来爽翻天。第一次用qiime就发觉这个软件果然酸爽。

公司测序回来发了好几个文件回来，我选择质量控制并拼接好的数据直接上手。

P.S.文件暂缺，后面会补上

质控并拼接好的样品为：trim.2.fq, trim.4.fq, trim.7.fq, trim.11.fq, trim.14.fq, trim.18.fq, trim.19.fq, trim.23.fq, trim.25.fq, trim.28.fq, trim.31.fq, trim.35.fq, trim.37.fq, trim.41.fq, trim.44.fq，合计15个样品。样品如果使用qiime拆分的话，序列如下（以2号样为例）：

@2_241533    M02419:322:000000000-BDP8B:1:2119:24295:14257 1:N:0:CAAGAC    orig_bc=CTGTA    new_bc=CTGTA    bc_diffs=0
TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGTGCGTAGGTGGCTAGGTAAGTCAGATGTGAAAACCCAGGGCTTAACTTTGGGACTGCATTTGAAACTGCTTGGCTAGAGTGCAGGAGAGGTAAGTGGAATTCCTAGTGTAGCGGTGAAATGCGTAGATATTAGGAGGAACACCAGTGGCGAAGGCGGCTTACTGGACTGTAACTGACACTGAGGCACGAAAGCGTGGGGAGCAAACAGGATTAGAAACCCTTGTAGTCC
+
GGGGGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGEGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFE;FGGGGGGGGGGFGGGGGGGGGGGGGFGGGG
...

序列命名必须以@样品名_数字开头。下面开始正式处理：

1. 文件合并

> cat 文件1 文件2 ... > 合并文件

cat trim.2.fq trim.4.fq trim.7.fq trim.11.fq trim.14.fq trim.18.fq trim.19.fq trim.23.fq trim.25.fq trim.28.fq trim.31.fq trim.35.fq trim.37.fq trim.41.fq trim.44.fq > seq.fq

也可以利用通配符一下搞定

cat *.fq > seq.fq

P.S.图形操作界面：选中文件夹中所有文件，右键选择属性，复制所文件名到新文本中，替换“,”。

2.利用Mothur拆分fq文件

mothur > fastq.info(fastq=seq.fq)

3.打开qiime软件

export PATH=~/miniconda2/bin:$PATH #设定conda路径，不知道为啥我这个每次用都要重设一下
source activate qiime1

打开后显示如下

(qiime1) ruminant@ruminant-PowerEdge-T630[qiimeOTU] 

 

4. 看一下序列分布

 

count_seqs.py -i seq.fasta

结果显示为：

551360  : seq.fasta (Sequence lengths (mean +/- std): 272.8898 +/- 0.5639)
551360  : Total

长度基本分布在273，和我们预期长度一致。P.S.看序列分布和长度还是用Mothur比较爽，后面讲到Mothur再说。

5.OTUpicking

先准备一个配置文件：params.txt，正反向比对都打开，因为测序的可能是负链也可能是正链。文件内容如下：

pick_otus:enable_rev_strand_match True

用官方推荐的open_reference的方法做聚类，默认的reference是greengene13.8（老掉牙了，看到Mothur和Usearch在可劲的喷该数据库）,在/usr/share/qiime/data/gg_13_8_otus/rep_set。-r修改参考基因，-m修改聚类方法。15个样的数据量大概需要20分钟，简直是飞快。

 

pick_open_reference_otus.py -o otus1/ -i seq.fasta -p params.txt -r /usr/share/qiime/data/gg_13_8_otus/rep_set/97_otus.fasta

生成biom格式otu表2个:otu_table_mc2.biom, otu_table_mc2_w_tax.biom，前一个没有注释，后一个有注释

 

6.转换biom文件

 

biom格式不是给人看的，简直累死。我们用biom程序转化他为我们熟悉的二维表

 

cd otus1
biom convert -i otu_table_mc2.biom -o otutable.txt --to-tsv
biom convert -i otu_table_mc2_w_tax.biom -o otutable_tax.txt --to-tsv #注释结果没在表里
biom convert -i otu_table_mc2_w_tax.biom -o otutable_tax.txt --to-tsv --header-key taxonomy #注释结果在表里
 

代表序列在文件rep_set.fna，树文件在rep_set.tre。

可选步骤. 如果没有注释出来，可以用Mothur注释OTU

mothur> classify.seqs(fasta=rep_set.fna,  template=/home/ruminant/Documents/database/tax&seqs/greengene/gg_13_8_99.fasta, taxonomy=/home/ruminant/Documents/database/tax&seqs/greengene/gg_13_8_99.gg.tax)

在此生成了一个rep_set.gg.wang.taxonomy的文件，内容就是OTU和对应的物种。本次有5000多个OTU，我决定用libreoffice calc来处理，还是可视化的界面比较爽。OTU表格中新建一列Taxonomy，Q3中输入=VLOOKUP(A3,[rep_set.gg.wang.taxonomy]rep_set.gg.wang!$A:$B,2,0)，填充手柄填充，大功告成。OTU表长这个样子：