企业官网建设流程全解析

搞微生物组研究的人，几乎都绕不开扩增子分析。不管是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试图把测序错误和真实的生物学序列区分开，得到的是单核苷酸精度的ASV（Amplicon 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 特征表过滤：删掉那些没有意义的ASV

DADA2跑完得到的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、Faith's 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对小样本相对友好，但要求原始计数矩阵，不能是抽平后的数据。具体用哪个，要看你手头的数据情况。

如果你不想切换到R，QIIME 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的社区活跃度很高，大多数报错信息都能搜到前人遇到过的解决方案。搜索时直接用报错原文去搜，不要自己翻译成自然语言，效果会好很多。

扩增子分析这个领域变化挺快的，新的插件、新的参考数据库、新的统计方法层出不穷。但底层逻辑没变：理解每个步骤在干什么、参数为什么这么设、数据为什么长这样。把基本功打扎实了，不管工具怎么换，你都能快速上手。