16S扩增子属水平分析完整流程：从数据质控到注释与可视化

一测序数据拿到手，报告里从门到属列了一大串名字，真正想回答的问题——比如"这个菌属到底在两组之间有没有差异""哪种属撑起了样本之间的差异"——却往往卡在注释这一步。做微生物组分析的人都会经历这个过程：从原始下机数据一路跑到属水平结果，中间任何一步不小心都会让前面所有工作白费。这篇文章就是把"genus综合流程"从质控、注释到下游可视化完整串一遍，适合刚入门生信的研究生、检测机构里天天跑流程的实验员，也适合那些自学了很久但总在某个环节踩坑的人。

1. 为什么要卡在"属"这一级：分类分辨率的现实选择

1.1 从界门纲目科属种说起

生物分类系统从界到种有七个层级，16S扩增子测序在门、纲、目、科这几个层级上通常都能给出比较稳定的注释结果，但到了属这个级别，问题就开始有意思了。属是"能够被现有数据库相对可靠地注释、又能够保留足够生物学差异信息"的最小分类单元。换句话说，属水平分析是分辨率和可靠性的一个平衡点。

很多初学者上来就问："能不能直接跑到种？"理论上可以，但16S的保守区域（比如V4区或者V3-V4区）序列长度只有250~500bp，信息量有限，很多亲缘关系很近的种在16S序列上根本分不开。强行注释到种，会得到大量"未分类种"或者不准确的结果。而如果真的需要种水平分辨率，应该考虑全长度16S或宏基因组测序，而不是指望扩增子数据分析去硬撑。

在实际项目中，属水平分析能直接回答三类问题：群落中有哪些核心属、不同样本之间优势属和稀有属的结构差异、以及特定条件（疾病、处理组）下有哪些属发生了显著丰度变化。绝大多数已发表的微生物组研究，核心结论也是落在属这一层。

1.2 属水平分析能回答什么问题

举个例子。做肠道菌群与2型糖尿病关联研究时，你可能会发现厚壁菌门和拟杆菌门的丰度比例（F/B ratio）在两组之间有差异——但门的层级太粗了，门水平的变化往往掩盖了门内不同属的此消彼长。可能厚壁菌门整体没有变化，但里面的梭菌属（Clostridium）多了、乳杆菌属（Lactobacillus）少了。这种"内部置换"是微生物组的常见现象，只有到属水平才能看出来。

再做临床样本分析，医生最关心的是"哪个属和疾病最相关"。如果只停在科甚至门的层级，给出的结论会让临床合作者觉得太笼统。属水平可以落到具体的益生菌或致病菌候选属，后续支持动物实验或者临床验证时才会有方向。

另外，像普雷沃氏菌属（Prevotella）和拟杆菌属（Bacteroides）的比值，在很多营养学和代谢病研究中是直接拿来当作"肠型"判别的指标。这些指标从定义上就要求属水平分类注释。

1.3 什么时候不要死磕属水平

属水平分析也有它的局限。第一，如果参考数据库里某些属的16S序列本身很少，注释器给出的置信度很低，这时候硬要"属"的结果就会引入噪声。第二，在宏基因组数据里如果用了属水平的标记基因分析，基因拷贝数差异会造成丰度的系统偏差，直接比较属间丰度会失真。第三，如果样本采集或DNA提取环节引入了污染，属水平会把污染源（比如水里的假单胞菌属）也识别出来，分析解读时要格外小心。

所以我的习惯建议是：正式分析以属水平为主，但必须同时检查门、纲、科级别的结果是否与属水平一致。如果某个属的变化在更高分类层级中完全看不到——甚至与门级趋势相反——就要重新审视这个属的注释质量。

2. 流程总览和软硬件准备

2.1 一条主线：从原始数据到属水平丰度表

完整的genus综合流程可以概括为八个环节：原始数据质检、引物切除与质控过滤、合并双端序列、去噪聚类生成特征表、代表性序列分类注释、生成属水平丰度表、下游差异分析与可视化、结果解读验证。

每一环的输出都是下一环的输入。这个链条中，我特别想强调一个容易被忽略的问题：你选择的不同参数会沿着链条传导下去。比如质控阶段切掉太长或者太短的序列，直接影响后续特征表里的真实成员；参考数据库版本变了，注释结果就可能在同一份数据上出现明显差异。

