简介这份资源面向缺乏生物信息学背景的微生物组研究者与初学者提供使用QIIME 2流程分析16S rRNA基因扩增子测序数据的完整学习材料帮助读者快速掌握从数据导入到可视化呈现的标准分析路径。资源包为单个docx文档约578KB内容涵盖软件安装环境选择、特征表生成、α与β多样性分析、物种组成与差异物种分析等核心环节并附有配套视频、分析代码、测序数据与预期结果便于对照复现。文档还系统总结了安装与使用中的常见问题及解决策略对参数优化和硬件配置给出具体建议例如推荐4核CPU、16GB内存及大于原始数据三倍的硬盘空间。目前已有2792人学习下载适合希望系统入门微生物组扩增子分析、需要可复现操作指引的研究人员参考。1. 从一份能跑通的 QIIME 2 流程说起16S 扩增子分析到底难在哪如果你刚拿到一摞 Illumina 双端 fastq样本表还没整理导师又催着要 α 多样性箱线图和门水平堆叠柱状图那 QIIME 2 大概率是你绕不开的一站。它用 Python 3 重写插件化、可交互、结果可复现是目前 16S rRNA 基因扩增子分析里引用量最高的流程之一。但真上手你会发现两个现实问题一是它压根不支持 Windows 直接装二是官方文档十万字起步新手光看安装就能劝退。这份资源的价值就在这——它把从 Miniconda 装环境、GreenGenes 建分类器到 DADA2 降噪、Alpha/Beta 多样性、ANCOM 差异分析整条链路串成了一份可复现的脚本还配了示例数据和视频。适合谁做微生物组、植物根系、肠道菌群这类课题手里有双端测序数据、想自己跑一遍完整流程的人。下面我按实际拆包顺序讲重点放在参数怎么定、哪一步最容易翻车。2. 环境部署与数据库准备Miniconda、QIIME 2 与 GreenGenes 分类器2.1 为什么优先选 Linux 服务器或 WSL而不是虚拟机QIIME 2 官方只发 Linux 和 Mac 的 conda 包Windows 用户只有两条路WSLWindows Subsystem for Linux或 VirtualBox 虚拟机。资源里明确推荐 WSL 的 Ubuntu 20.04 LTS虚拟机标了“不推荐效率低”。原因很直接DADA2 降噪和分类器训练都是吃 CPU 和内存的活虚拟机多一层硬件抽象同样的数据可能多花一倍时间。我一般建议数据量小于 20G、样本几十个WSL 足够上百样本或要反复训练 V 区分类器直接上 4 核 16G 起的 Linux 服务器硬盘留原始数据 3 倍以上因为中间 qza 文件会膨胀。WSL 安装本身在微软商店搜 Ubuntu 20.04 LTS 即可装完在开始菜单启动就是命令行。Mac 用户系统自带终端能直接 ssh 远程服务器但资源里对 Mac 本地安装标了“兼容性问题较多”我的经验是 Mac 本地跑小数据可以大数据还是连服务器稳。2.2 Miniconda 与 QIIME 2 2020.6 的安装命令逐行拆解整个安装的核心是 conda 环境隔离。先装 Miniconda3再下载 QIIME 2 的 yml 环境文件最后 create 环境。命令如下# 下载并安装 Miniconda3已装 conda 可跳过 wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh ~/miniconda3/bin/conda init # 下载 QIIME 2 2020.6 的环境依赖列表 wget -c https://data.qiime2.org/distro/core/qiime2-2020.6-py36-linux-conda.yml # 新建名为 qiime2-2020.6 的环境并安装 conda env create -n qiime2-2020.6 --file qiime2-2020.6-py36-linux-conda.yml # 每次分析前激活环境 conda activate qiime2-2020.6逻辑说明conda init把 conda 写进 shell 配置之后新开终端才能直接用conda activate。conda env create会按 yml 里的版本号拉取所有依赖这一步网络不稳容易断断了重跑即可conda 有缓存。参数上环境名qiime2-2020.6建议和版本号一致以后装新版本不会互相覆盖。注意激活环境后命令行提示符前会出现(qiime2-2020.6)没出现说明没激活成功后面所有qiime命令都会报 command not found。2.3 GreenGenes 数据库导入与全长/V 区分类器训练分类器是物种注释的“字典”训练一次可以反复用。资源里给了两条路全长通用分类器和指定 V 区分类器。全长耗时约半小时V 区以 V5-V7 为例约 9 分钟提取序列加 8 分钟训练。先导入参考序列和物种分类# 下载并解压 GreenGenes 13_8 wget -c ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz tar -zxvf gg_13_8_otus.tar.gz # 导入 99% 聚类代表序列 qiime tools import \ --type FeatureData[Sequence] \ --input-path gg_13_8_otus/rep_set/99_otus.fasta \ --output-path 99_otus.qza # 导入物种分类信息 qiime tools import \ --type FeatureData[Taxonomy] \ --input-format HeaderlessTSVTaxonomyFormat \ --input-path gg_13_8_otus/taxonomy/99_otu_taxonomy.txt \ --output-path ref-taxonomy.qza--type指定数据类型FeatureData[Sequence]是序列FeatureData[Taxonomy]是分类表。--input-format HeaderlessTSVTaxonomyFormat是因为 GreenGenes 的分类文件没有表头不写这个参数导入会报格式错。导入完成后训练全长分类器time qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads 99_otus.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier classifier_gg_13_8_99.qzatime是看耗时的不是必须。fit-classifier-naive-bayes用的是朴素贝叶斯分类器这是 QIIME 2 默认也是 16S 最常用的。如果你的引物是特定 V 区强烈建议训练特异分类器精度会明显提升。以 V5-V7 为例# 按引物提取对应区段序列 time qiime feature-classifier extract-reads \ --i-sequences 99_otus.qza \ --p-f-primer AACMGGATTAGATACCCKG \ --p-r-primer ACGTCATCCCCACCTTCC \ --o-reads ref-seqs.qza # 基于提取序列训练特异分类器 time qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier classifier_gg_13_8_99_V5-V7.qza--p-f-primer和--p-r-primer必须和你实验用的引物完全一致这里 V5-V7 用的是 799F/1193R。写错引物会导致提取不到序列或提取错区段分类结果全乱。这一步的坑我后面单独讲。3. 从原始 fastq 到特征表数据导入与 DADA2 降噪3.1 manifest 文件的生成与数据导入QIIME 2 不直接吃一堆 fastq它要一个 manifest 文件里面三列sample-id、forward-absolute-filepath、reverse-absolute-filepath。资源里用 awk 从 metadata 批量生成# 根据 metadata 生成 manifest注意 $PWD 是当前绝对路径 awk NR1{print sample-id\tforward-absolute-filepath\treverse-absolute-filepath} \ NR1{print $1\t$PWD/seq/$1_1.fq.gz\t$PWD/seq/$1_2.fq.gz} \ metadata.txt manifestNR1处理表头NR1处理数据行。$PWD必须用绝对路径QIIME 2 对相对路径支持不好容易报找不到文件。生成后导入qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2--input-format选PairedEndFastqManifestPhred33V2对应 Illumina 的 Phred33 质量值编码这是目前最常见的。资源里提到 1G 的 fq 导入要 7 分钟但 fq.gz 压缩格式只要 34 秒所以原始数据别解压直接导 gz。注意如果你的数据是混池测序没拆样本得先让测序公司拆分或者用 QIIME 1 的脚本拆QIIME 2 不负责拆分。3.2 DADA2 降噪参数trim 和 trunc 怎么定DADA2 是 QIIME 2 里生成特征表ASV 表的核心它不聚类直接去噪。关键参数是--p-trim-left-f/r和--p-trunc-len-f/r。资源里给的示例是 trim-left-f 29、trim-left-r 18、trunc-len 都设 0time qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-n-threads 8 \ --p-trim-left-f 29 --p-trim-left-r 18 \ --p-trunc-len-f 0 --p-trunc-len-r 0 \ --o-table dada2-table.qza \ --o-representative-sequences dada2-rep-seqs.qza \ --o-denoising-stats denoising-stats.qza--p-trim-left-f 29是切掉正向引物及前面碱基--p-trim-left-r 18切反向。这两个值取决于你引物长度不是固定值。--p-trunc-len设 0 表示不截断但实际项目中我一般会看质量分布图在质量掉到 Q20 以下的位置截断否则末端错误碱基会进 ASV。--p-n-threads是线程数资源里给了实测96 线程 34 分钟24 线程 44 分钟8 线程 77 分钟1 线程 462 分钟。线程不是越多越快超过物理核数反而有调度开销一般设物理核数的 80% 左右。跑完把 dada2 结果复制成主流程用的名字cp dada2-table.qza table.qza cp dada2-rep-seqs.qza rep-seqs.qza3.3 特征表统计与抽平阈值的确定feature-table summarize生成 table.qzv用 view.qiime2.org 打开看每个样本的测序量分布qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata.txt资源里示例数据最小值 27060第一分位数 28581中位数 30867分布很均匀所以抽平阈值直接取最小值 27060。但如果最小值是个别样本的异常低值和第一分位数差很多就得在最小值和第一分位数之间选尽量保留更多样本和更多数据。注意低于阈值的样本会被丢弃不参与多样性分析所以阈值定太高会丢样本定太低会浪费数据。资源里特别提醒抽平最小值 1000 是 454 时代的标准现在 Illumina 通量高最小值一般不低于 5000推荐 1 万起。4. 多样性与物种组成分析Alpha、Beta、注释与差异4.1 进化树构建与 core-metrics 多样性计算Alpha 和 Beta 多样性里Faiths PD 需要进化树所以先建树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.qzaalign-to-tree-mafft-fasttree一步完成比对、屏蔽高变区、建树、生根。输出里rooted-tree.qza后面要用。然后跑核心多样性qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 27060 \ --m-metadata-file metadata.txt \ --output-dir core-metrics-results--p-sampling-depth就是抽平阈值这里用 27060。这一步会生成 4 种 Alpha 指数faith_pd、shannon、observed_features、evenness和 4 种 Beta 距离矩阵unweighted_unifrac、bray_curtis、weighted_unifrac、jaccard还有对应的 PCoA 结果。抽平是随机重采样所以每次跑结果会有微小差异但样本量够大时不影响结论。4.2 Alpha 多样性组间显著性检验以 observed_features 为例indexobserved_features qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/${index}_vector.qza \ --m-metadata-file metadata.txt \ --o-visualization core-metrics-results/${index}-group-significance.qzv结果里有箱线图和 Kruskal-Wallis 两两比较的 p 值和 q 值。资源示例里 KO vs WT 的 p0.010139q0.030417显著OE vs WT 的 p0.054241q0.081362不显著。这里要注意 q 值是 FDR 校正后的比 p 值严格报告时优先看 q 值。index变量可以换成 faith_pd、shannon、evenness命令结构一样。4.3 Beta 多样性 PCoA 与 PERMANOVA 检验Beta 多样性先看 PCoA 图再做成对 PERMANOVAdistanceweighted_unifrac columnGroup qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/${distance}_distance_matrix.qza \ --m-metadata-file metadata.txt \ --m-metadata-column ${column} \ --o-visualization core-metrics-results/${distance}-${column}-significance.qzv \ --p-pairwise--p-pairwise开启两两比较不开启只做整体检验。资源示例里三组两两比较 p 值都小于 0.05q 值也小于 0.05说明组间群落结构差异显著。PCoA 图在weighted_unifrac_emperor.qzv里可以调颜色、形状、透明度导出 SVG 矢量图。注意PERMANOVA 的置换检验很耗时指定--m-metadata-column只比较关心的分组能省不少时间。4.4 物种注释、堆叠柱状图与 ANCOM 差异分析物种注释用之前训练的分类器qiime feature-classifier classify-sklearn \ --i-classifier classifier_gg_13_8_99_V5-V7.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzvclassify-sklearn对每个代表序列做分类输出特征 ID、分类结果和置信度。堆叠柱状图qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file metadata.txt \ --o-visualization taxa-bar-plots.qzv在 qzv 里可以切分类级别、改配色、按分组和丰度排序导出 SVG 和 CSV。差异分析用 ANCOM# 添加伪计数ANCOM 要求无零值 qiime composition add-pseudocount \ --i-table table.qza \ --o-composition-table comp-table.qza # 在属水平合并后再做 ANCOM qiime taxa collapse \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-level 6 \ --o-collapsed-table table-l6.qza qiime composition add-pseudocount \ --i-table table-l6.qza \ --o-composition-table comp-table-l6.qza time qiime composition ancom \ --i-table comp-table-l6.qza \ --m-metadata-file metadata.txt \ --m-metadata-column Group \ --o-visualization ancom-Group-l6.qzv--p-level 6是属水平界门纲目科属第 6 级。ANCOM 输出火山图显著 ASV 或属可以在 taxonomy.qzv 里查分类。资源示例里组间差异小只有 1 个显著 ASV属水平结果更适合结合生物学讨论。5. 避坑与常见问题排查从安装到分析的 5 个血泪教训5.1 conda 环境创建卡在 solving environment现象conda env create跑很久不动或者报Solving environment: failed。原因conda 默认源在国外网络不稳或者 yml 里某个包版本冲突。解决换国内镜像源或者用mamba替代 conda 解依赖mamba 快很多。实在不行下载 yml 后手动删掉版本号太死的包让 conda 自己解。5.2 分类器训练报内存不足现象fit-classifier-naive-bayes跑到一半被 kill日志显示 OOM。原因全长 GreenGenes 99% 序列约 100 万条朴素贝叶斯训练吃内存16G 机器跑全长容易爆。解决优先训练 V 区特异分类器序列数少很多或者加内存到 32G或者用 97% 聚类而非 99%序列数少一个量级但分辨率会降。5.3 DADA2 输出特征数异常少或异常多现象denoising-stats 里大部分 reads 没进 ASV或者 ASV 数量几千上万。原因trim/trunc 参数不对。trim-left 切多了切到保守区切少了引物没去干净trunc-len 设 0 不截断末端错误碱基被当成真实变异ASV 虚高。解决先看 demux.qzv 里的质量分布图正向在质量掉到 Q20 的位置截断反向同理trim-left 按引物长度加 1-2 个碱基设。资源里示例 trim-left-f 29、r 18 是匹配它引物的别照抄。5.4 抽平后样本丢失现象core-metrics 跑完发现样本数比原来少。原因--p-sampling-depth设得比某些样本的测序量高那些样本被丢弃。解决在 table.qzv 的 Interactive Sample Detail 里看每个样本的测序量阈值设在最小值和第一分位数之间兼顾样本保留和数据利用。如果丢的样本是关键组宁可降低阈值。5.5 qzv 文件打不开或显示空白现象view.qiime2.org 上传 qzv 后一直转圈或空白。原因qzv 本质是 zip里面是网页和图表浏览器兼容性或文件损坏。解决先确认文件大小正常没下完的重新生成换 Chrome 或 Firefox本地解压 qzv 看里面 index.html 能不能打开。注意qzv 不要用文本编辑器改会破坏结构。6. 进阶技巧用 V 区特异分类器把物种注释精度再提一档如果你已经跑通全长分类器下一步最值得投入的就是训练实验特异的 V 区分类器。原理不复杂GreenGenes 全长序列里不同 V 区的保守程度不一样用全长训练时分类器学到的特征有一部分来自你根本没扩增的区域反而引入噪声。用extract-reads按你的引物把参考序列切成对应区段再训练分类器只学你测到的那段精度通常能提升几个百分点尤其是属水平。具体操作上引物序列必须和实验完全一致正向反向都不能错。资源里 V5-V7 用的是 799FAACMGGATTAGATACCCKG和 1193RACGTCATCCCCACCTTCC。提取时如果引物有简并碱基QIIME 2 支持 IUPAC 简并符号直接写就行。提取完先看ref-seqs.qza的序列长度分布正常应该在目标区段长度附近如果长度分布很散说明引物匹配有问题得检查引物方向或序列。训练完的分类器建议用一小批已知物种的序列做个 sanity check拿几条代表序列跑classify-sklearn看注释结果和预期是否一致。我一般会留一个全长分类器做对照如果 V 区分类器结果反而更差说明引物或提取参数有问题别硬用。还有一个容易被忽略的点ANCOM 在属水平做差异分析时taxa collapse的--p-level 6对应属但 GreenGenes 的分类层级里界门纲目科属种是 7 级第 6 级是属第 7 级是种。如果你要种水平设 7但种水平注释置信度通常很低不建议。ANCOM 的伪计数步骤不能省表里有零值会直接报错。最后说个习惯每次跑完 DADA2我都会把 denoising-stats.qzv 打开看 input、filtered、denoised、merged、non-chimeric 五个数字的漏斗。如果 filtered 掉太多说明 trim/trunc 太狠如果 non-chimeric 掉太多说明嵌合体比例高可能是 PCR 循环数太多。这个漏斗比最终 ASV 数更能反映数据质量。从那以后我每次拿到新数据都强制先跑一遍 denoising-stats 再决定后续参数希望帮到你。本文还有配套的精品资源点击获取
