微生物多样性qiime2分析流程(2）使用qiime2和dada2处理16S序列

我不太确定这篇文档里还有没有打错的命令，代码建议参考：
一文完不成qiime2微生物多样性分析(2021版上)这篇文档，应该没错误命令了,希望对大家有帮助，喜欢的小伙伴可以关注我的公众号R语言数据分析指南持续分享更多优质资源

1. 构建样本文件清单 sample.txt (注：以制表符分割)

sample-id   forward-absolute-filepath   reverse-absolute-filepath
A1  $PWD/data/A1_16S_R1.fastq   $PWD/data/A1_16S_R2.fastq
A2  $PWD/data/A2_16S_R1.fastq   $PWD/data/A2_16S_R2.fastq
A3  $PWD/data/A3_16S_R1.fastq   $PWD/data/A3_16S_R2.fastq

2. 导入文件的`QIIME2`命令

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

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

2.2 去除引物

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

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

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

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

4. `dada2` 序列质量控制和构建特征表

截断位点可以通过figaro软件得到dada2插件正反向截断位点信息
参考：https://www.jianshu.com/p/66ebd1d4558a
也可以设置--p-trunc-len-f 0，--p-trunc-len-r 0 表示不进行修剪，

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

dada2插件的输出包括ASV表，代表序列以及有关该过程的一些统计信息，所有内容均为压缩格式(此过程特别耗费时间）测试41个样本耗费4小时

4.1 导出特征表

mkdir phyloseq
qiime tools export \
--input-path table.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
cd ../

在这里导出otu_table.tsv是因为我们需要对ASV ID进行修改；

image.png

dada2输出的ASV ID如图所示，这显然不是我们想看到的，因此我们需要将其改为ASV1,ASV2.....等；可以通过以下R代码进行修改：

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

rown <- paste0("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 = F)

数据导出的时候我们是按如下步骤进行：

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

现在我们更改了otu_table.tsv内容后再按照如下步骤将其还原回去：

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

biom convert -i otu_table.tsv -o feature-table.biom --to-hdf5 --table-type="OTU table"
qiime tools import \
  --input-path feature-table.biom \
  --type 'FeatureTable[Frequency]' \
  --input-format BIOMV210Format \
  --output-path otu_table.qza

4.2 导出代表序列

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

4.2.1 更改代表序列名称：

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)

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.qza

4.3 特征表统计

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

4.4 代表序列统计

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

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

time qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences rep-seqs.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

注：构建进化树耗费20分钟
通过以上步骤我们得到了ASV表，代表序列及进化树，接下了可以根据参考数据库训练特质分类器进行物种注释，最终可通过phyloseq与MicrobiotaProcess等R包进行数据可视化。