基于qiime2的16S数据分析全流程：从导入数据到下游分析一条龙

目录

创建metadata

把数据导入qiime2

去除引物序列

双端合并 （dada2不需要）

质控 （dada2不需要）

使用deblur获得特征序列

使用dada2生成代表序列与特征表

物种鉴定

可视化物种鉴定结果

构建进化树（ITS一般不构建进化树）

生成多样性分析所需数据

α多样性分析

beta多样性分析

绘制稀释曲线

差异物种分析

如何把.qza格式文件导出？比如把特征表导出

演示单端数据导入qiime2

---

基于qiime2的16s测序数据分析分析

本文演示的是目前最长用的双端测序数据

激活qiime环境

alias activate_qiime='conda activate qiime2-amplicon-2024.2'

我已经把上面这句代码放在了.bashrc里面，因此激活qiime2环境只需要运行activate_qiime即可

按照下图策略分析

演示Illumina公司的16S rRNA基因区的V3到V4区

数据导入，需要首先准备好一个manifest.tsv文件，这个文件有三列，分别是 sample-id forward-absolute-filepath reverse-absolute-filepath，如图所示

生成manifest.tsv可以使用下面的代码，这个代码是my_script文件夹里面的create_manifest.sh

#!/bin/bash

# 设置目录路径

dir="/data3/sunjs/chenlina/69contorl47PDAC"

# 输出表头 echo -e "sample-id\tforward-absolute-filepath\treverse-absolute-filepath"

# 遍历_R1和_R2文件