对刚接触的朋友，我建议不要在分析刚开始就把所有"高级选项"全部打开，先拿默认参数跑通整个流程，再用少量样本做参数敏感性测试。等流程稳定了再铺开全部样本。

2.2 软件选型：QIIME2、DADA2还是自写pipeline

目前主流方案有三种。

第一种是QIIME2完整流程，把质控、去噪、注释、多样性分析全部集成在一个环境里，对新手最友好，操作也比较标准化。但它也有问题：QIIME2的版本更新很快，不同版本之间的数据格式（尤其是FeatureTable和分类器）并不完全兼容，而且它的中间产物格式比较特殊，不喜欢跟外部工具互通的人可能会觉得别扭。

第二种是R语言跑DADA2全流程。DADA2的优势在于R生态特别丰富，从DADA2到phyloseq再到ggplot2可视化，一条链下来全在R里完成。而且DADA2去噪得到了ASV（Amplicon Sequence Variant）而不是传统的OTU，从原理上说更能反映真实的生物学序列差异。缺点就是数据库格式需要自己准备，训练分类器也是个额外步骤，对不熟悉R的人有门槛。

第三种是混合方案：用fastp或cutadapt做质控、用usearch或vsearch做聚类、再用自己训练的分类器注释。这种方式灵活性最高，能针对特殊引物或特殊样本类型做精细调整，但所有环节都得自己组装，出错概率也高。

我的个人建议是：如果你从零开始，选QIIME2或DADA2二选一，先跑通再考虑优化。如果你有大量样本而且需要完全可控的分析，选混合方案。

2.3 参考数据库怎么选：Greengenes2、SILVA、RDP

数据库是注释的"词典"，没有匹配的序列就永远注释不出来。目前常用的三个是Greengenes（经典版和2024年新版Greengenes2）、SILVA（132/138版本）和RDP。

三个数据库各有特点。Greengenes在16S分析里用得最久，很多历史文献的OTU表格也是按Greengenes注释的，如果你要对标旧文献结果，用Greengenes比较方便。SILVA的覆盖率更高，覆盖细菌、古菌和真核微生物，而且它的分类体系相对稳定，是目前扩增子分析的首选数据库之一。RDP提供的分类器训练接口比较友好，但整体覆盖率不如前两者。

需要特别提醒的是：数据库的选择会直接影响属水平注释率。我做土壤样本时用过Greengenes和SILVA对比，同样的数据，SILVA 138版本注释到属的比例比Greengenes 13_8高大概5到8个百分点。所以在正式跑全量分析之前，先拿一组代表性样本做注释率对比实验，选注释率较高的那个数据库作为正式分析基础。

下表是我在实践中总结的数据库选择参考：

3. 从下机数据到OTU/ASV表的关键步骤

3.1 质量控制与引物切除

拿到下机数据的第一步永远是看质量。用fastqc跑一批样本，看每个位置的碱基质量得分（Q score），重点关注两条信息：序列长度的分布和末端质量是否急剧下降。如果用的是双端250bp测序，而你的插入片段长度在460bp左右，拼接后刚好能覆盖V3-V4区域，这个设计是最理想的；如果插入片段太短，两端的引物序列会占据很大比例，直接导致有效序列变少。

引物切除我习惯用cutadapt，切完引物之后再用trimmomatic或fastp做质量修剪。这里有个细节容易被忽略：切除引物必须在质控之前还是之后？我的顺序是先切引物再做质量修剪。原因是如果不先把引物切掉，质量修剪可能会把引物所在的位置当作低质量区域直接裁掉，反而把下游有效序列切短。

设计引物切分时要注意允许一定比例的错配。cutadapt示例参数如下：

cutadapt -g ^GTGYCAGCMGCCGCGGTAA -G ^GGACTACNVGGGTWTCTAAT \ --discard-untrimmed -e 0.2 -m 200 -M 480 \ -o R1_trimmed.fastq.gz -p R2_trimmed.fastq.gz \ raw_R1.fastq.gz raw_R2.fastq.gz

