企业官网建设流程全解析

别人拿到测序公司交付的fastq数据包，解压后几百个.gz文件，第一反应通常是找教程从导入数据开始跑。但QIIME2和传统分析工具最大的不同是：它不直接吃fastq，你得先理解它的artifact机制，否则会在导入这一步反复报错，连数据都进不了门。这篇博文我按一条龙流程完整走一遍——导入数据、质控去噪、多样性分析、分类学注释、差异丰度与功能预测——每一步我都会给出命令、讲解参数取舍，顺便把容易坑人的细节一次性交代清楚。

博文正文

1. 开工之前：先理解QIIME2的数据组织和运行环境

1.1 artifact不是格式洁癖，是为了数据追溯

先不谈命令，我把话放前面：QIIME2（发音是“chime two”）之所以要引入.qza这种封装格式，不是故意增加学习门槛，而是要解决生物信息项目里最让人头疼的“溯源问题”。举个实际例子：我把一个table.qza发给合作者，里面嵌入了完整的产生历史，包含上游去了多少噪音、用了哪些参数、输入文件是什么，谁都可以沿着provenance（数据谱系）往前追溯。这一点在写论文Methods部分时特别爽——你不需要单独维护一张Excel记录每一步参数，qiime tools界面里直接能查看完整操作日志。

理解了这层逻辑，开篇那些看似繁琐的概念就顺了：.qza是“分析工件”的容器，.qzv是可视化的容器，两者类似一个带说明书的数据包。每次操作都会生成新的对象，而你之前的原始数据始终不会被覆盖，中间每一步都有迹可循。所以这里的建议是：给项目建一个清晰的目录结构，建议按input、scripts、output分成三层，避免后续分析时被一堆.qza淹没。

1.2 环境准备：conda、版本选择与资源预判

QIIME2官方推荐用Miniconda安装，避免污染系统Python环境，这一点我强烈建议照做。下面是安装的大致流程：

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh conda update conda conda create -n qiime2-2023.5 python=3.8 conda activate qiime2-2023.5

这里的版本号重要。每个版本对应的Python版本和插件版本捆绑较死，我不建议你随便用pip去补装插件，尽量整体安装官方release。以我常用的2023.5版本为例，整个环境装完约占用几个GB空间，跑dada2去噪时内存消耗看样本量，常规几十个样本8-16GB内存也能跑，但如果数据量上了百级PE reads，建议准备32GB以上，并配合多核CPU用。

conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/qiime2/ qiime2

如果在国内环境，把官方conda源替换成清华镜像会省很多时间，但要注意镜像列表里必须包含qiime2这个channel。另外一定记得确认环境中没有同时存在另一个Python版的q2cli，否则导入命令会提示找不到可执行文件。

1.3 拿到测序数据后的第一件事：摸清R1/R2与目录结构

Illumina双端测序下机后，每个样品会有R1和R2两个fastq.gz文件，R1是正向读长，R2是反向读长，二者中间是插入片段。很多第一次做16S的人容易犯一个错误：直接把两个文件单独导入，当作单端数据处理，这样会把双端信息割裂，后续dada2拼接就没有意义。正确思路是：一对R1/R2在导入后代表同一个样本的两个方向，样本数量应该是文件对数，而不是文件个数。

先解压并看一眼目录：

ls rawdata/ | head

正常情况你会看到Sample001_R1.fastq.gz、Sample001_R2.fastq.gz这种成对文件，第三方公司还常带barcode和adapter信息。这里务必先把样本名统一，只保留字母、数字和下划线。如果样本名下划线太多，后续分组信息会写得很痛苦。我最开始跑项目时，公司文件叫A-1_S1_L001_R1_001.fastq.gz，光是把名字改成A1.R1.fastq.gz就折腾了半天。

2. 数据导入：manifest文件、常见报错与验证

2.1 两种导入姿势：Casava目录格式与manifest清单

QIIME2支持多种导入格式，最常见的是CasavaOneEightSingleLanePerSampleDirFmt和manifest。Casava格式适合测序公司交付的原始目录还保留着标准命名的场景——目录层级固定为SampleID_SampleNumber_L001_R1_001.fastq.gz，此时可以直接：

