做微生物组数据分析的人，十有八九都经历过这样的场景：QIIME2或者DADA2跑完，满怀期待地打开特征表，结果屏幕上躺着一万多行OTU/ASV。这其中有大量只出现过一条reads的幽灵序列，有不知道怎么就混进来的叶绿体和线粒体，还有每个样本之间相差悬殊的测序深度。第一次拿到这种数据的人，往往在“过滤”和“相对丰度转换”这两个看似简单的环节上卡住，不知道该删哪些、阈值定多少、标准化到底怎么选。这篇文章就把这条从OTU/ASV过滤到相对丰度转换的标准化流程完整拆解一遍，把每一步的取舍逻辑和坑都讲清楚，适合刚入门微生物组研究的学生，也适合想把手头多样本数据重新规范处理的研究者。整个流程做下来，你会发现后续的多样性分析和差异比较都稳得多。
1. 项目概述：这套流程解决了什么问题
1.1 从OTU到ASV：特征表的进化史
先把这个最基础的概念捋清楚，因为后面所有过滤操作都建立在“你手里这张表是怎么来的”之上。OTU，也就是操作分类单元（Operational Taxonomic Unit），经典做法是拿所有序列按97%相似度聚类，然后把每个簇当作一个“物种类群”。这种做法的初衷是容忍测序错误——反正PCR和测序都会引入噪声，聚一下类可以把单个碱基错误造成的假序列合并掉。代价也很明显：聚类阈值是人为定的，不同研究、不同聚类算法之间，OTU的边界并不完全可比，你今天聚出的OTU1和隔壁实验室聚出的OTU1，可能根本不是同一个东西。
ASV是这几年更主流的做法，全称是扩增子序列变体（Amplicon Sequence Variant），由DADA2、Deblur这类去噪算法得到，理论上能达到单碱基分辨率。它不依赖97%聚类这种人为设定，能把样本里真实存在的序列变异直接分出来，而且可重复性很高，哪怕是不同研究的数据，只要都用同样的去噪流程，特征是可以直接对齐比较的。所以在条件允许的情况下，我一般建议优先用ASV流程。
但实际项目里，很多存量数据还是OTU格式，或者别人给了你一张已经生成好的OTU表，这个时候没必要推倒重来。好在从数据分析的角度看，OTU表和ASV表在phyloseq里结构完全一样，都是“特征×样本”的计数矩阵加分类学注释表，下面这套标准化流程两种表都能直接用。
1.2 标准化流程的三个关键节点
整套流程浓缩起来就是三件事：去掉污染序列、去掉低质量低频特征、把计数矩阵标准化到相对丰度。三者缺一不可。
污染序列不去掉，你分析出来的可能根本不是微生物群落，而是叶绿体和线粒体的狂欢。在植物根际土或者叶片样本里，叶绿体序列占比超过一半都很常见，不处理的话，后续的多样性指数、PCoA排序全都会被带偏。低质量低频特征不去掉，几万个只有一条reads的ASV会让alpha多样性虚高，也会让beta多样性里的距离计算被大量噪声主导。标准化不做，你没法公平比较不同样本之间“同样一个物种”有多少，因为每个样本测序深度根本不一样，绝对reads数不能直接比较。
这三步做完，得到的那张“相对丰度表”，才是后续做PCA、PCoA、热图、LEfSe这些分析的正确输入。
2. 数据准备与过滤策略设计
2.1 拿到数据先别急着过滤：先读懂特征表
很多同学拿到BIOM文件或者TSV格式的OTU表，直接就套用网上的过滤代码，结果要么过滤得太狠，一整套分析下来只剩几百个特征，要么过滤了个寂寞，跑出来的图和没过滤一个样。问题就是没有先花十分钟把数据结构看清楚。
我建议拿到数据后，第一件事是用R读进来，确认三样东西：特征表的维度、样本测序深度的分布、分类学注释各层级的情况。这一步不是走过场，它能直接告诉你后面过滤参数该怎么定。
library(phyloseq) library(tidyverse) # 如果你有BIOM文件 biom <- import_biom("feature-table.biom") # 样本元数据 meta <- read.delim("metadata.txt", row.names = 1) ps <- merge_phyloseq(biom, sample_data(meta)) # 三维体检 ps # 看特征数和样本数 sample_sums(ps) %>% summary() # 看每个样本的reads总数分布 View(as.data.frame(tax_table(ps)) %>% head(20)) # 看分类注释长什么样一个比较理想的状态是：样本数几十个，特征数几千到几万，每个样本的reads数在一万到二十万之间。如果你发现某个样本只有几百条reads，大概率是建库或测序出了问题，这种样本在过滤之前就应该考虑是否剔除，而不是指望靠过滤来救。
2.2 污染序列过滤要优先处理
污染序列是数据里最“硬”的噪声，因为它们来自明确的生物学来源，只是不是你关心的微生物。最常见的两类：
第一类是叶绿体序列。16S扩增子所用的引物本来就是针对细菌16S rRNA基因设计的，但叶绿体16S rRNA基因和细菌同源，植物组织残渣一旦进入样本，就会被一起扩增出来。在分类注释结果里，这类序列在Order层级经常被注释为Chloroplast。第二类是线粒体序列，动物线粒体16S和12S区域也很容易被扩增，在Family层级常见注释为Mitochondria。
我在处理水稻根系样本时见过最夸张的情况，叶绿体reads占比超过60%。如果不做过滤，整个样本的群落结构展示出来基本就是“叶绿体丰度排名第一”，没有任何生物学意义。
另外还有一种容易被忽略的污染是宿主序列。动物样本里如果混入宿主DNA，即使不是16S区域，有时也会被非特异性扩增。这类污染通常无法通过分类学注释直接识别，需要比对参考基因组来排除。在流程里我们主要处理前两类，宿主污染一般依靠建库前的去除步骤来控制。
2.3 低丰度过滤的两个维度：频率和丰度
低丰度过滤是整个流程里最需要“手感”的环节。常见的错误是只用单一维度去切，比如只看总reads数小于10就删，结果忽略了一个关键问题：一个特征在100个样本里有90个都检测到，只是每个样本里都很少，那它可能是一个普遍存在的低丰度共生菌，直接删掉就可惜了。
所以我一般把过滤拆成两个维度：
频率过滤（Prevalence filter）。统计这个特征在多少个样本中出现过（reads>0）。如果只在极少数样本里出现，那它大概率是扩增错误、测序错误或者样本间交叉污染。常用经验阈值是“至少在10%的样本中出现”，样本量比较大时可以提高到20%，样本量很小（比如少于20个样本）时这个阈值基本等于放行所有特征，作用不大。
丰度过滤（Abundance filter）。统计这个特征在所有样本中的总reads数或平均reads数。DADA2输出里有一大堆singleton——整个数据集中只出现过一条reads的ASV，这些基本可以确定是测序噪声，直接删掉不会丢失真实的生物学信号。常用阈值是总reads数大于等于10，或者相对丰度大于0.01%。
两个维度要配合使用，不要只选一个。我的习惯是先在频率维度上做一个较温和的过滤，保证常见的低丰度共生菌能留下来，再做一次丰度过滤，把那些总reads极低的噪声序列清掉。
3. 过滤实操：基于phyloseq的完整流程
3.1 数据导入与预检查
上一步已经提到导入了，这里直接进入预检查后的实操。拿到phyloseq对象之后，强烈建议先输出几个关键数字，作为后续过滤效果的“基线”。
# 过滤前基线 baseline_taxa <- ntaxa(ps) baseline_reads <- sum(sample_sums(ps)) # 看看每个样本的reads分布 sink("filter_log.txt") # 建议开一个日志文件 cat("=== 过滤前 ===\n", file = "filter_log.txt") cat("特征数:", baseline_taxa, "\n", file = "filter_log.txt") cat("总reads:", baseline_reads, "\n", file = "filter_log.txt") cat("样本reads范围:", range(sample_sums(ps)), "\n", file = "filter_log.txt") sink()这一步的核心目的是留痕。等文章写好后，方法学部分写“过滤条件：去除仅在少于10%样本中检出的特征；去除总reads数小于10的特征”，审稿人如果质疑，你能直接把日志调出来证明每一步都做了。
3.2 第一步：过滤污染序列与未分类特征
先处理叶绿体和线粒体。这里有一个需要特别小心的点：分类学注释表里不是所有特征都有完整的分类信息，很多特征的Order或Family层级是NA。如果你直接写tax_table(ps)[, "Order"] != "Chloroplast"，NA值会被当成TRUE保留，这没问题；但如果用subset_taxa(ps, Order != "Chloroplast")，某些版本的phyloseq在对NA执行比较时会产生NA，导致这些特征被误删。
稳妥的写法是显式加上is.na判断：
# 过滤叶绿体和线粒体，保留缺失注释的特征 ps1 <- subset_taxa( ps, (is.na(Order) | Order != "Chloroplast") & (is.na(Family) | Family != "Mitochondria") ) # 顺便看一下过滤结果 ntaxa(ps1) sum(sample_sums(ps1)) / baseline_reads # 看看保留了多大比例的总reads至于那些Kingdom或Phylum级别就注释为NA的特征，要不要直接删掉？我的建议是先不要急着删。很多NA只是数据库里没有近缘参考序列，不代表它不是真实的微生物。如果删掉后特征数减少得特别多，反而说明你的数据库或分类注释这一步可能有问题。把这些NA暂时保留，留到后续流程里结合丰度和频率来综合判断。
3.3 第二步：基于频率与丰度的过滤
污染过滤之后，接着做低丰度过滤。这里我展示一套比较常用且可调整的组合：
# 频率过滤：至少在10%样本中出现 min_presence <- ceiling(0.1 * nsamples(ps1)) ps2 <- filter_taxa( ps1, function(x) sum(x > 0) >= min_presence, prune = TRUE ) # 丰度过滤：总reads数至少为10 ps3 <- prune_taxa( rowSums(otu_table(ps2)) >= 10, ps2 )filter_taxa里的匿名函数接收每个特征的计数向量x，sum(x > 0)统计的是“阳性样本数”，prune = TRUE表示直接删除不满足条件的特征。这里用了ceiling()向上取整，避免样本数乘以0.1后出现小数。
丰度过滤我用的是prune_taxa加rowSums，逻辑更直白。两条命令执行完后，再对比一下那个“基线”：
ntaxa(ps3) # 和baseline_taxa对比，看删掉了多少特征 sum(sample_sums(ps3)) / baseline_reads # 看总reads保留比例一个比较健康的过滤结果是：特征数从一两万降到三四千，但总reads保留比例还在70%~90%。如果总reads只剩不到40%，说明你的阈值设得太激进，或者原始数据质量确实很差，这时候不要盲目继续分析，先回去检查样本质量。
3.4 过滤效果评估
评估过滤效果，不能只看“删掉了多少特征”，更要看“删掉的这些reads占总reads多少”。大量低频特征虽然数量很多，但它们在总reads里的占比通常很小。如果过滤后总reads保留比例很高，说明删掉的主要是噪声；如果保留比例很低，说明你可能把真实的稀有物种也一起删掉了。
另外还要检查样本reads的分布变化：
summary(sample_sums(ps3)) min(sample_sums(ps3))这里有一个很常见的坑：过滤前每个样本reads都挺高，过滤后某个样本只剩几百甚至几十reads。原因通常是这个样本本来质量就差，大量reads分布在低频特征里。遇到这种情况，我在后面第5章会专门讨论处理方案，但你在评估阶段就要把这个信号识别出来。
4. 标准化与相对丰度转换的原理与实操
4.1 为什么不能直接用reads数相互比较
过滤做完，进入标准化环节。先讲一个最基本的道理：同一个物种在A样本里测到500条reads，在B样本里测到2000条reads，能直接说这个物种在B样本里更多吗？不能。因为可能B样本的测序深度就是A样本的四倍，所有物种的reads都被等比例放大了。
打个比方，两个班级考试，A班50人，B班200人，A班有5个人考了90分以上，B班有20个人考了90分以上。人数不同，不能直接说B班考得更好。相对丰度就是把每个样本的总数都缩放到同一个基准上，相当于把两个班的人数都换算成百分比，这样才公平。
所以标准化到相对丰度，本质上是把“绝对计数”转换成“组成比例”，让样本之间变得可比。
4.2 相对丰度转换的两种实现
在phyloseq里，相对丰度转换最简单的方式是用transform_sample_counts：
# 方法一：phyloseq原生 ps_ra <- transform_sample_counts(ps3, function(x) x / sum(x)) # 验证：每个样本的丰度总和应该等于1 colSums(otu_table(ps_ra)) %>% head()transform_sample_counts的作用是对每个样本的计数向量执行给定的函数，x / sum(x)就是把每一个reads数除以样本总reads数，得到比例。转换后otu_table里的值就是相对丰度，范围0到1。
如果你想把结果导出成表格手动检查，可以用：
# 方法二：手动计算，方便导出 otu_mat <- as.data.frame(otu_table(ps3)) otu_ra <- sweep(otu_mat, 1, rowSums(otu_mat), "/") otu_ra <- otu_ra[order(rowSums(otu_ra), decreasing = TRUE), ] write.csv(otu_ra, "relative_abundance.csv")sweep函数第一遍看可能有点绕，它的意思就是：对每一行（样本），除以该行所有值的总和。结果和方法一完全一样。
4.3 不同标准化方法的横向对比
相对丰度（TSS）是最常用的标准化，但绝不是唯一的选择。做微生物组分析时，到底该用哪种标准化，取决于你后续要跑什么分析。我在下面整理了一张对照表，建议收藏：
| 标准化方法 | 核心思想 | 适用场景 | 主要注意事项 |
|---|---|---|---|
| TSS（相对丰度） | 每个特征除以样本总reads数 | 群落组成描述、Bray-Curtis距离、PCoA | 实现简单，但对差异丰度检测可能产生假阳性 |
| CSS | 基于分位数缩放，受高丰度特征影响小 | 样本间测序深度差异大时的beta多样性 | 来自metagenomeSeq包，转换后可能出现负值占位 |
| TMM | 基于M值截尾均值，来自edgeR | 差异丰度分析前，特征间文库大小差异大时 | 通常配合edgeR的负二项模型使用 |
| CLR | 对数比变换，消除成分数据的闭合效应 | 相关性网络、CoDA分析 | 需要先处理零值，常用伪计数加1 |
| Rarefying（抽平） | 统一重抽样到相同深度 | 老牌流程常用 | 会丢弃有效数据，近年不太推荐，使用时要设随机种子 |
你可能会问，为什么我看很多文章都在用相对丰度？因为对于大多数“描述性”分析——比如展示门水平群落柱状图、算Bray-Curtis距离、做PCoA排序——TSS简单直观，结果解读成本低，审稿人也不会质疑。但如果你要做的就是看哪些物种在两组之间有显著差异，直接用TSS后的比例去做t检验，很容易出现假阳性，因为比例数据是“闭合”的，一个物种比例上升，另一个物种比例必然下降，这种结构会导致检验失效。
4.4 按分析目的选择合适的标准化方案
我习惯给建议时按分析目标来分，而不是按照标准答案背一遍。给出三类最常见的场景：
如果你做的是alpha多样性（Shannon、Chao1这些）或者beta多样性（PCoA、NMDS），直接对过滤后的原始reads数计算就行，不需要先转相对丰度。phyloseq的plot_richness和ordinate都接受原始计数，内部会正确处理。
如果你要画群落组成柱状图、热图、堆叠面积图，那就必须用相对丰度。这类图本质上展示的是“组成比例”，用原始reads画出来的图会被高深度样本完全主导，看不出比例差异。
如果你要做相关性网络分析，比如推断物种间共现关系，建议使用CLR转换，甚至用SparCC等专门为成分数据设计的方法。因为相对丰度是闭合数据，直接用TSS矩阵算Spearman相关会产生大量伪相关。这几年成分数据分析（CoDA）的论文很火，主要就是在解决这个问题。
# CLR转换示例，需要在过滤后先处理零值 library(microbiome) ps_clr <- transform(ps3, "clr") # 或者手动做：log(x + 1)后减去行均值 otu_clr <- t(apply(log(otu_table(ps3) + 1), 1, function(x) x - mean(x)))注意CLR对零值非常敏感，加1伪计数只是一个临时方案，学术上还有更严格的零值替代方法（如乘法替换），不过做探索性分析时加1足够用了。
5. 常见问题与排查技巧
5.1 过滤后样本reads数参差不齐怎么办
这是我被问到最多的问题。过滤前所有样本都在三万reads以上，过滤后有个别样本只剩800，甚至有个样本变成了0，这种情况怎么处理？
先不要慌，排查步骤分三步。第一步，看这个样本过滤前是什么水平。如果过滤前就只有5000 reads，那它本身就偏少，过滤掉低频特征后自然雪上加霜。第二步，看这个样本的reads主要被谁吃掉了。用plot(sample_sums(ps), sample_sums(ps3))把两个向量做个散点图，能直观看到哪些样本受影响最大。第三步，判断是不是污染序列占比过高。如果一个样本本来一半reads都是叶绿体，过滤后它的有效reads骤减是正常的，这个样本的真实微生物信号本来就很弱。
处理方式上，我个人的偏好是：如果有一个样本过滤后有效reads低于总reads的10%，那它在后续多元分析里大概率只是个离群点，我会直接剔除并说明原因。如果只是相对偏低，保留，因为TSS标准化本身就已经把测序深度差异缩放了。
5.2 相对丰度转换后出现全零样本或特征
这种情况多数不是代码写错，而是过滤后样本里本身就没有任何特征剩下来——通常是样本质量太差，或者过滤阈值设得太高。而如果某个特征变成全零，说明它只在被删掉的样本里出现过，或者总reads太低被丰度过滤清掉了。
遇到全零样本，我的建议是回到过滤前检查这个样本的原始reads数和物种组成，判断是建库问题还是真实存在的低生物量样本。不要试着用伪计数去“救”它，没有任何一个标准化方法能把一个空样本变成有生物学意义的样本。遇到全零特征就直接删掉即可，反正不影响任何分析，只会让矩阵更干净。
5.3 分类学注释大量NA怎么办
注释NA的问题在环境样本里非常常见，尤其是深层的门级别。有些同学一看到NA就直接过滤，觉得“用不上的信息就是垃圾”，这个做法往往会把真实稀有微生物全删光。
我的经验是分三层判断：如果NA出现在Kingdom这种最高级别，特征可能是嵌合体或未知真核污染，可以谨慎考虑删除。如果NA出现在Phylum级别但Kingdom已经被注释为Bacteria，说明它可能是一个新的、数据库里没有近缘序列的分类群，先保留，看它的丰度和频率。如果NA出现在Genus或Species这种细分类级别，完全正常，很多环境菌属本来就注释不到那么细，保留并标记为“Unclassified”即可。
5.4 与QIIME2或其他平台衔接时参数不一致
很多项目是先在QIIME2里做完DADA2，然后导出BIOM文件到R里做下游。这种情况下，QIIME2的filter-features和R里用phyloseq过滤，本质功能类似，但参数表达完全不同，容易出现“两套流程叠加过滤”的问题。
我的建议是“责任划分清楚”：如果在QIIME2里做DADA2时已经用默认参数过滤了低质量序列，那R里的过滤就只做污染序列过滤和最低限度的频率/丰度过滤，不要重复过滤两遍。记录参数尤其重要，QIIME2的filter-features里--p-min-frequency和R里rowSums >= 10不是同一个概念，混用时容易把数据清理到面目全非。
5.5 下游分析前的最终自检清单
完成过滤和标准化之后，我每次都会按这张表检查一遍，确认没有问题再跑下游分析：
| 检查项目 | 健康区间 | 为什么重要 |
|---|---|---|
| 总reads保留比例 | 过滤后/过滤前 > 70% | 太低说明过滤过度或原始数据质量差 |
| 保留特征数 | 通常在1000~10000 | 太少说明删过头，太多说明过滤不够充分 |
| 样本最小reads | 不小于总reads中位数的5%~10% | 过低的样本会成为下游离群点 |
| 是否存在全零样本 | 不应存在 | 存在则说明特定样本无有效数据 |
| 过滤前后前10优势物种 | 基本一致或略有调整 | 变化太大说明过滤把主要信号也删掉了 |
| 分类学各层级注释率 | Phylum注释率 > 80% | 过低说明数据库或聚类参数有问题 |
这张表是我自己踩过无数次坑之后总结出来的，尤其是前10优势物种这一项。有一次我调阈值调得过于激进，过滤后特征数只剩800多个，但前10优势物种几乎全变了，再一看总reads保留比例才55%，明显是把很多真实存在的优势类群因为“在部分样本里丰度低”给误删了。后来我把阈值调低，结果立刻恢复合理。
最后再分享一点个人体会：微生物组数据分析里，过滤和标准化这个环节，看起来只是“删几行、除个总数”，但实际上它决定了你后续所有分析结论的可靠性。很多人在alpha多样性和差异物种上得出矛盾结论，根子往往不在统计方法，而是前面这步没做扎实。我自己养成的习惯是每次过滤都写一个日志文件，把基线、阈值、每一步后的特征数和reads数都记录下来，这套流程最早就是为“可回溯”这件事设计的。按这套流程走完，你得到的那张相对丰度表，才是能真正支撑起后续分析的好数据。