这个命令的含义是：正向引物和反向引物必须在序列开头（^表示锚定在5'端），允许20%的碱基错配率，丢弃没能匹配引物的序列（--discard-untrimmed），且保留的序列长度范围在200~480bp。这些参数在正式上生产流程之前，应该拿少量样本试跑并检查结果。

3.2 去噪与聚类：ASV和OTU怎么选

QC完成后，下一步就是把质量合格的序列聚成特征单元。传统方法是按97%相似度聚成OTU，现在的DADA2方法是经过去噪（denoising）后直接得到ASV，ASV并不依赖于任意设定的相似度阈值，理论上每个ASV只代表一个真实的生物学序列变体。

在实际操作上，DADA2去噪有两个关键参数：maxEE（期望错误数）和truncLen（截断长度）。maxEE=2是常用的默认值，对于双端250bp的数据，V4区域的truncLen我一般设置成240和200，也就是每个读段在末端低质量区域被截掉10~50bp。这两个参数的组合需要根据你的数据质量去调，没有万能答案。一个简单原则：如果末尾20bp的质量图掉到了Q20以下，就应该适当增大截断长度。

library(dada2) fnFs <- snakemake@input$fq1 # 示例占位，实际用你的文件路径 fnRs <- snakemake@input$fq2 filtFs <- file.path("filtered", paste0(sample_names, "_F_filt.fastq.gz")) filtRs <- file.path("filtered", paste0(sample_names, "_R_filt.fastq.gz")) out <- filterAndTrim(fnFs, filtFs, fnRs, filtRs, truncLen=c(240,200), maxN=0, maxEE=c(2,2), truncQ=2, rm.phix=TRUE, compress=TRUE, multithread=TRUE) errF <- learnErrors(filtFs, multithread=TRUE) errR <- learnErrors(filtRs, multithread=TRUE) dadaFs <- dada(filtFs, err=errF, pool="pseudo", multithread=TRUE) dadaRs <- dada(filtRs, err=errR, pool="pseudo", multithread=TRUE)

pool="pseudo"这个参数在多批次数据合并分析时特别有用，它能在不显著增加计算量记忆量的前提下，让来自不同样本的低丰度序列仍有被识别为同一变体的机会。

ASV的优势是分辨率高、可重复性好，不同研究之间的ASV可以直接比。但它也有代价：ASV数量通常比97% OTU多很多，在属水平注释时会产生大量低频ASV，这些低频ASV往往连属都注释不上，就变成了"数据里的黑洞"。实践上我的处理策略是：正式分析以ASV为主，但在做属水平差异分析前，先按属水平合并，并过滤掉那些总丰度非常低（比如在全部样本中相对丰度总和低于0.01%）的ASV。

3.3 生成特征表和代表序列

无论走DADA2还是vsearch聚类路线，最后都会生成一张特征表（Feature Table，OTU/ASV × 样本的计数矩阵）和代表序列（每个OTU/ASV的一条代表性序列）。这张特征表就是后续所有分析的输入基础。

生成特征表时要记录清楚每个步骤过滤了多少条序列。我见过不少人最后只汇报"最终得到N条有效序列"，但是如果中间QC步骤把80%的数据都扔掉了，却没有记录下来，等于把质量问题藏住了。建议每个环节都保留日志，从原始序列开始，每一步记录序列数变化，形成一个清晰的"数据漏斗"。这不仅是为了自己检查，投稿时审稿人问起数据质量时，这个日志就是最有说服力的材料。

4. 属水平注释的实现与坑

4.1 训练朴素贝叶斯分类器

注释这一步是整个流程中最影响最终结果质量的地方。QIIME2内置了对Greengenes和SILVA的预训练分类器，可以直接下载使用，但如果你用的是自己设计的目标区域引物（比如只扩增ITS2或者特定的功能基因），就必须自己训练分类器。

训练分类器的原理是：用参考数据库里的序列和对应的分类学信息，训练一个朴素贝叶斯分类器，然后对每个ASV/OTU的代表序列做分类。QIIME2的q2-feature-classifier模块是主流选择。

训练分类器时有一个细节：参考序列应该按你的目标区域来截取，而不是直接用全长16S序列。原因是分类器学习的是目标区域内的序列特征，如果你用全长序列训练，而待分类的序列只是V3-V4区域，两者特征分布不匹配，分类准确率会下降。

