QIIME 2扩增子分析全流程详解：从数据导入到差异分析

搞微生物组研究的人几乎都绕不开扩增子分析。不管是16S、ITS还是18S只要是想看群落组成、多样性差异QIIME 2基本就是绕不过去的那道坎。这个流程我前前后后跑了几百次从最开始对着报错干瞪眼到后来能闭着眼把上游到下游一条龙跑完中间踩过的坑确实不少。这篇东西我想把整个qiime2扩增子分析流程从头到尾捋一遍包括为什么这么设计、每一步在干什么、参数怎么选、报错怎么排查尽量让刚接触的人能少走点弯路也让已经跑过流程的人能回头查漏补缺。说实话QIIME 2的官方教程已经写得很详细了但教程讲的是标准路径实际拿到自己的数据时总会遇到各种意外。比如导入报错、样本拆分后序列量对不上、DADA2跑完一看大部分序列被滤掉了这些问题教程里往往一笔带过但实际处理起来真要命。所以我这篇东西打算换个写法不按官方教程的顺序平铺直叙而是以拿到一批原始下机数据怎么一步步得到能发文章的结果为主线把整个流程拆开揉碎了讲重点放在那些官方文档不会明说的细节上。1. 为什么是QIIME 2这套流程的设计逻辑1.1 从QIIME到QIIME 2变化的不只是版本号很多人第一次接触QIIME 2时会下意识地拿它和老的QIIME也就是QIIME 1做对比。QIIME 1时代整个分析基本就是Python脚本串起来的一条流水线输入是文本格式的OTU表中间经过uclust、usearch这些工具聚类最后输出一堆表格和图。用起来简单粗暴但问题也很明显中间产物格式不统一分析步骤无法追溯换个环境跑结果可能就不一样。QIIME 2的架构几乎是推倒重来的。它把所有分析数据都封装成一种叫artifact的东西这个东西本质上是一个zip压缩包里面除了数据本身还带着完整的元数据和 provenance也就是这些数据是怎么一步步变出来的。这意味着每一步操作都会被记录任何人拿到你的artifact都能知道它经历了哪些处理这在科研可重复性上是一个质的提升。我最早用QIIME 2时其实很不习惯总觉得多了一层封装反而麻烦。但用久了就发现这种设计在合作和投稿时非常占便宜审稿人问你某个结果怎么来的你直接把provenance导出来发过去就行不用翻聊天记录和脚本历史。1.2 插件化架构你需要什么就装什么QIIME 2的另一个核心设计是插件化。它的功能不是一坨大杂烩而是拆分成一个个插件比如q2-dada2负责去噪q2-feature-table负责特征表操作q2-diversity负责多样性分析q2-taxa负责物种注释。每个插件各管一摊互不干扰还能单独更新。这种设计的实际好处在换环境时特别明显。比如你的服务器上装好了QIIME 2后来官方发布了新版本你不需要把整个环境推倒重装只需要更新某一个插件就行。再比如你只需要做16S分析完全不用装那些用不上的宏基因组插件环境能清爽不少。同时也带来一个需要适应的地方你得搞清楚每个步骤由哪个插件负责、参数格式是怎样的。刚开始会觉得命令复杂但写熟了之后反而觉得思路清晰每一步在做什么一目了然。1.3 整个分析流程的框架与我的推荐路线扩增子分析看起来步骤多但核心框架其实就三块上游质控导入、质量过滤、去噪、核心分析多样性、物种组成、下游统计差异分析、可视化。我跑的最多的标准路线是原始数据导入q2-tools-import→ 质控与去噪DADA2→ 多样性分析q2-diversity→ 物种注释q2-feature-classifier→ 差异物种分析ANCOM或DESeq2。其中DADA2在QIIME 2里集成了q2-dada2插件这也是目前主流的去噪方式替代了老流程里的OTU聚类。两者思路不同OTU聚类是按序列相似度通常是97%把reads归成一类而DADA2试图把测序错误和真实的生物学序列区分开得到的是单核苷酸精度的ASVAmplicon Sequence Variant。也就是说DADA2能区分只差一个碱基的序列聚类只能区分差异超过3%的序列分辨率上有本质差别。所以现在发文章的趋势也是ASV越来越主流OTU相对变少了。但要注意ASV和OTU各有适用范围你的样本量、测序深度、物种分辨率要求不同选择也会不一样。后面我会详细展开怎么选、怎么调参数。2. 准备工作从环境搭建到数据格式检查2.1 用conda装环境省一半的心除非你用的是QIIME 2官方提供的Docker镜像否则我强烈建议用conda或者mamba来安装。原因很简单QIIME 2的依赖关系非常复杂直接用pip装大概率会装出个跑不起来的残废环境来。安装时有一个重要的选择装哪个版本。不同版本的QIIME 2对应不同的Python版本和conda源这里一定要对照官方文档的版本说明来选。以我常用的2024.5版本为例创建环境的命令大概是这样的wget https://data.qiime2.org/distro/core/qiime2-2024.5-py38-linux-conda.yml mamba env create -n qiime2-2024.5 --file qiime2-2024.5-py38-linux-conda.yml注意文件里会写明Python版本比如py38下载时一定要跟你的系统架构匹配。如果你用的是Apple Silicon芯片的Mac需要找osx-arm64版本如果用的是Linux服务器一般是linux-64。装好之后每次使用前先激活环境conda activate qiime2-2024.5 qiime --help能正常打出帮助信息说明核心环境就绪了。这里有一个经常有人踩的坑conda装完后直接输qiime提示找不到命令大概率是没激活环境或者安装过程没走完看最后几行日志。2.2 输入数据盘点你的数据到底属于哪种情况做扩增子分析的数据来源很杂有公司返回的、有自有测序仪下机的、还有从公共数据库下载的。不管哪来的在导入之前都必须搞清楚一个关键问题你的数据是拆分好的还是没拆分的。所谓拆分好的demultiplexed指的是每个样本已经单独有一个fastq文件文件名里带了样本ID所谓没拆分的multiplexed是所有样本混在同一个fastq文件里需要靠barcode信息来拆分。这两种情况在QIIME 2里的导入方式是截然不同的。常见的是第一种也就是公司已经按样本返回了双端fastq文件。这时候你需要准备一个manifest文件这是一个CSV格式的表格列出每个样本的样本ID、正向序列文件路径、反向序列文件路径。格式长这样sample-id,forward-absolute-filepath,reverse-absolute-filepath sample1,/path/to/sample1_R1.fastq.gz,/path/to/sample1_R2.fastq.gz sample2,/path/to/sample2_R1.fastq.gz,/path/to/sample2_R2.fastq.gz注意路径必须是绝对路径相对路径在QIIME 2里经常出问题这是很多新手卡住的地方。还有一种情况是公司的数据直接按照Casava 1.8命名规则给好了文件名比如Sample1_S1_L001_R1_001.fastq.gz这种QIIME 2支持用CasavaOneEightSingleLanePerSampleDirFmt格式直接导入不用写manifest文件。但前提是目录结构必须规范一个样本一个子目录目录名就是样本名fastq文件放在里面。2.3 导入命令和导入后的数据验证以最常见的双端fastq加manifest为例导入命令长这样qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest.csv \ --output-path demux-paired-end.qza \ --input-format PairedEndFastqManifestPhred33--input-format这里有个重要细节Phred33还是Phred64。绝大多数Illumina下机数据都是Phred33编码但如果你用的是老旧的GA平台或者某些特殊流程可能输出Phred64。判断方法很简单用less打开fastq文件看质量值那行的字符如果出现!到I范围ASCII 33-73就是Phred33如果出现到h范围ASCII 64-104就是Phred64。弄错了后面所有质控都会出问题这个一定要确认。导入完成后强烈建议立刻跑一下可视化检查数据质量qiime demux summarize \ --i-data demux-paired-end.qza \ --o-visualization demux-paired-end.qzv qiime tools view demux-paired-end.qzv这一步会生成一个交互式网页能看到每个样本的reads数、序列长度分布、质量分数分布。拿到这个可视化结果后你需要做三件事第一确认样本总数和实际样本数一致第二看每个样本的reads量是否正常有没有特别少的第三看质量曲线判断后面DADA2的截断参数应该怎么设。这一眼检查能避免后面跑完了才发现某个样本数据有问题返工成本高得多。3. 核心分析步骤每一步的实操与参数逻辑3.1 DADA2去噪参数不是随手填的DADA2是整个流程里最核心、也最容易出问题的一步。它做的事情可以粗略理解为把测序错误和真实序列变异区分开去除嵌合体得到ASV表和代表序列。它在QIIME 2里的调用形式是qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux-paired-end.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f 240 \ --p-trunc-len-r 200 \ --p-n-threads 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qza--p-trunc-len是根据质量曲线来定的不是随便填的。基本原则是把质量分掉到20以下的区域截掉保留序列质量高的区域。比如正向序列在240bp以后质量开始下滑那--p-trunc-len-f就设240反向序列读长通常更短、质量更差所以--p-trunc-len-r往往要比正向短。如果你不确定先看demux summarize生成的质量图再决定。这里还有一个非常关键的参数容易被忽略--p-n-threads。默认-1表示用全部核心但如果你在共享服务器上跑最好指定一个合理的值避免把整台机器挤爆。我一般设成--p-n-threads 8左右不耽误事也不影响别人。还有--p-trim-left它代表从序列5端截掉多少个碱基主要是为了去掉引物序列。如果你的测序文件里没有去除引物而引物长度大约是20bp那就应该--p-trim-left-f 20。但大部分公司给的数据已经去掉了引物所以这里设为0通常没问题——前提是你确认过公司的数据说明。我见过很多人想当然地设一个值结果把有效序列截掉一大截后面分析直接没数据了非常可惜。跑完DADA2后会生成denoising-stats.qza这个文件记录了每个样本经过质量过滤、合并、去嵌合体后的reads数量。拿到它一定要看一眼如果发现某个样本的最终reads数只剩最初的10%那说明前面的参数设置或者样本本身质量有问题需要调整参数重跑。这里分享一个实操心得第一遍跑DADA2时可以先不追求最优参数用相对宽松的截断长度比如都设220快速跑一遍拿到大致结果看看每个样本的保留比例和总ASV数量再根据结果决定要不要调整参数。这样可以避免一遍遍试参数、每次等几小时的尴尬。尤其样本量大时DADA2跑一次很耗时先快后慢的节奏更合算。3.2 特征表过滤删掉那些没有意义的ASVDADA2跑完得到的table.qza是原始的ASV表里面会有很多仅在个别样本里出现个位数的ASV。这些ASV有可能是测序噪声或者极低丰度的稀有物种保留它们会干扰后续的多样性分析和差异分析所以一般要做一步过滤。QIIME 2提供了q2-feature-table插件最常用的过滤策略有两个按频率过滤和按样本元数据过滤。按频率过滤的做法是把在所有样本中总丰度低于某个阈值的ASV删掉。常见命令是qiime feature-table filter-features \ --i-table table.qza \ --p-min-frequency 10 \ --o-filtered-table table-filtered.qza--p-min-frequency 10的意思是保留在所有样本里总计数不低于10的ASV。这个值怎么选要看你的测序深度和样本量如果每个样本测了5万reads那总丰度至少为10的ASV是有意义的如果测序深度很浅比如每个样本只有5000 reads那阈值可以降到5甚至2。按样本元数据过滤则是在样本层面做筛选比如剔除某些异常样本。这个前提是你已经准备了sample-metadata.tsv文件里面记录了每个样本的分组信息。如果你想剔除掉特定几个样本可以用--p-sample-metadata-file配合过滤条件来实现。提醒过滤操作一定要保留一份过滤前的table作为备份后面如果要重新算不同阈值的结果直接回到备份接着做就行。我见过不少人过滤完就把原始表覆盖了事后要重新分析只能重跑DADA2白白浪费时间。3.3 多样性分析抽平深度的选择是门学问多样性分析包括alpha多样性样本内部的多样性和beta多样性样本间的差异QIIME 2里用q2-diversity插件完成。典型流程是qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table-filtered.qza \ --p-sampling-depth 11000 \ --m-metadata-file sample-metadata.tsv \ --output-dir core-metrics-results这个命令会一次输出所有核心多样性指标Observed Features、Shannon、Faiths PD、Unweighted UniFrac、Weighted UniFrac、PCoA图等。这里最关键的参数是--p-sampling-depth也就是抽平深度。它的含义是从每个样本中随机抽取相同数量的reads然后再计算多样性指标。为什么要抽平因为不同样本的测序深度不一样如果不抽平测序量大的样本会天然看起来多样性更高这是不公平的。抽平深度怎么选通常看demux summarize的结果找到所有样本中reads数最少的那个样本选择小于等于这个最小值的某个数。但也不能太低太低会丢失信息把测序较深的样本浪费掉。我的经验是找所有样本reads数的分布取大约10%到90%分位数之间的某个值或者用官方教程推荐的常见值11000。但注意如果某几个样本的reads数特别少比如只有3000那必须把抽平深度降到3000以下否则这些样本会在抽平后变成空样本直接报错。还有一个容易忽略的地方core-metrics-phylogenetic需要在之前构建好系统发育树也就是rooted-tree.qza。这个树的构建在q2-phylogeny插件里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 \ --p-n-threads 8如果跳过了建树这一步直接用table去跑core-metrics-phylogenetic会报错提示找不到树文件。很多人第一次跑都会卡在这里其实只要回到流程里把树建好就行。3.4 物种注释用对分类器结果差出一大截物种注释把ASV比对到参考数据库上得到每个ASV对应的分类学信息界、门、纲、目、科、属、种。QIIME 2里的标准做法是用q2-feature-classifier但它需要一个预先训练好的分类器classifier。官方提供了针对Greengenes 13_8和Silva 138等主流数据库的预训练分类器可以直接下载使用。但有一个问题官方分类器是用全长16S序列训练的如果你用的是某个特定高变区比如V3-V4区直接用它来分类精度会打折扣。正确的做法是下载参考序列和分类学信息用你自己扩增区域的引物来提取对应片段然后重新训练分类器。训练分类器的流程大致是# 提取引物对应区域 qiime feature-classifier extract-reads \ --i-sequences 85_otu_taxonomy_reference_seqs.qza \ --p-f-primer GTGCCAGCMGCCGCGGTAA \ --p-r-primer GGACTACHVGGGTWTCTAAT \ --o-reads ref-seqs-v3v4.qza # 训练分类器 qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs-v3v4.qza \ --i-reference-taxonomy 85_otu_taxonomy.qza \ --o-classifier v3v4-classifier.qza # 用训练好的分类器注释 qiime feature-classifier classify-sklearn \ --i-classifier v3v4-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza \ --p-n-jobs 4这一步的引物序列必须和你实际实验用的引物一致。很多人在这一步图省事直接拿官方预训练的分类器用结果注释出来的属水平一大半是unclassified然后又反过来质疑数据质量不好。其实问题往往出在分类器没有针对自己的扩增区域校准。提示如果你的样本是古菌、真菌ITS或18S请选择对应的参考数据库和分类器。用16S分类器去注释ITS数据结果基本是废的。3.5 差异丰度分析和可视化把结果变成能讲的故事前面这些步骤做完你已经有了多样性结果、物种组成结果接下来就是差异分析。QIIME 2自带ANCOM和ANCOM-BC但说实话很多人在实际投稿时会选择把ASV表导出来用R里的DESeq2、edgeR或者LEfSe做差异分析毕竟R生态系统里的工具更丰富、出图更灵活。把QIIME 2的结果导出到R用的命令是qiime tools export --input-path table.qza --output-path exported-table导出后会得到一个BIOM格式的文件feature-table.biom在R里用biomformat包读入或者用qiime2R包直接读取.qza文件方便得多。需要提醒的是做差异分析前要明确你的实验设计有几个分组、每组几个样本、是配对还是独立样本。ANCOM假阳性控制比较好但样本量太少时容易报错DESeq2对小样本相对友好但要求原始计数矩阵不能是抽平后的数据。具体用哪个要看你手头的数据情况。如果你不想切换到RQIIME 2也提供了qiime taxa barplot来做物种组成堆叠图qiime taxa barplot \ --i-table table-filtered.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv生成的可视化文件是交互式的可以直接按分组着色用来快速浏览数据非常方便。但它不是出版级图表发文时还是建议用R出图。4. 常见问题与排查技巧实录4.1 导入报错bad magic number和格式不匹配导入数据时报错bad magic number或者Failed to import是非常常见的情况。bad magic number通常意味着你的环境版本和.qza文件版本不匹配。比如你用2024.5版跑出的结果再用2021.8版去读大概率会报这个错。解决方案很简单统一版本或者在新环境里重新导入原始数据。再有一种情况是--input-format选错了。比如你给的是PairedEndFastqManifestPhred33但实际fastq质量编码是Phred64或者你把单端数据用双端格式去导入。这种报错信息往往不会直接告诉你是格式问题而是提示某个字段解析失败。排查思路是先确认数据是单端还是双端再确认质量编码类型最后再检查manifest文件路径和格式。4.2 样本拆分后reads数对不上导入时一切正常但demux summarize一看某个样本的reads数比公司报告里写的少很多甚至为0。这种情况多发生在用multiplexed方式导入时barcode和样本对应关系搞错了。解决办法是重新检查barcode序列和样本的对应表确认有没有barcode缺失、barcode序列写反正向反向的情况。如果你用的是公司已经拆分好的数据样本reads数骤减的原因更可能是文件本身不完整。比如某些公司按样本返回数据时样本量大的会拆成多个lane文件夹你没把所有的fastq文件都合并进来。这种情况需要把同一个样本的多个lane文件合并后再导入。4.3 DADA2跑完序列怎么全被滤掉了这是最让人崩溃的报错之一DADA2跑完了denoising-stats显示每个样本的最终reads数都是0或者接近0。出现这种情况最常见的原因是--p-trunc-len设置得太长导致正向和反向序列拼接后重叠区太短甚至为负数DADA2合并时全部失败。比如你的读长是150bp但你把--p-trunc-len-f设成了150--p-trunc-len-r也设成了150两条150bp的序列如果重叠区少于一定的阈值默认12bp就没办法合并。解决办法是把截断长度调短比如正向设140、反向设130或者设成0不截断看合并率是否恢复。其次要检查--p-trim-left是否设得过大。如果引物已经去除了你再截掉20bp序列太短也会导致合并失败。有一招很实用先用小样本量快速测试不同参数组合比较denoising-stats里的final reads数选保留率最高的那组参数。样本多的时候能省不少时间。4.4 可视化文件打不开或交互图加载慢qiime tools view生成的.qzv文件其实是网页格式默认需要用浏览器打开。如果你在服务器上跑分析本机没有图形界面需要把.qzv文件下载到本地再用浏览器打开。如果你的.qzv文件特别大比如样本量几千个浏览器加载时可能会卡住。这时候建议缩小可视化范围比如只对其中一个分组做barplot或者用非交互式的静态图片方式导出。另外新版QIIME 2的可视化文件对浏览器版本也有要求太旧的浏览器可能打不开交互组件。遇到这种情况换用Chrome或者Firefox最新版通常能解决。4.5 环境里多个QIIME版本共存这个情况很常见服务器上之前装了2021.8后来又装了2023.9两个版本都还在用。问题是两个版本的.qza格式可能不兼容有时候你明明激活的是新环境却读了旧环境跑出来的.qza就会报格式错误。我的建议是每个项目固定用一个环境并且在整个项目周期内不要升级。如果你确实需要切换版本务必把项目的所有中间文件table.qza、rep-seqs.qza、taxonomy.qza都用同一个版本重新跑一遍避免混用。5. 我踩过坑之后总结出的几条经验跑这套流程久了最深的体会是扩增子分析流程本身并不复杂难的是每一步都要对数据保持敏感不能机械地执行命令。第一拿到数据先花半小时做数据体检。看看每个样本的reads量、序列长度分布、质量分布养成习惯。这一步能提前发现很多问题比如某个样本测序失败、部分样本被污染、引物没有去除干净等。第二保存每一步的中间产物。QIIME 2本身有provenance记录确实很好但实际操作中还是建议把你觉得重要的中间文件单独复制到一个文件夹里做好命名和日期标记。毕竟.qza文件是有版本依赖的万一你重装了环境旧文件可能就废了。第三样本元数据文件一定要提前准备规范。sample-metadata.tsv的格式要求其实比较严格第一列必须是sample-id第二列必须有一个非空列可以包含分类信息和数值信息。格式不对的话后面很多分析和可视化都会报错。第四不要一口气把所有分析跑完再去检查结果每一步都停下来看一眼输出。DADA2跑完看denoising-stats多样性跑完看PCoA图注释完看barplot。只有每一步的结果都合理最终分析才站得住脚。第五善用官方论坛和GitHub Issues。QIIME 2的社区活跃度很高大多数报错信息都能搜到前人遇到过的解决方案。搜索时直接用报错原文去搜不要自己翻译成自然语言效果会好很多。扩增子分析这个领域变化挺快的新的插件、新的参考数据库、新的统计方法层出不穷。但底层逻辑没变理解每个步骤在干什么、参数为什么这么设、数据为什么长这样。把基本功打扎实了不管工具怎么换你都能快速上手。