大型网站建设报价方案深圳建站软件
目录
创建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