用SILVA 138.1训练V3-V4区域分类器的示例：

qiime feature-classifier extract-reads \ --i-sequences silva-138.1-ssu-nr99-seqs.qza \ --p-f-primer ACTCCTACGGGAGGCAGCA \ --p-r-primer GGACTACHVGGGTWTCTAAT \ --p-trunc-len 469 \ --o-reads ref_seqs_v34.qza qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref_seqs_v34.qza \ --i-reference-taxonomy silva-138.1-ssu-nr99-tax.qza \ --o-classifier silva_v138_v34_classifier.qza

这个过程在全长序列上可能要跑一两个小时，但在目标区域截取之后的数据集上通常很快。如果样本是特殊环境（如高盐、热泉），建议在训练集中额外加入环境中新测序获得的序列，能明显提升这些陌生序列的注释率。

4.2 注释置信度阈值怎么定

分类器输出每个分类层级时都会附带一个置信度（confidence），典型范围是0到1。默认的置信度阈值是0.7，意味着当分类器判定"这个序列属于某个属"的置信度低于0.7时，结果会被标为未分类。

我实测过不同阈值的效果：阈值设到0.8以上，注释到属的序列比例下降10%左右，但留下来的都是高置信度结果；阈值降到0.5，注释率上去了很多，但明显会出现错分——比如把埃希菌属（Escherichia）跟志贺菌属（Shigella）相互掺混，这种错误在16S短片段上本来就难以避免，降低阈值等于给这类错误开了后门。

我的经验是：环境样本（土壤、水体）里有很多参考数据库里没有的序列，阈值保持0.7合适；临床肠道样本中已知菌属较多，可以适当提高到0.8。所谓"提高阈值提高准确性"并不是绝对的，因为阈值越高，留给你做差异分析的"有效属"就越少，统计功效会受影响。

4.3 一个经常被忽略的问题：多属注释

当我们说"这个ASV注释到属X"，默认是说它的最佳匹配是属X。但现实中经常出现一个ASV的序列既与属X匹配、又与属Y匹配得分相同或相近，尤其当两个属在16S目标区域完全相似时。这种情况下，分类器给出的单一最佳匹配就变得不可靠。

我在分析时习惯统计注释结果中那些"ambiguous"（模棱两可）的条目。如果这一类条目占比超过5%，就说明当前的目标区域或数据库确实无法支撑属水平的判断，需要重新考虑是否把分析层级降低到科，或者在该区域引入更加特异的参考序列。

处理多属注释的常见做法是设置confidences-based filtering，在给属水平丰度表算占比时，把所有质量不够高的ASV归入"未注释属"这一类，而不是强行分配。这样虽然会丢掉部分信息量，至少在结论层面是诚实的。

5. 下游分析：属水平有哪些能做的

5.1 从特征表到属水平丰度表

拿到ASV/OTU特征表之后，第一个动作就是把每个特征映射到属，并按属水平合并。这一步在QIIME2里叫taxonomic collapse，在R里用phyloseq的tax_glom函数就能完成。

library(phyloseq) ps <- readRDS("ps_ASV.rds") # 你的phyloseq对象 ps_genus <- tax_glom(ps, taxrank="Genus", NArm=FALSE)

注意这个NArm参数。默认是TRUE，会把那些注释不到属的ASV直接丢弃，但如果你关心的是"注释率到底有多高"，建议先设成FALSE，看看未注释的部分占了多少比例，再决定要不要过滤。

合并之后还有一个关键步骤是归一化。最常用的是总丰度归一化（Total Sum Scaling，TSS），也就是把每个样本里的所有属的计数转换为相对丰度（百分比）。但要注意：如果样本测序深度差异很大，TSS归一化后低深度样本的相对丰度会不稳定——就像抽样件数太少时比例估算的波动会很大。更稳健的做法是使用CSS（Cumulative Sum Scaling）或rarefaction（抽平），omegaQA。于比较简单和主流的表达习惯，我一般先给出TSS归一化结果，再针对关键属做均值-方差检查，如果方差随均值变化明显，就用DESeq2或ANCOM-BC的模型内置归一化方式重新算。

5.2 堆叠柱状图与热图：属水平可视化的两种常见形式