qiime tools import \ --type CasavaOneEightSingleLanePerSampleDirFmt \ --input-path rawdata \ --output-path demux.qza

但我遇到更多的情况是：公司交付的文件名已经改成了自己的命名习惯，或样本来自多个批次。这时老老实实用manifest清单。先在文本编辑器里生成一个manifest.tsv，三列：sample-id、forward-absolute-filepath、reverse-absolute-filepath。

sample-id forward-absolute-filepath reverse-absolute-filepath A1 /absolute/path/rawdata/A1_R1.fastq.gz /absolute/path/rawdata/A1_R2.fastq.gz A2 /absolute/path/rawdata/A2_R1.fastq.gz /absolute/path/rawdata/A2_R2.fastq.gz B1 /absolute/path/rawdata/B1_R1.fastq.gz /absolute/path/rawdata/B1_R2.fastq.gz

然后导入：

qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path manifest.tsv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2

这里有个细节需要说明：Phred33和Phred64两种格式容易混淆。现在主流Illumina下机数据基本都是Phred33，默认选PairedEndFastqManifestPhred33V2不会错。如果分析老数据，不确定时用seqkit head检查质量值字符范围，ASCII码43以下对应Phred33，43以上则是Phred64，再对应修改格式。

2.2 manifest文件的隐形坑：绝对路径、制表符和CRLF

很多人在导入时报File not found或解析错误，多半不是数据问题而是路径问题。manifest里必须写绝对路径，而且要保证路径里的目录和文件真的有读权限。服务器环境下，如果是从Windows下编辑的manifest，记事本默认换行符是CRLF，而Linux下解析时会残留\r字符，导致样本ID后面多个不可见符号。解决办法是：

sed -i 's/\r$//' manifest.tsv

还有一种坑是：字段分隔符到底用tab还是逗号。PairedEndFastqManifestPhred33V2严格识别TSV（制表符分隔），而旧的PairedEndFastqManifestPhred33虽也是TSV，但如果从Excel复制粘贴会把tab换成空格，导致解析字段数不对。我的习惯是：用VSCode或Notepad++先开启“显示所有字符”来确认分隔符，再上传服务器。

2.3 导入后第一件事：查看每个样本的reads数

导入完成不代表数据好用，你还需要马上验证：

qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv

然后把demux.qzv拖到QIIME2 View网页里查看。这里主要看两个信息：一是每个样本的序列条数，判断有没有样本明显偏少，二是质量分数图。如果看到某个样本reads数只有几千，而其他样本十几万，尽早记录样本名，后续多样性分析时要么剔除，要么做好心理准备。这一步不需要复杂统计，纯粹是“提前认清数据”。

质量分数图同样值得看：坐标横轴是读长位置，纵轴是质量分数，会分别显示forward和reverse。我通常在导入后只判断一件事：总体质量是不是断崖式下降，如果是，下一步DADA2的截断参数就要重点照顾这段位置。

3. 质量控制与DADA2去噪：参数别乱抄，要看图

3.1 质量图到底怎么读：剪哪里留哪里

从demux.qzv里可以看到两条曲线：Forward（R1）和Reverse（R2）。R1的质量通常前几十个碱基不太好（测序仪启动不稳定），中间稳定，末端开始下降；R2整体质量低于R1，且末端下降更早更快。你的任务就是确定：左端要剪掉多少（trim-left），右端从哪个位置开始切（trunc-len）。

不要照抄网上的--p-trunc-len-f 240 --p-trunc-len-r 200。不同引物、不同测序平台和不同插入片段长度会导致最优参数差很多。我自己的判断逻辑是：

• 观察质量曲线，找到质量分数中位数开始低于30的那个位置。
• 给“重叠区”留足余量。双端reads后续dada2要合并，如果R1/R2截得太短，两个方向没有足够重叠，会直接导致拼接失败。一般要求截断后至少留有20-50bp重叠，长度低于100bp会使dada2默认拒绝合并。
• 不要把曲线尾部硬保，质量差的碱基会显著增加错误ASV。