for r1 in ${dir}/*_R1.fastq.gz; do # 检查R1文件是否存在

if [[ -f $r1 ]]; then

# 获取样本ID（去掉_R1.fastq.gz）

sample_id=$(basename "$r1" _R1.fastq.gz)

# 构建R2文件路径 r2="${dir}/${sample_id}_R2.fastq.gz"

# 检查R2文件是否存在 if [[ -f $r2 ]]; then

# 输出结果 echo -e "${sample_id}\t${r1}\t${r2}" fi fi done

这个脚本会在终端直接显示生成的内容，我们一般会把这些内容重定向到另一个文件去，所以通常运行代码

bash create_manifest.sh>manifest.tsv

这样就会在当前工作目录下生成一个manifest.tsv文件

创建metadata

还需要创建一个metadata，比如有这些列

不过我创建的meta.tsv如图所示，只有样本名和分组两列

把数据导入qiime2

export TMPDIR=/data3/sunjs/ cd ~/test/qiime2/

input_path=~/test/qiime2/manifest.tsv

qiime tools import \

--type 'SampleData[PairedEndSequencesWithQuality]' \

--input-path "${input_path}" \

--output-path "PDAC.qza" \

--input-format PairedEndFastqManifestPhred33V2

去除引物序列

v3到v4区域的引物序列是固定的，就是这两个参数指定的序列

cd /data3/sunjs/chenlina/69contorl_47PDAC_res

qiime cutadapt trim-paired \

--i-demultiplexed-sequences 69contorl_47PDAC.qza \

--p-cores 10 \

--p-no-indels \

--p-front-f ACTCCTACGGGAGGCAGCAG \

--p-front-r GGACTACHVGGGTWTCTAAT \

--o-trimmed-sequences primer-trimmed-demux.qza

双端的数据不需要拆分barcode，直接往下进行即可

双端合并 （dada2不需要）

如果接下来打算 使用deblur或OTU聚类方法，现在就合并序列。如果计划使用 dada2对序列进行去噪，则不要合并——dada2会在对每个序列进行去噪之后自动执行序列合并。

qiime vsearch merge-pairs \

--i-demultiplexed-seqs primer-trimmed-demux.qza \

--o-unmerged-sequences unmerged_demux-joined.qza \

--p-threads 16 \

--o-merged-sequences demux-joined.qza

质控 （dada2不需要）

cd /data3/sunjs/chenlina/test qiime quality-filter q-score \

--p-min-quality 20 \ --i-demux demux-joined.qza \

--o-filtered-sequences demux-joined-filtered.qza \

--o-filter-stats demux-joined-filter-stats.qza

使用deblur获得特征序列

qiime deblur denoise-16S \

--i-demultiplexed-seqs demux-joined-filtered.qza \

--p-trim-length 400 \ --p-sample-stats \

--p-jobs-to-start 2 \ --p-min-reads 1 \

--o-representative-sequences repset-seqs.qza \

--o-table feature-table.qza \

--o-stats deblur-stats.qza

或者使用data2得到特征表和代表序列，可以直接替换掉质控+使用 deblur获得特征序列两步，但是不管用deblur还是dada2都要完成导入数据，去引物

使用dada2之前需要先获得 --p-trim四个参数，把下面这段代码生成的数据拖到qiime2的可视化网站里面可以看到两张图

qiime demux summarize \

--i-data primer-trimmed-demux.qza \

--o-visualization demux-summary.qzv

使用dada2生成代表序列与特征表

#p-trim-left-f正向序列左边

#p-trim-left-r 反向去列左边

#p-trunc-len-f正向序列x轴的值

#p-trunc-len-r反向序列x轴的值

# --p-retain-all-samples True保留所有样本

qiime dada2 denoise-paired \

--i-demultiplexed-seqs primer-trimmed-demux.qza \

--p-n-threads 20 \

--p-trim-left-f 0 --p-trim-left-r 0 \

--p-trunc-len-f 270 --p-trunc-len-r 187 \

--p-retain-all-samples True \

-o-table dada2-table.qza \

--o-representative-sequences dada2-rep-seqs.qza \

--o-denoising-stats denoising-stats.qza

其中

#p-trim-left-f正向序列左边

#p-trim-left-r 反向去列左边

#p-trunc-len-f正向序列x轴的值

#p-trunc-len-r反向序列x轴的值

这四个参数的取值要通过看上面那两张图来得到

一般来讲会去掉质量中位数小于20的位点。

上面两种方法选一个就行，推荐dada2

物种鉴定

生成代表序列之后可以进行物种鉴定了

export TMPDIR=/data3/sunjs/

input_path=/data3/sunjs/chenlina/qiime2/PDAC

output_path=/data3/sunjs/chenlina/qiime2/PDAC

db_path=~/Metagenome/database/silva-138-99-nb-classifier.qza

cd "${input_path}"

#全长引物注释,使用的数据库日期是2024.2

#要求输入代表序列，也就是dada2结果生成的那个

#--p-n-jobs如果设置为 -1，那么程序会尝试使用所有可用的 CPU 核心来执行任务

#--p-n-jobs指定的是进程而非线程，和threads不是一回事

qiime feature-classifier classify-sklearn \

--i-classifier "${db_path}" \

--i-reads dada2-rep-seqs.qza \

--p-n-jobs 20 \

--o-classification "${output_path}/taxonomy.qza" #过滤污染 叶绿体 线粒体 qiime taxa filter-table \ --i-table dada2-table.qza \

--i-taxonomy taxonomy.qza \

--p-exclude mitochondria,chloroplast \

--p-include p__ \

--o-filtered-table "${output_path}/feature-table-filt-contam.qza"

#过滤稀有ASV qiime feature-table filter-features \

--i-table feature-table-filt-contam.qza \

--p-min-frequency 50 \ --p-min-samples 1 \

--o-filtered-table "${output_path}/feature-table-final.qza"

#过滤掉ASV对应的序列 qiime feature-table filter-seqs \

--i-data dada2-rep-seqs.qza \

--i-table feature-table-final.qza \

--o-filtered-data "${output_path}/repset-seqs-final.qza"

这里需要提前下载好对应的数据库，下载的数据库理论上来讲是越新越好，但是需要和你用的qiime2软件版本一致，比如qiime2是2024.2的，那数据库也要下载2024.2的，下载数据库可以从这个网址下载 QIIME 2 Resources ，运行这段代码之后就得到了物种鉴定表也就是 taxonomy.qza，同时上面的代码还进行了一些过滤，然后我们可以使用下面的代码来对物种鉴定的结果进行可视化

可视化物种鉴定结果

qiime metadata tabulate \

--m-input-file taxonomy.qza \

--o-visualization taxonomy.qzv

#绘制柱状图

qiime taxa barplot \

--i-table feature-table-final.qza \

--i-taxonomy taxonomy.qza \

--m-metadata-file meta.tsv \

--o-visualization taxa-bar-plots.qzv

其中taxonomy.qzv文件拖到网站里面长这样子

taxa-bar-plots.qzv长这样子，网站允许把柱状图的数据导出成csv格式，推荐导出之后用R来绘制

构建进化树（ITS一般不构建进化树）

#构建进化树，ITS一般不构建进化树

export TMPDIR=/data3/sunjs/ input_path=/data3/sunjs/chenlina/qiime2/PDAC output_path=/data3/sunjs/chenlina/qiime2/PDAC

cd "${input_path}"

qiime phylogeny align-to-tree-mafft-fasttree \

--i-sequences repset-seqs-final.qza \

--o-alignment aligned-repset-seqs.qza \

--p-n-threads 40 \

--o-masked-alignment masked-aligned-repset-seqs.qza \

--o-tree ${output_path}/unrooted-tree.qza \

--o-rooted-tree "${output_path}/rooted-tree.qza"

下面这段代码用于生成代表序列和特征表的qzv格式文件，在网址 view.qiime2.org可以可视化，需要先导出然后拖进去

cd /data3/sunjs/chenlina/84contorl_50PDAC_res

# 特征表和特征序列汇总

qiime feature-table summarize \

--i-table feature-table-final.qza \

--m-sample-metadata-file meta.tsv \

--o-visualization table.qzv qiime feature-table tabulate-seqs \

--i-data repset-seqs-final.qza \

--o-visualization rep-seqs.qzv

把table.qzv拖到网站里面是这样的

点击detail这个，拖到最下边，出现这样的界面

后续将选择53这个数量进行抽平，当然51也可以

生成多样性分析所需数据

#53是把上一步生成的table.qzv文件拖到qimme2网站查看来确定的

qiime diversity core-metrics-phylogenetic \

--i-phylogeny rooted-tree.qza \

--i-table feature-table-final.qza \

--p-sampling-depth 53 \

--m-metadata-file meta.tsv \

--output-dir ./diversity

生成一个名为diversity的目录，里面有一堆文件，其中有四个是qzv格式的，这四个qzv格式的可以直接拖到网站进行可视化，都是beta多样性的比如PCA图，PCoA图这种的

比如下面这样子

α多样性分析

cd /data3/sunjs/chenlina/84contorl_50PDAC_res/

#α多样性分析

qiime diversity alpha-group-significance \

--i-alpha-diversity diversity/faith_pd_vector.qza \

--m-metadata-file meta.tsv \

--o-visualization faith-group-significance.qzv

输入是上一步生成的faith_pd_vector.qza，生成文件如图所示

beta多样性分析

#beta多样性分析

cd /data3/sunjs/chenlina/67contorl_50PDAC_res

qiime diversity beta-group-significance \

--i-distance-matrix diversity/unweighted_unifrac_distance_matrix.qza \

--m-metadata-file meta.tsv \

--m-metadata-column group \

--p-pairwise \

--o-visualization unweighted-unifrac-group-significance.qzv

生成了一个文件，拖到网页里面是这样的，代表了其中某个分组与其他分组有没有差别，使用的方法是置换多因素方差分析

画图，实际上前面已经有过这个图了

#qiime diversity core-metrics-phylogenetic这一步已经生成的结果

cd /data3/sunjs/chenlina/69contorl_47PDAC_res

qiime emperor plot \

--i-pcoa diversity/unweighted_unifrac_pcoa_results.qza \

--m-metadata-file meta.tsv \ --p-custom-axes days-since-experiment-start \

--o-visualization unweighted-unifrac-emperor-group.qzv

这里--i输入的是前面

qiime diversity core-metrics-phylogenetic

这一步已经生成的结果，最后得到的结果也是跟前面

qiime diversity core-metrics-phylogenetic运行得到的那个结果是一样的，尚不清楚这个有啥用处。

绘制稀释曲线

#稀释曲线

cd /data3/sunjs/chenlina/69contorl_47PDAC_res

qiime diversity alpha-rarefaction \

--i-table feature-table-final.qza \

--i-phylogeny rooted-tree.qza \

--p-max-depth 9000 \

--m-metadata-file meta.tsv \

--o-visualization alpha-rarefaction.qzv

我这里选择了取样是9000条，结果如图所示

可以看到后面基本平了，说明测序饱和了，如果发现后面还在往上升，说明测序没饱和，很多东西测不到，这种数据发出去会被拒稿，因为说不定再多测点就结果不一样了。

差异物种分析

有时候可能会先过滤部分数据，比如我先按照频率和特征数过滤 掉一些低质量数据

#过滤部分特征表的数据

cd /data3/sunjs/chenlina/69contorl_47PDAC_res

#--p-min-frequency：样本必须具有的最小总频率，才能被保留

#--p-min-features：样本必须具有的最小特征数，才能被保留

qiime feature-table filter-samples \

--i-table feature-table-final.qza \

--m-metadata-file meta.tsv \

--p-min-features 10 \

--p-min-frequency 1500 \

--o-filtered-table filter-table.qza

可以查看qiime feature-table filter-samples的帮助文档来查看其他过滤参数，比如我可以使用参数--p-where只保留meta.tsv中的某一个分组。

接下来使用 ANCOM-BC方法来进行差异物种分析

#差异分析 cd /data3/sunjs/chenlina/69contorl_47PDAC_res

#--p-formula描述了微生物的绝对丰度是如何依赖于元数据中的变量的

#使用ANCOM-BC方法来进行差异物种分析

qiime composition ancombc \

--i-table filter-table.qza \

--m-metadata-file meta.tsv \

--p-formula group \

--o-differentials ancombc-group.qza

#生成差异丰富度分析结果的条形图可视化

qiime composition da-barplot \

--i-data ancombc-group.qza \

--p-significance-threshold 0.001 \

--o-visualization da-barplot-group.qzv

如何把.qza格式文件导出？比如把特征表导出

cd /data3/sunjs/chenlina/84contorl_50PDAC_res

qiime tools export \

--input-path feature-table-final.qza \

--output-path exported-feature-table biom convert \

-i exported-feature-table/feature-table.biom \

-o feature-table.tsv \

--to-tsv

这就把特征表导出成了tsv这种表格形式的文件

演示单端数据导入qiime2

先使用fastp对下载的数据进行质控

input_path=/data3/sunjs/chenlina/contorl_data output_path=/data3/sunjs/chenlina/contorl_clean_data

cd "$input_path" #SRR15884889.fastq for file in *.fastq; do

sample_id=$(echo "$file" | cut -d'.' -f1)

echo "sample_id是$sample_id

输入文件是$file"

echo "输出文件为${sample_id}.fastp.json和${sample_id}.fastp.html"

fastp -i "${file}" -o "${output_path}/${file}" \ -

-qualified_quality_phred 20 \

--length_required 50 \

--cut_front \

--cut_tail \

--trim_poly_g \

--html "${output_path}/${sample_id}.html" \

--json "${output_path}/${sample_id}.json" done

然后运行代码自动构建一个manifest.tsv文件

cd /data3/sunjs/chenlina/contorl_data

echo -e "sample-id\tabsolute-filepath" > contorl_manifest.tsv for file in *.fastq; do

echo -e "${file%.fastq}\t$(pwd)/$file" >> contorl_manifest.tsv done

解释下上面这段代码

上面的命令是用来生成一个清单文件（manifest file），这个文件将每个样本的ID和对应的绝对文件路径关联起来。这个文件随后可以被用来将单端FASTQ文件导入到QIIME 2中。

数据来源于文章 Dysbiosis of human gut microbiome in young onset colorectal cancer，他的数据中扩增子测序文件里面是单端测序文件，所以构建的manifest长这样子

然后把健康人的feces数据导入qiime2，由于是单端的所以

--type 参数是'SampleData[SequencesWithQuality]'

export TMPDIR=/data3/sunjs/

input_path=/data3/sunjs/chenlina/contorl_data/contorl_manifest.tsv

qiime tools import \

--type 'SampleData[SequencesWithQuality]' \

--input-path "${input_path}" \

--output-path contorl_demux-single-end.qza \

--input-format SingleEndFastqManifestPhred33V2