属水平结果最经典的呈现是堆叠柱状图，每个样本一根柱子，不同颜色代表不同属，柱子的高度表示相对丰度。制作时有一个容易踩的坑：当属的数量超过20个时，图右的图例就会变成一堆密密麻麻的小色块，根本读不了。实际建议是把丰度低于某个阈值的所有属合并成"Others（其余）"，通常以5%或2%为界。

另一个更好用的图是分组堆叠柱状图，每个处理组用多根柱子展示，相当于把每个组的属水平组成并列对比。用ggplot2画这种图时，排序很重要——如果按属名首字母排序，图会非常杂乱；建议先对属按照整体平均丰度排序，或者按照主成分分析排序，让视觉上最明显的变化方向（如对照组到处理组）自然凸显出来。

热图适合展示样本和属的关系。行是属，列是样本，颜色代表相对丰度。用pheatmap画热图时，我建议将丰度数据做log10(x+1)转换，否则丰度差异大的少数属会压掉所有其他属的颜色梯度，极差大看得很累。在热图旁边一定要标注聚类树，并且标注属的分类学信息——至少标注门或纲，这样读者才能快速理解"哪些门内部的属聚类在一起"。

5.3 差异分析：LEfSe与ANCOM-BC

找到"哪些属在组间显著变化"是属水平分析的核心输出。最常用的工具有LEfSe、DESeq2、edgeR、ANCOM-BC。

LEfSe是经典的LDA Effect Size方法，它在微生物组文献里用得非常多，适合两组或多组比较。它的漂亮之处是先做非参数Kruskal-Wallis检验，再做LDA评分，把"有显著差异"和"效应大小"两个维度一起呈现出来。实际运行时的阈值通常是LDA score > 2.0且p < 0.05。LEfSe的输出图（条形图加 cladogram）很直观，投稿时很多期刊都喜欢。

DESeq2和edgeR本来是从转录组学引入的方法，它们在处理高维稀疏数据时有统计模型支撑，性能比简单的非参数检验更能处理小样本。但它们在微生物组数据上有过拟合风险，尤其当样本量很少（每组少于5个）时容易把不太稳的差异也筛出来。

ANCOM-BC是最近几年比较被看好的方法，它显式地考虑了微生物组数据的组成性质，并通过偏置校正来部分缓解相对丰度分析中"此消彼长"带来的伪差异。如果你的样本量足够（每组10个以上），我建议直接上ANCOM-BC，结果的可信度会更高。

实际操作中我建议至少跑两种方法，比如ANCOM-BC + LEfSe，取交集；只被一种方法检出的属，在结论里标注"某一种方法支持"即可，不要当作确凿结论。

5.4 多样性分析注意分辨率问题

属水平分析同样可以算alpha多样性（如Shannon指数、Chao1指数）和beta多样性（如Bray-Curtis距离、UniFrac距离）。

alpha多样性方面，属水平会损失一部分信息量，因为很多稀有ASV在属水平合并后看不出多样性贡献。一般做法是在ASV水平算alpha多样性，再考虑属水平单独展示。beta多样性方面，属水平矩阵的PCoA图也常见发表，但要注意：如果属数量太少（比如某些极端的单菌属样本），PCoA结果的解释力会明显下降。

一个实用建议：在做PCoA时，同时算weighted UniFrac和unweighted UniFrac，两者结果可能相差很大——这本身就是一个重要信息，说明样本间的差异到底是来自优势属的丰度变化，还是来自稀有属的有无变化。这两种生物学解读在结论上可能完全不同。

6. 常见问题与排查

6.1 属水平注释率很低怎么办

注释率低有几个常见原因。一是参考数据库覆盖不足，尤其在土壤、深海或特殊宿主样本中；二是目标区域本身保守，不同属之间差异太小；三是质控时过滤太严格，剩下有效序列长度不足，分类信息被切掉了。

遇到注释率低，第一步检查数据漏斗，看看是不是在质控阶段把大量序引物序列被误认为接头序列处理掉了。第二步换个数据库试试，Greengenes2有时能力挽狂澜。第三步如果还是不行，考虑降低分类层级到科，或者去更针对性的数据库做BLAST验证，而不是强行用低置信度的属注释结果。

6.2 样本量小但属数量太多