3.2 DADA2参数选择：trim-left与trunc-len的执行逻辑

执行去噪的命令如下：

qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.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

trim-left-f和trim-left-r是把每个read前几个碱基剪掉，目的是去除引物序列。如果测序完成后公司已经去引物，那这里就用0。如果没去引物，按你所用引物的长度来设置，比如515F/806R这对V4区引物，通常正向剪掉19bp、反向剪掉20bp，但具体以引物序列的实际长度为准。注意这里有个原则：剪错左端比截错右端更难受，因为左端是结构性序列，剪多了会整体偏移，后续比对注释全乱。

trunc-len-f和trunc-len-r是从read的3‘端截断到指定长度。这个参数越接近质量图显示的“下降点”越好。比如forward在240bp后质量明显下滑，就设240；reverse在200bp后下滑，就设200。设得太大会保留大量错误碱基，设得太小又会损失重叠区信息。我见过一个项目，reverse只有180bp有效长度，有人硬按教程设了250，结果合并成功率只有20%，重跑后设置成180，成功率回到85%。

3.3 去噪后必须看的两个数量指标

DADA2跑完会产出三个文件：table.qza（特征表）、rep-seqs.qza（代表序列）、denoising-stats.qza（统计信息）。先查看denoising-stats：

qiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv

这个表会列出每个样本的input（输入reads数）、filtered（过滤后）、denoised（去噪后）、merged（合并后）、non-chimeric（去嵌合体后）。我一般会算一个比例：non-chimeric除以input，正常情况在50%-80%之间。如果低于40%，先考虑截断参数是否太保守，再考虑数据质量是不是本身有问题。

紧接着查看特征表：

qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv

这里最关键是看每个样本的“Frequency per sample”分布。如果样本之间测序深度差一个数量级，比如最少的才3000，最多的30万，那下一步多样性分析的采样深度选择会变得很难受。我会在这个环节记录两个数：最小样本reads数，和所有样本频率分布的中位数，后续core-metrics选择采样深度就用得上。

4. Alpha与Beta多样性分析：样品组的差异从哪里看出来

4.1 先建系统发育树再做多样性，不是多此一举

有些做过旧流程的人会疑惑：为什么QIIME2的多样性分析之前要先跑q2-phylogeny？因为UniFrac距离需要系统发育树提供物种间的亲缘信息，单纯Bray-Curtis不需要树，但要做加权/非加权UniFrac就必须有树。一条龙流程里，这一步几乎是必跑的：

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-unrooted-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza

这里我不建议换参数。MAFFT用来多序列比对，FastTree用来构建进化树，对于16S这种保守区域组合，这组工具配合得最稳。构建完成后，rooted-tree.qza在UniFrac计算中会被用到，而unrooted-tree.qza在需根分析中同样有用。全部保留即可。

4.2 Alpha多样性：香农指数、Chao1与抽平深度的逻辑

执行核心多样性pipeline：

qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 10000 \ --m-metadata-file sample-metadata.tsv \ --output-dir core-metrics-results

这里最容易犯错的是sampling-depth——抽平深度。这个值的意思是：为了让每个样本的测序深度可比，把每个样本随机抽到同样数量的reads。如果设成10000，那所有低于10000的样本都会被丢弃。选太高，样本量骤减；选太低，数据浪费严重。我的基本原则是：先看table.qzv里的频率分布，取一个能保留80%以上样本的深度值。如果样本量少或者实验设计不允许丢样本，宁可选低一点的深度，也要保留更多生物学重复。

Alpha多样性结果重点关注这几个指标：observed_features（观察到的ASV数）、shannon（香农指数）、faith_pd（基于系统发育的多样性）、evenness（均匀度）。我一般先看shannon，因为它兼顾物种丰富度和均匀度。如果想看每组之间的差异是否显著，可以继续用Kruskal-Wallis检验：

qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/shannon_vector.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization alpha-shannon-group-significance.qzv

4.3 Beta多样性：PCoA图怎么讲出生物学故事

