一文完不成qiime2微生物多样性分析(2021版上)

之前写过一系列qiim2的教程都是零零散散的记录，里面也有不少命令打错的步骤给各位小伙伴造成了不少的疑惑，恰逢qiime2-2020.11发布，所以就在重新总结一次，这次稍后会整理上传分析数据，希望对大家有所帮助，喜欢请关注个人公众号R语言数据分析指南获取，在此先行拜谢了

conda activate qiime2 
# conda deactivate

1. 更改样本名称

ls *.1.fq|awk -F "." '{print "mv "$_" ""_16s_R1.fastq"}' >> R1.sh
ls *.2.fq|awk -F "." '{print "mv "$_" ""_16s_R2.fastq"}' >> R2.sh
sh R1.sh;sh R2.sh

2. 去除barcode(根据自身情况设置参数)

mkdir trim
for i in *_16s_R2.fastq;do echo seqtk trimfq -b 7 -e 7 $i echo trim/$i;done|sed 's/echo/>/g' > trim2.sh
for i in *_16s_R1.fastq;do echo seqtk trimfq -b 7 -e 7 $i echo trim/$i;done|sed 's/echo/>/g' > trim1.sh
sh trim2.sh;sh trim1.sh

3. 构建样本文件清单(sample.txt)

less $PWD/trim/*fastq|perl -e 'while(<>){chomp;@a=split/\//;@b=split/\_16s/,$a[-1];$inf=<>;print"$b[0]\t$_\t$inf"}'|sed '1i sample-id\tforward-absolute-filepath\treverse-absolute-filepath' > sample.txt

4.1 导入数据

time qiime tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-path sample.txt \
--output-path paired-demux.qza \
--input-format PairedEndFastqManifestPhred33V2

注：如果导入的原始数据已经去掉引物则直接执行下一步，
如果未去除引物则执行以下命令去除引物（替换为自己的引物序列）

4.2 去除引物

time qiime cutadapt trim-paired \
--i-demultiplexed-sequences paired-demux.qza --p-front-f GTGCCAGCMGCCGCGG \
--p-front-r CCGTCAATTCMTTTRAGTTT \
--o-trimmed-sequences paired-end-demux.qza \
--verbose \
&> primer_trimming.log

4.3.创建可视化文件查看质量

time qiime demux summarize --i-data paired-end-demux.qza --o-visualization demux.qzv

.qzv文件 可以通过https://view.qiime2.org/查看

5 figaro获得dada2截断位点

dada2插件对原始数据paired-end-demux.qza进行质量过滤需要两个参数，trunc-len-f和trunc-len-r。这个想法是通过尽可能多地删除较低质量部分来优化正向和反向的读取与合并，并且仍然保留足够的重叠部分。我们可以通过检查原始数据的质量值来获得过滤参数。现在可以通过figaro程序来得到合适的参数.

5.1 安装figaro (https://github.com/Zymo-Research/figaro)

wget http://john-quensen.com/wp-content/uploads/2020/03/figaro.yml
conda env create -n figaro --file figaro.yml

5.2. figaro的使用方法

wget https://github.com/Zymo-Research/figaro/archive/master.zip
unzip master.zip 
mv figaro-master figaro
cd figaro
chmod 777 figaro.py

5.3. 运行figaro

conda activate figaro 
python $PWD/figaro/figaro.py -i trim -o trim -f 16 -r 20 -a 300 -F zymo

输出文件示例

trimParameters.json
forwardExpectedError.png
reverseExpectedError.png
要获取推荐的截断参数，请查看trimParameters.json:如下所示
[
    {
        "trimPosition": [
            141,
            181
        ],
        "maxExpectedError": [
            2,
            3
        ],
        "readRetentionPercent": 87.24,
        "score": 82.24249607692565
    },
  ]
推荐的正向截断位置为141，
推荐的反向截断位置为181。
修剪和截断后，正向读取的预期错误数为2，
反向读取的预期错误数为3,
qiime2 dada2插件合并87.24％的读数

• -i 输入文件文件夹
• -o 结果输出文件夹
• -f 正向引物的长度。如果已除去引物，请输入1。
• -r 反向引物的长度。如果已除去引物，请输入1。
• -a 预期的合并扩增子长度(即引物扩增长度），可以保守一些，并提供稍大的价值。(此处以515F-806R为示例）
• -F 文件名格式(请严格按照代码所示填写）

6.1 dada2序列质量控制和构建特征表

time qiime dada2 denoise-paired \
--i-demultiplexed-seqs paired-end-demux.qza \
--p-trunc-len-f 141 \
--p-trunc-len-r 181 \
--o-table table.qza \
--o-representative-sequences rep-seqs.qza \
--o-denoising-stats denoising-stats.qza

6.2 vsearch对特征表进行聚类

qiime vsearch cluster-features-de-novo \
  --i-table table.qza \
  --i-sequences rep-seqs.qza \
  --p-perc-identity 0.97 \
  --o-clustered-table table-dn-97.qza \
  --o-clustered-sequences rep-seqs-dn-97.qza

注此处需要使用vsearch对特征表进行聚类原因请参考:
https://toolshed.g2.bx.psu.edu/repository/display_tool?repository_id=a43ea4908816827c&tool_config=%2Fsrv%2Ftoolshed%2Fmain%2Fvar%2Fdata%2Frepos%2F003%2Frepo_3930%2Fqiime2%2Fqiime_vsearch_cluster-features-open-reference.xml&changeset_revision=71f124e02000&render_repository_actions_for=tool_shed

https://forum.qiime2.org/t/clustering-denovo-specifying-input-sequence-processing/16306
https://docs.qiime2.org/2020.11/plugins/available/vsearch/cluster-features-de-novo/

6.3 导出特征表

mkdir phyloseq
qiime tools export --input-path table-dn-97.qza --output-path phyloseq

biom convert -i phyloseq/feature-table.biom -o phyloseq/otu_table.tsv --to-tsv
cd phyloseq; sed -i '1d' otu_table.tsv
sed -i 's/#OTU ID/ASV/' otu_table.tsv

6.4 R更改特征表名称

library(pacman)
pacman::p_load(tidyverse,magrittr,stringr)
otu <- "otu_table.tsv" %>%
read.delim(check.names = FALSE,header = T,sep="\t")

rown <- patse0("ASV",seq_len(nrow(otu)))
otu[,1] <- rown
colnames(otu)[1] <- paste0("ASV",colnames(data)[1])
write.table (otu,file ="otu_table.tsv",sep ="\t",row.names = T,quote = F)

# otu_table.tsv --> feature-table.biom ---> otu_table.qza

将格式转化为qza

biom convert -i otu_table.tsv -o feature-table.biom --to-hdf5 --table-type="OTU table"

time qiime tools import \
  --input-path feature-table.biom \
  --type 'FeatureTable[Frequency]' \
  --input-format BIOMV210Format \
  --output-path otu_table-rename.qza

6.5 导出代表序列

qiime tools export --input-path rep-seqs-dn-97.qza --output-path phyloseq

6.6 更改代表序列名称

less dna-sequences.fasta |paste - -|sed '1i ASVID,seq' > rep.fa

library(pacman)
pacman::p_load(tidyverse,magrittr,stringr)
rep <- "rep.fa" %>%
read.delim(check.names = FALSE, row.names = 1) %>%
set_rownames(paste0(">ASV", seq_len(nrow(.))))
write.table (rep,file ="rep.xls", sep ="\t", row.names = T,quote = F)

less rep.xls|sed '1d'|sed 's/"//g'|sed 's/\r//g'|tr "\t" "\n" > rep-seqs.fasta

将代表序列转换成qza格式

time qiime tools import
--type 'FeatureData[Sequence]'
--input-path rep-seqs.fasta
--output-path rep-seqs-rename.qza

6.7 特征表统计

time qiime feature-table summarize
--i-table table-rename.qza
--o-visualization table.qzv
--m-sample-metadata-file sample-metadata.tsv

6.8 代表序列统计

time qiime feature-table tabulate-seqs
--i-data rep-seqs-rename.qza
--o-visualization raw.fq.list

7.比对代表性序列，并构建系统发育树

time qiime phylogeny align-to-tree-mafft-fasttree \
 --i-sequences rep-seqs-rename.qza \
 --o-alignment aligned-rep-seqs.qza \
 --o-masked-alignment masked-aligned-rep-seqs.qza \
 --o-tree unrooted-tree.qza \
 --o-rooted-tree rooted-tree.qza

导出无根进化树

mkdir phyloseq

qiime tools export \
--input-path unrooted-tree.qza \
--output-path phyloseq
cd phyloseq; mv tree.nwk unrooted_tree.nwk; cd ..

导出有根进化树

qiime tools export \
--input-path rooted-tree.qza \
--output-path phyloseq
cd phyloseq; mv tree.nwk rooted_tree.nwk;cd .. 

8.训练特征分类器

8.1.导入参考序列数据库

time qiime tools import \
--type 'FeatureData[Sequence]' \
--input-path silva.16s-bacteria.fasta \
--output-path silva.16s-bacteria.qza

8.2.导入物种分类注释数据库

time qiime tools import \
--type 'FeatureData[Taxonomy]' \
--input-format HeaderlessTSVTaxonomyFormat \
--input-path silva.16s-bacteria.tax \
--output-path silva.16s-bacteria.tax.qza

8.3. 训练分类器

time qiime feature-classifier fit-classifier-naive-bayes \
--i-reference-reads silva.16s-bacteria.qza \
--i-reference-taxonomy silva.16s-bacteria.tax.qza \
--o-classifier classifier.qza

注：此步骤特别耗费时间

8.4 数据注释方法(1)

这一步输入我们得到的代表序列文件，对其进行分类注释

time qiime feature-classifier classify-sklearn \
--i-classifier classifier.qza \
--i-reads rep-seqs-rename.qza \
--o-classification taxonomy.qza

导出注释文件

qiime tools export \
--input-path taxonomy.qza \
--output-path phyloseq

注:注释步骤也特别耗费时间

8.5 数据注释方法(2)

此处也可以运行qiime1中的assign_taxonomy.py脚本进行数据注释
因此需要安装qiime1,或者下载我上传的此脚本使用时调用即可

assign_taxonomy.py -m uclust \
-i rep-seqs.fasta \
-r silva.16s-bacteria.fasta  \
-t silva.16s-bacteria.tax -o taxonomy > tax.log

-m 可选uclust、blast 、rdp
-e 默认为0.001，仅适用于BLAST法

assign_taxonomy.py -m blast -e 0.01 \
-i rep-seqs.fasta \
-r silva.16s-bacteria.fasta  \
-t silva.16s-bacteria.tax -o taxonomy > tax.log

具体参考:http://qiime.org/scripts/assign_taxonomy.html

qiime1的优秀脚本清单：http://qiime.org/scripts/index.html

对注释结果进行排序

less dna-sequences_tax_assignments.txt|sort -k1.4n|\
awk 'BEGIN{OFS=FS="\t"}{print , }' > otu_taxa.xls

经测试此方法大大提高了注释速度

9. 整合分析结果

此处有2种方式，一种是以.qza的压缩格式整理数据，另一类则是以导出的文本格式整理

9.1 .qza法

rm(list = ls())
pacman::p_load(tidyverse,MicrobiotaProcess,phyloseq)
otu <- "otu_table.qza"
rep <- "rep-seqs.qza"
tree <- "rooted-tree.qza"
tax <- "taxonomy.qza"
sample <- "group.txt"
dada2 <- import_qiime2(otuqza=otu, taxaqza=tax,refseqqza=rep,
                          mapfilename=sample,treeqza=tree)

9.2 整合文本文件

将otu表，注释信息表，进化树，样本信息整合，即以下文件
由于qiime2进行数据注释较慢所以才采用此种方法

otu_table.tsv
group.xls
rooted_tree.tre
otu_taxa.xls

对otu_taxa.xls进行数据格式处理

head otu_taxa.xls

ASV1    Unassigned
ASV2    Unassigned
ASV3    Unassigned
ASV4    d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Burkholderiales;f__Nitrosomonadaceae;g__Ellin6067;s__uncultured bacterium
ASV5    d__Bacteria;p__Verrucomicrobiota;c__Verrucomicrobiae;o__Chthoniobacterales;f__Chthoniobacteraceae;g__Candidatus Udaeobacter

less otu_taxa.xls|sed '1i ASV,domain,phylumt,class,order,family,genus,species'|sed 's/,/\t/g'|sed 's/;/\t/g' |sed 's/[a-z]__//g' > taxa.xls 

修改为如下格式

head taxa.xls

ASV domain  phylumt class   order   family  genus   species
ASV1    Unassigned
ASV2    Unassigned
ASV3    Unassigned
ASV4    Bacteria    Proteobacteria  Gammaproteobacteria Burkholderiales Nitrosomonadaceae   Ellin6067   uncultured bacterium
ASV5    Bacteria    Verrucomicrobiota   Verrucomicrobiae    Chthoniobacterales  Chthoniobacteraceae Candidatus Udaeobacter
ASV6    Unassigned

phyloseq包进行数据整合

经过上面的一系列处理我们得到了最终文件，接下来通过phyloseq包进行数据整合

安装phyloseq包

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("phyloseq")

rm(list=ls())
pacman::p_load(tidyverse,phyloseq,MicrobiotaProcess,ape)
otu_mat <- read.delim2("otu_table.tsv",header=T,
                       sep="\t",check.names = F,row.names = 1) %>%
  as.matrix()

tax_mat <- read.delim("taxa.xls",header=T,row.names = 1,
                      sep="\t",check.names = F) %>% as.matrix()
samples_df <- read.delim("group.xls",header = T,row.names = 1,
                         sep="\t",check.names = F)
tree <- read.tree("rooted_tree.tre")

OTU = otu_table(otu_mat,taxa_are_rows =T)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)

ps <- phyloseq(OTU,TAX,samples,tree)
ps

注：此处需注意group.xls文件中样本名称不能是纯数字

> ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3761 taxa and 12 samples ]
sample_data() Sample Data:       [ 12 samples by 1 sample variables ]
tax_table()   Taxonomy Table:    [ 3761 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 3761 tips and 3760 internal nodes ]