这是高频问题：样本每组只有3~5个，但属水平特征有100多个，直接做差异分析会产生多重检验问题，假阳性率很高。解决方案有两条路。一条是预先过滤掉丰度过低、出现频率过低的属。比如设定"在超过20%的样本中相对丰度都低于0.1%的属直接过滤"。另一条是采用多重检验校正，比如BH-FDR校正，并把显著性阈值适当放宽。两者要结合使用，单纯靠p.adjust在样本量很小时效果有限。

6.3 同一属内不同的种差异很大

有时候属水平分析会掩盖真实信号。比如拟杆菌属（Bacteroides）里，脆弱拟杆菌（B. fragilis）和多形拟杆菌（B. thetaiotaomicron）的功能完全不同。如果你只在属水平观察到拟杆菌属无差异，但临床特征有明显表现，要怀疑是不是种水平发生了置换。

此时有两个选择：一是把关注的属单独提取出来，调出该属内所有ASV的注释结果，做种水平的展示和检验；二是重新用更高的分辨率工具分析原始数据。很多研究中"属水平无差异但种水平有差异"就是这样被发现的。

6.4 参考数据库版本不同导致结果对不上

同一个样本用Greengenes 13_8和SILVA 138注释，结果对不上的概率极大。不是因为你分析错了，而是因为数据库的分类学修订本身在变。如果你的研究需要与历史文献对比，建议用相同的数据库版本和相同分类器，才能做真正的比较。

如果已经用不同数据库跑了，要对比结果，最好的办法是只对比有"同义名"关系的属，比如查一下SILVA里的"Clostridium sensu stricto 1"在旧Greengenes里是否标注为"Clostridium"。这些命名差异很容易让人误以为结果不一致，但其实只是数据库重命名的问题。

6.5 一个容易被忽略的额外问题：基因拷贝数

16S rRNA基因在细菌基因组中不总是一个拷贝，有些菌属（比如芽孢杆菌属）可以有10个以上拷贝，这会让基于16S的相对丰度在属之间产生偏差。对于肠道菌群领域，通常不太做拷贝数校正，因为校正需要种水平的拷贝数信息，而这些信息在扩增子数据里本来就难以获得。但如果你做的样本以芽孢杆菌属等高频多拷贝菌为主，建议到rrnDB数据库查一下相关属的拷贝数中位数，在结果解读时说明潜在偏差方向。

我在实际项目中遇到过这类情况：先前某环境样本中芽孢杆菌属的相对丰度是30%，拷贝数校正后可能只有8%，这直接影响到了"优势属"的判定。如果读者对这类偏差敏感，建议方法学部分明确写明"未进行16S拷贝数校正，结果以相对丰度呈现"。

7. 实操经验与内容扩展方向

先分享几个我踩过不止一次坑的小技巧。

第一，流程跑完后一定要回到原始数据抽样检查几样本。我当时遇到过某个样本因为PCR扩增效率低，短片段占比较高，与正常样本明显不同，但流程本身并不会报错。这种数据混在分析里，属水平结果就是解读不出规律。手动抽看原始测序质量能救回一批"伪坏样本"。

第二，别把样本全放到一个筐里。同一研究里如果包含不同批次、不同提取方法、不同建库方式的样本，这些技术差异会在属水平形成很强的批次效应。如果条件允许，把样本按照实验批次单独跑一遍流程，比较各批次中"稳定属"与"差异属"，再决定是否合并分析。

第三，属水平不是终点，可能是起跳板。在很多情况下，属水平分析发现的关键菌属，后续可以用实时荧光定量PCR或宏基因组测序做靶向验证。这是我目前项目里常用的做法：先用扩增子测序找到候选属，再用qPCR对特定属进行绝对定量验证，投稿时审稿人对这套"相对丰度+绝对定量"组合的认可度远高于单独一个相对丰度结果。

最后再分享一个小技巧：即便是属水平的结果，一定不要忘记保存每个ASV的代表序列和注释信息。很多人在完成分析后只保留一个"属×样本"的矩阵，后续如果想补充种水平分析、做进化树或者提交序列到公共数据库，才发现原始的代表序列文件早就丢了。保存好这些中间产物，对复现和深入分析都是最简单的保障。