Beta多样性的核心是样本间差异距离。QIIME2默认输出四个距离矩阵：unweighted_unifrac（只看物种有无差异）、weighted_unifrac（同时考虑物种相对丰度）、bray_curtis、jaccard。实操里我更关注weighted和unweighted两个UniFrac：前者反映群落组成的量变，后者反映群落组成的质变。如果两组分开但两个指标结论不一致，往往说明差异来自丰度低的稀有物种，而不是优势物种变化。

直接用主坐标分析（PCoA）画图太常见，但你得会看：

qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization beta-weight-unifrac-group-significance.qzv

在QIIME2 View里看PCoA图时，我常做三件事：第一，按分组着色，判断组间是否形成明显聚类；第二，旋转不同主坐标轴，不能只盯着PC1/PC2，有些数据在PC3/PC4上分离得更明显；第三，配合PERMANOVA检验，也就是beta-group-significance输出的p-value，p<0.05才谈得上组间显著差异。如果PCoA图上看着有分离趋势但p值不显著，可能是样本量太少或组内异质性太大，不能简单下结论说“没有差异”。

5. 分类学注释：数据库选择与分类器训练的取舍

5.1 三个主流数据库怎么选

分类学注释就是把ASV代表序列与参考数据库比对，猜出它们属于哪个门纲目科属种。目前最常用的三个16S数据库：

我个人的建议：如果课题组之前大量使用了Greengenes，为了批阅对比一致性继续用Greengenes也说得通，但新项目建议直接用SILVA或GTDB。尤其临床样本或环境样本中有不少未培养类群，SILVA在这块的注释能力更完备。

5.2 预训练分类器 vs 自己训练

QIIME2官方提供了基于Greengenes的预训练分类器（515F/806R区域），但针对SILVA或GTDB，预训练分类器有时不是专门针对你引物区域训练的，这时我更建议自己训练。过程分两步：

第一步，按引物区域从参考序列里提取片段：

qiime feature-classifier extract-reads \ --i-sequences silva-138-99-seqs.qza \ --p-f-primer GTGCCAGCMGCCGCGGTAA \ --p-r-primer GGACTACHVGGGTWTCTAAT \ --p-trunc-len 300 \ --o-reads ref-seqs.qza

第二步，训练朴素贝叶斯分类器：

qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier silva-138-515-806-classifier.qza

自己训练的耗时通常几十分钟到一两个小时，完全可以接受。训练好的分类器要妥善保存，同一个引物区域的多个项目都能复用。如果你用的是公司交付的通用V3-V4区引物（比如341F/806R），提取片段时一定换成对应的正向和反向引物序列，不要拿V4区的515F/806R硬套。

5.3 注释结果的可视化与分层柱状图

分类注释的核心命令：

qiime feature-classifier classify-sklearn \ --i-classifier silva-138-515-806-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv

注释结果里你会看到每个ASV对应的一系列分类层级，以及置信度。实操里我有一个检查习惯：看属水平上有没有“未分类”的比例过高的样本，如果某个样本Unclassified占比超过30%，先不要急着怀疑注释，回去看看该样本reads数是不是太少，或者质控后留下的序列是不是质量偏低。

从这里还能生成分类学堆叠柱状图：

qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv

在View里可以切换不同分类层级（门/纲/目/科/属），按分组着色查看相对丰度结构。我一般会先用“门水平”看宏观趋势，再用“属水平”看具体差异类群。

6. 下游差异分析与R生态衔接：一条龙的最后一段路程

6.1 ANCOM与ANCOM-BC的选型逻辑

差异丰度分析用来回答“哪些物种在组间显著差异”。QIIME2内置了ANCOM，命令是：

qiime composition ancom \ --i-table table.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization ancom.qzv

注意ANCOM有严格的内部假设：它假设绝大多数分类单元在组间没有差异（通常少于25%的OTU发生差异）。如果实际数据里差异物种占比很高，ANCOM的结果会偏保守，很多真正有差异的物种被检测不出来。因此如果样本量允许，我更推荐用较新版本的QIIME2中提供的q2-composition扩展或直接在R里跑ANCOM-BC。ANCOM-BC加入了偏差校正，可以输出完整的p值和置信区间，处理稀疏数据的能力也更强。

实操时我还有一个提醒：差异分析之前一定要把特征表里丰度太低、出现于太少样本的ASV过滤掉。我个人常用过滤规则：至少在某一个分组内超过20%的样本中出现，且总丰度大于等于5。过滤命令：

qiime feature-table filter-features \ --i-table table.qza \ --p-min-frequency 5 \ --p-min-samples 2 \ --o-filtered-table filtered-table.qza

不事先过滤，少数ASV的极端丰度变化会干扰整体分析的稳健性。

6.2 功能预测PICRUSt2：注意它不在内置插件里

16S本身只提供物种组成信息，功能预测则是利用参考基因组推断群落可能具备的代谢通路。PICRUSt2是目前的主流工具。需要注意，PICRUSt2已经不再作为QIIME2内置插件发布，需单独安装并作为独立命令行使用：

conda install -c bioconda -c conda-forge picrust2

使用流程大致是把QIIME2的table和rep-seqs导出为PICRUSt2能接受的BIOM格式和fasta格式：

qiime tools export --input-path table.qza --output-path picrust-input qiime tools export --input-path rep-seqs.qza --output-path picrust-input biom convert -i picrust-input/feature-table.biom \ -o picrust-input/feature-table.tsv --to-tsv

然后运行PICRUSt2的一个核心命令：

picrust2_pipeline.py \ -s picrust-input/dna-sequences.fna \ -i picrust-input/feature-table.biom \ -o picrust-output \ -p 4

生成的预测结果包含KEGG Orthologs、MetaCyc pathway等。做这一步前务必想清楚：16S的功能预测是基于系统发育推断的“潜在功能”，不是实测的宏基因组结果。你可以在讨论里描述为“群落代谢潜能的初步推测”，不能直接当宏转录组或宏基因组结论使用。

6.3 把数据导出给R生态：phyloseq与ggplot2的衔接

不少下游可视化最终要离开QIIME2界面去R里完成，比如个性化的PCoA图、热图或LEfSe分析。最常用的做法是导出BIOM和分类学表：

qiime tools export --input-path filtered-table.qza --output-path export/table qiime tools export --input-path taxonomy.qza --output-path export/taxonomy qiime tools export --input-path rooted-tree.qza --output-path export/tree

导出后的feature-table.biom用R的phyloseq读取：

library(phyloseq) ps <- import_biom("export/table/feature-table.biom") tax <- read.delim("export/taxonomy/taxonomy.tsv", row.names = 1)

然后再做alpha多样性箱线图、NMDS/PCA散点图都相对自由。这里我想额外提一句：在QIIME2 View里直接导出的PCoA图适合快速看趋势，但论文投稿建议还是用R重新绘制，ggplot2对点的大小、颜色、标签、置信椭圆的控制要灵活得多。

6.4 LEfSe之类工具的定位

很多教程会把LEfSe（线性判别分析效应量）也塞进QIIME2流程，但严格说LEfSe是独立工具，通常需要从Galaxy或本地运行。它适合筛选组间生物学差异的标记类群，输出LDA得分与柱状图。实际项目中，我会先用ANCOM或ANCOM-BC做初步筛选，再用LEfSe丰富结果的可视化。如果你有分组且样本量足够，LEfSe对故事线的构建很有帮助，但它同样对样本量敏感，样本太小时结果容易被个别极端值驱动。

说到最后，我还是想强调“一条龙”并不等于“一路默认参数跑到底”。从导入数据开始，每一步的取舍实际上都在影响最后的生物学结论。最典型的就是DADA2截断长度和抽平深度选择，这两个参数在不同项目里必须分别判断。作为从业者，我建议你每个环节都留下日志、保存好qzv可视化文件，并在Methods里准确记录每个参数。这样就算半年后审稿人要求补分析，你也能迅速复现整个流程，而不是对着文件夹里几百个.qza发呆。

如果在实战中遇到导入报错、注释结果异常或者多样性图与预期不符，欢迎在评论区带上你的具体报错信息或qzv截图来讨论，我看到后会结合自己的项目经验帮你一起排查。16S分析这条路上，每一个参数背后都是可积累的细节，愿你少走弯路。