企业官网建设流程全解析

简介：微生物α多样性是评估局域生境群落物种丰富度与均匀度的核心指标，也是微生物组研究中最常见的生态统计量之一。这份R项目代码包面向微生物组科研人员与生物信息学初学者，提供了一套可直接运行的α多样性分析流程，涵盖数据导入、多样性指数与丰富度指数计算、Shannon指数箱线图绘制等完整步骤，并给出了Sobs、Chao、ACE等指数含义的计算说明。资源共十个文件，以两个R脚本为主要执行程序，另有txt数据与结果文件、Markdown分析报告及配置文件，便于对照项目日志和分组信息进行复现与二次修改，压缩包仅16KB，结构紧凑易用。运行该代码后，可得到标准化结果表格与分组箱线图，能直接用于论文或报告配图。目前已有261人学习，适合需要快速上手群落生态统计分析、同时希望结合真实项目代码理解指数原理的研究者。

1. 微生物 α 多样性分析是什么：用 R 把 OTU 表变成一张能放进论文的图

拿到扩增子测序的下机数据时，很多人第一反应是打开某个在线分析平台，点几下鼠标得到一张 OTU 丰度表和几个多样性指数。问题是平台不会告诉你这些指数是怎么算的，换一组参数结果就变，写方法部分时更是一头雾水。微生物 α 多样性分析正是用来回答「单个样本内部有多少物种、物种分布是否均匀」这个问题，它给出的是 Observed、Shannon、Simpson、Chao1 这样一组指数，是微生物组文章里几乎必放的图。R 语言这条路用 vegan 和 phyloseq 两个包就能把从 OTU 表到指数、稀释曲线、组间比较、出图的完整链路打通，代码可控、可复现，适合所有需要自己掌控分析过程、要应对审稿人追问的人。这篇文章按真实分析流程讲透每一步的参数和坑，新手能跟着跑通，老手能对照排查自己项目里的玄学问题。

2. 先立住数据结构：OTU 表、元数据和 phyloseq 对象怎么搭

2.1 三张核心数据：OTU 表、元数据表、系统发育树

α 多样性分析的第一步不是写代码，而是把手上的文件理清楚。绝大多数扩增子项目里你会拿到三类数据：第一类是 OTU 丰度表或 ASV 表，行是 OTU/ASV 编号，列是样本名，值是序列数或相对丰度；第二类是元数据表，记录每个样本的分组、时间点、处理浓度等信息；第三类是系统发育树，通常是 Newick 格式，只有在计算 Faith's PD 这类系统发育多样性指数时才必须用。Observed、Shannon、Simpson、Chao1、ACE 这些最常用的指数只需要前两张表，这个认知能帮你避免一开始就被树文件拖住。

我一般会先把 OTU 表转成纯数值矩阵，行名是 OTU ID，列名是样本 ID，矩阵里不能残留分类信息列。平台导出的表经常自带 taxonomy 列，读进来之后必须单独拆出去，否则后面 vegan 的函数会把它当成数值列处理，结果全是 NA。

# 读取平台导出的 OTU 表，注意跳过隐藏的注释行 otu_raw <- read.csv("otutab.csv", row.names = 1, check.names = FALSE) # 分离分类信息：假设 taxonomy 列名叫 "taxonomy" tax_col <- otu_raw$taxonomy otu_mat <- otu_raw[, setdiff(colnames(otu_raw), "taxonomy")] otu_mat[is.na(otu_mat)] <- 0 # 只保留纯数值列并转成矩阵 otu_mat <- as.matrix(otu_mat) # 读取元数据，第一列是样本 ID meta <- read.csv("metadata.csv", row.names = 1, check.names = FALSE) # 分组列转成因子，levels 顺序决定后面所有图例顺序 meta$group <- factor(meta$group, levels = c("Control", "Treat"))

逻辑说明：read.csv 里的 row.names=1 把第一列设为行名，check.names=FALSE 避免 R 自动把样本名里的横线、空格改写成点，这个参数是很多人没注意就开始踩坑的地方。tax_col 暂时保存起来，后面做物种注释相关分析时还要用。把 taxonomy 分离出来后，otu_mat 必须是一个纯数值矩阵，as.matrix 之后可以用 is.numeric() 验证一下，我曾经在一份数据里看到某个 OTU 的丰度值出现了一个不可见字符，体检时完全看不出来，算指数时才暴露。

再说回元数据表。group 列转成因子之后，levels 的顺序就是后面所有图里箱线图从左到右的顺序。很多人在这里犯了第一个错误：分组名是 High、Low、Medium 这种字母序和逻辑序不一致的名字时，R 默认按字母序排成 High、Low、Medium，和论文要求的「对照在前」冲突。养成习惯：levels 变量单独定义一次，后面所有图、检验、字母标记都复用同一个变量。

2.2 构建 phyloseq 对象：把三张表装进一个容器

phyloseq 包的价值在于把 OTU 表、元数据、物种注释、系统发育树封装成一个对象，后续所有算指数、抽平、画图都是针对这一个对象操作，不用反复匹配样本名和 OTU 名。构造代码如下：

library(phyloseq) # tax_df 由 tax_col 拆成"OTU ID + 各分类级别"的数据框，行名与 otu_mat 行名一致 tax_df <- data.frame(row.names = rownames(otu_mat), Kingdom = "Bacteria", Phylum = "NA") # 构建 phyloseq 对象，三部分可以缺省 ps <- phyloseq( otu_table(otu_mat, taxa_are_rows = TRUE), sample_data(meta), tax_table(as.matrix(tax_df)) ) # 检查对象是否完整 ps sample_names(ps)[1:5] taxa_names(ps)[1:5]

逻辑说明：otu_table() 的 taxa_are_rows=TRUE 声明行是物种、列是样本，这个参数会决定后面所有匹配方向，设反了会导致样本数和物种数混乱。sample_data() 要求传入的 data.frame 行名与 OTU 表列名完全一致，这里是我早期翻车最频繁的地方：平台导出的样本名带了批次前缀，元数据里没有，结果 sample_data 里全是 NA，后面一算指数就报错。tax_table() 接受矩阵或数据框，行名必须等于 taxa_names 的子集，构建之前先跑一次 intersect() 确认匹配情况。

参数说明：如果暂时不需要物种注释和系统发育树，可以只构建 otu_table 和 sample_data 两部分，phyloseq 允许缺省组分，后面用到时再补。tax_df 里我用了一列占位的 Kingdom 和 Phylum，实际分析中应该用数据库注释结果填写，这里只是演示结构。构建完成后 ps 这个对象可以保存成 RDS 文件，下次直接读入，省去每次重复读表的步骤。

saveRDS(ps, "ps_object.rds") # 下次分析时 ps <- readRDS("ps_object.rds")

2.3 用 vegan 还是 phyloseq 计算多样性指数

两条路我都常走：vegan 的函数更贴近生态学教科书，适合理解每个指数的计算公式；phyloseq 的 estimate_richness 是封装好的批量接口，适合快速出表。先看 vegan 的写法，理解它你才知道每个指数到底在算什么。

library(vegan) # diversity 要求样本 x 物种矩阵，所以先转置 t_otu <- t(otu_mat) shannon <- diversity(t_otu, index = "shannon") simpson <- diversity(t_otu, index = "simpson") invsimpson <- diversity(t_otu, index = "invsimpson") observed <- specnumber(t_otu) # Chao1 和 ACE 在 estimateR 里 richness <- estimateR(t_otu) chao1 <- richness["S.chao1", ] ace <- richness["S.ACE", ] # 拼成一张指数表 alpha_df <- data.frame( sample = names(shannon), shannon = shannon, simpson = simpson, observed = observed, chao1 = chao1, ace = ace )

逻辑说明：diversity() 的 index 参数支持 shannon、simpson、invsimpson 三种，其中 InvSimpson 是 1/Simpson，数值越大代表多样性越高，画图时比原始 Simpson 更直观，因为 Simpson 数值越大反而代表多样性越低，读者容易读反。specnumber() 计算的是观察到的物种数，逻辑是统计每个样本中丰度大于 0 的物种列数。estimateR() 返回一个矩阵，行名是 S.chao1、S.ACE、S.ACE.1 等，取数时用字符串行名而不是位置索引，否则换一个数据版本就错位。

参数说明：三个函数的输入都是样本 x 物种矩阵，所以必须转置。这里转置是 α 多样性分析里最容易搞错的一步，vegan 的生态学包惯例是样本在行，而测序平台的输出格式通常是物种在行，两者的差异导致代码在别人机器上跑通、在自己数据上报错的情况非常多见。

phyloseq 的写法可以把所有指数一次性算完，更适合批量探索：

alpha_ps <- estimate_richness(ps, measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "InvSimpson")) alpha_ps$sample_id <- rownames(alpha_ps) # 合并元数据，为组间比较做准备 alpha_plot <- merge(alpha_ps, sample_data(ps), by.x = "sample_id", by.y = "row.names")

参数说明：measures 参数名要写对，Observed、Chao1、Shannon、Simpson、InvSimpson 的首字母大写是 phyloseq 的约定，和 vegan 的小写写法不一致，混着用会直接报错。estimate_richness 输出里 Observed 用的就是 specnumber 的逻辑，Chao1 和 ACE 基于一阶和二阶丰度估计，对低频 OTU 非常敏感，这也是后面「抽平深度不一致时不能直接比」这个坑的根源。

2.4 四个指数怎么选：一张表说清楚

短句展开：Observed 就是数个数，深度不够时数值直接偏低，它不能用来跨样本比较，只能配合稀释曲线说明测序深度是否充分。Chao1 试图用「只出现过一次的 OTU 数量」外推总丰富度，对低频 OTU 非常敏感，因此最怕测序深度不一致。Shannon 对均匀度敏感，样本内各物种丰度越平均数值越高，适合描述群落整体格局。Simpson 反过来，优势种占比越大数值越小，如果研究场景是污染或胁迫导致某个菌暴增，Simpson 比 Shannon 更敏感。

我的习惯是正文里选两个指数讲清楚就够，Shannon 代表整体多样性格局，Chao1 代表丰富度，Observed 辅助判断测序深度，Simpson 在有明显优势种变化的场景才单独提。

3. 稀释曲线和物种累积曲线：测序深度够不够，能不能比

3.1 稀释曲线的原理和 rarecurve 的用法

α 多样性可比性的前提是各样本测序深度一致。现实中样本的序列数从几千到几十万不等，直接比较 Chao1 会被深度差异干扰。稀释曲线解决的是这个问题：把每个样本随机抽到不同的序列数，看观测物种数如何随抽取量变化。曲线进入平台期，说明当前测序深度下物种数基本被观测完整；没进入平台期，说明再测更多数据还会发现新 OTU。

library(vegan) # 按所有样本的最小总序列数作为抽样上限 min_depth <- min(colSums(otu_mat)) # step 控制抽样间隔，越小曲线越平滑但越慢 set.seed(123) rc <- rarecurve(t(otu_mat), step = 100, sample = min_depth, xlab = "Number of sequences", ylab = "Observed OTUs", label = TRUE)

逻辑说明：rarecurve 接收样本 x OTU 矩阵，所以这里传的是 t(otu_mat)。它内部对每个样本从少量序列开始逐步抽样，步长由 step 控制，每个抽样深度上计算观测物种数，最终画出曲线。sample 参数决定抽样终点，设成 min_depth 意味着所有样本都抽到同一个深度，这是组间可比的前提。

参数说明：set.seed(123) 必须放在 rarecurve 前面。稀释曲线是随机抽样过程，不设随机种子的话每跑一次曲线形状都会微调，论文里同一根曲线两次运行结果不同，这是会被审稿人抓住的复现性问题。sample 参数设成 min_depth 会丢失高深度样本的尾部信息，如果某些样本深度特别大、而最小样本深度特别小，曲线会显得参差不齐。我在实际项目中习惯把 sample 设成所有样本深度的中位数，或一个自己圈定的目标深度，比如 10000 条，这样既保证可比性又不浪费数据。

3.2 用 iNEXT 做外推：不用干等平台期

稀释曲线只能回答「当前深度下观测到了多少物种」，不能回答「再测多少序列才能到平台期」。iNEXT 的方法是用已有的丰度数据外推物种积累曲线，并给出置信区间，这比一句「曲线趋于平缓」更有说服力，也常被审稿人要求补充。

library(iNEXT) # 输入要求是列表或矩阵，行是样本、列是物种 inext_data <- as.list(as.data.frame(t(otu_mat))) # q = 0 对应物种丰富度，endpoint 是外推到的最大序列数 inext_res <- iNEXT(inext_data, q = 0, datatype = "abundance", endpoint = max(colSums(otu_mat)) * 2, se = TRUE, nboot = 200) # 提取用于画图的长表格 curve_df <- fortify(inext_res)

逻辑说明：iNEXT 支持 abundance 和 incidence 两种数据类型，扩增子 OTU 表是丰度数据，用 abundance。q 参数切换多样性阶数，q=0 是物种丰富度，q=1 对应 Shannon 熵，q=2 对应 Simpson 多样性，三种都能画曲线，但最常用的是 q=0。endpoint 设为当前最大样本深度的 2 倍，意味着你在方法部分可以写「预测测序量翻倍时的物种饱和趋势」。

参数说明：nboot 是自举次数，200 次足够出置信区间，设到 1000 会明显变慢。se=TRUE 必须开，否则没有置信区间，论文里那根曲线就少了可信度。iNEXT 对零膨胀矩阵比较敏感，如果一个样本总序列数很小，外推结果会很飘，实际操作中我会先过滤掉总序列数低于 5000 的样本再跑。

3.3 读图的三条判断标准

曲线画出来之后，判断能不能用就看三点。第一，曲线末端是否进入平台期，如果所有曲线仍在上升，那也只能如实写「sequencing depth was sufficient for comparison」这类谨慎表述，不敢写「深度充分」。第二，多条曲线在末端如果有交叉缠绕，说明样本间的排序在抽稀前后不稳定，这时候重抽一次看看差异，如果每次结果顺序都变，说明测序深度不够的样本拉低了整体的可比性。第三，看置信区间的宽度，如果外推部分的区间宽到和实测部分一样大，说明外推不可靠，这里宁可把 endpoint 调小一点。

需要特别强调一个常被误解的点：稀释曲线平了并不等于可以直接用原始深度去做组间比较。单样本曲线进入平台期只能说明这个样本的物种数随深度增长变缓，样本与样本之间的深度差异仍然存在，Chao1 依然会被高深度样本中多出来的低频 OTU 抬高。标准做法仍然是先统一抽平再做组间比较，这个操作放在下一章详细展开。

4. 组间比较与可视化：箱线图、抖动点、显著性字母一次到位

4.1 比较之前先做三件数据卫生工作

拿到指数表先别急着画图，先做三件事。第一，剔除总序列数过低的样本，比如小于 1000 条的，这种样本的 Observed 和 Chao1 基本没有参考价值，留着只会引入噪声。第二，检查分组样本量，每组至少 3 个样本才能给出均值和误差线，少于 5 个时非参数检验的效能很低，p 值在 0.05 附近来回跳，结论经不起推敲。第三，把分组列明确设为因子并指定 levels，这一步定死后面所有图例顺序和检验的分组顺序。

# 按照总序列数过滤低深度样本 depth <- colSums(otu_mat) keep_samples <- names(depth[depth >= 1000]) otu_mat_filtered <- otu_mat[, keep_samples] meta <- meta[keep_samples, , drop = FALSE]

4.2 箱线图加抖动点的标准出图代码

ggplot2 画箱线图叠加抖动点是最稳的组合。箱线图给出分布轮廓，叠加的抖动点让每个样本的位置可见，审稿人一眼能看出每组的样本量和离群样本。

library(ggplot2) # alpha_plot 来自上一章合并后的数据，包含 group、shannon 列 p <- ggplot(alpha_plot, aes(x = group, y = shannon, fill = group)) + geom_boxplot(outlier.shape = NA, width = 0.5) + geom_jitter(width = 0.15, size = 2, alpha = 0.6) + scale_fill_manual(values = c("Control" = "#4DBBD5", "Treat" = "#E64B35")) + labs(x = "Group", y = "Shannon index") + theme_bw(base_size = 14) + theme(legend.position = "none") ggsave("shannon_boxplot.pdf", p, width = 4, height = 3.5, dpi = 300)

逻辑说明：outlier.shape=NA 让箱线图不再重复画出离群点，因为抖动点已经把离群样本表现出来了。width=0.15 的 jitter 让各样本点围绕组中心小幅散开，避免大量样本重叠在一起。alpha=0.6 是半透明效果，样本量大的时候减少遮挡，能看清点与点之间的重叠程度。

参数说明：如果每组样本量超过 30，把 size 降到 1.5、alpha 降到 0.4，否则点会糊成一团；如果每组样本量很小，比如 5 个以内，可以考虑去掉 outlier.shape 参数，保留箱线图的离群点展示。颜色值用的是蓝红搭配，具体按期刊偏好调整，我一般会把颜色值统一成一个向量放在脚本顶部，方便全局修改。

4.3 检验方法：为什么默认不用 t 检验

α 多样性指数的分布通常不满足正态性，Chao1 和 Observed 尤其明显，受少数稀有 OTU 影响会出现长尾。所以两组比较我用 Wilcoxon 秩和检验，多组比较用 Kruskal-Wallis 检验，整体显著之后再做 Dunn 事后检验。t 检验不是完全不能用，但要先通过 Shapiro 正态性检验，多数情况下过不了，直接换秩检验更省事，方法部分也更好写。

# 两组比较 wilcox.test(shannon ~ group, data = alpha_plot) # 多组比较 kruskal.test(shannon ~ group, data = alpha_plot) # 事后 Dunn 检验，需要 FSA 包 library(FSA) dunn_test <- dunnTest(shannon ~ group, data = alpha_plot, method = "bh") dunn_test$res

逻辑说明：wilcox.test 的公式接口要求 group 是两水平因子，多水平会报错。kruskal.test 不限制水平数，但只返回整体 p 值，告诉你至少有一组不同，却说不清哪些组之间有差异。dunnTest 的 method="bh" 使用 Benjamini-Hochberg 校正，比 Bonferroni 温和，多组比较场景更常用。

参数说明：如果比较的是 Simpson 指数，记得确认方向。Simpson 数值越小多样性越低，但画图时我会转成 InvSimpson，因为 InvSimpson 是「有效物种数」的语义，数值越大多样性越高，读图不会反。Dunn 检验输出的 p.adj 是两两比较的校正后 p 值，字母标记需要自己从这张表生成。

4.4 显著性字母标记的一种可靠做法

dunn_test 的输出里只有两两比较的 p 值，论文图要的是 a、b、c 这样的字母分组。常见做法是用 multcompView 包把 p 值矩阵转成字母，规则是：差异不显著的组共享同一个字母，显著的组之间标记不同字母。初学者经常在这里手动标字母导致错位，尤其是组数多的时候。

library(multcompView) # dunn_test$res 里包含 Comparison 和 p.adj 两列 pvals <- dunn_test$res$p.adj names(pvals) <- dunn_test$res$Comparison # 生成字母标记 letters_result <- multcompLetters(pvals, threshold = 0.05) print(letters_result$Letters)

代码说明：multcompLetters 的输入是带名字的 p 值向量，名字格式必须形如 "A-B" 的连接符，dunnTest 输出的 Comparison 列恰好是这个格式，可以直接用。threshold 是显著性阈值，默认 0.05，如果组数超过 5 个，字母组合会变得很复杂，可以考虑用 0.01 使结果更保守。生成字母之后，把字母放在箱线图顶部需要一个 annotate 层，不同组的字母位置由该组的最大值坐标决定，这一步可以写成循环自动生成。

5. 微生物 α 多样性分析避坑：七条血泪踩坑记录

5.1 抽平深度不一致，组间差异全是假的

现象：某模拟项目里处理组和对照组的 Chao1 差异 p=0.003，结论看起来非常漂亮。检查数据后发现处理组的平均序列数是对照组的两倍，再把两组统一抽平到同一深度重新计算，p 值变成 0.31。

原因：Chao1 关注的是低频 OTU，而低频 OTU 恰恰最容易被测序深度影响。序列数多的样本天然有更多稀有 OTU，Chao1 偏高是深度效应，不是真实生物学差异。这个问题在 Shannon 上影响小一些，但在 Chao1 和 Observed 上非常明显。

解决：统一样本深度之后再做任何组间比较。最保守的做法是按所有样本的最小序列数抽平，缺点是丢失高深度样本的信息；另一种做法是抽平到目标深度并手动剔除深度不足的样本。

library(vegan) set.seed(2024) # 按最小序列数抽平，之后矩阵变成样本 x OTU rarefied <- rrarefy(t(otu_mat), sample = min_depth) # 抽平后重新计算多样性 shannon_rare <- diversity(rarefied, index = "shannon") chao1_rare <- estimateR(rarefied)["S.chao1", ]

参数说明：rrarefy 是随机抽稀，结果受随机种子影响，所以 set.seed 必须写。如果担心一次抽平的偶然性，可以抽 100 次取均值，代价是计算时间。我在正式项目里通常抽 100 次取均值，探索阶段用一次抽平就够了。

5.2 diversity 返回一堆 NaN

现象：diversity() 跑完 shannon 列全是 NaN，simpson 也有 NA，但程序没有报错，只是输出里出现一堆警告。

原因：OTU 表里有全 0 的样本列，或者某个样本在过滤后变成了空样本。Shannon 的公式里出现 log(0)，R 返回 -Inf，接着参与其他运算就变成 NaN。

解决：先检查每个样本的总序列数，删掉零样本列再计算。

col_sums <- colSums(otu_mat) print(col_sums[col_sums == 0]) otu_mat <- otu_mat[, col_sums > 0]

实际操作中，全 0 列通常来自平台导出文件里的多余空白列，或者是抽平后某些低深度样本被抽成 0。删除之前先打印出来确认一下是哪些样本，避免误删有效数据。

5.3 图例顺序和显著性字母对不上

现象：图里箱线图的顺序是 Control、High、Low、Medium，而显著性字母是按处理逻辑顺序标的，图的左边和字母对不上，审稿人一眼看出标注错位。

原因：读取元数据时 group 列是字符型，R 自动按字母序排列，而手动标记字母时又按 Excel 里看到的顺序写，两套顺序不一致。

解决：从源头统一因子 levels，并在标记前打印 unique 确认实际顺序。

alpha_plot$group <- factor(alpha_plot$group, levels = c("Control", "Low", "Medium", "High")) unique(alpha_plot$group)

这个坑是最容易犯也最容易被审稿人抓的。图是好的，数据是好的，只有顺序是乱的。我的习惯是把 levels 写成一个变量放在脚本顶部，后面所有图、所有检验、所有字母标记都从同一个变量取，绝不手写第二次顺序。

5.4 rarecurve 报错或者画出一堆乱线

现象：rarecurve 报错提示矩阵里有非数值列，或者画出来的曲线数量远多于样本数量，线乱得没法看。

原因：读入的 OTU 表里含着 taxonomy 列，as.matrix 之后整个矩阵变成字符矩阵；或者列名有重复，R 自动加后缀导致样本名对不上。

解决：读入时把分类列拆出去，计算前加断言检查矩阵类型和列名唯一性。

stopifnot(is.numeric(otu_mat)) stopifnot(!any(duplicated(colnames(otu_mat))))

我把这两行断言放在所有计算之前。很多平台导出的 Excel 表会把数值列里混进一个空格或不可见字符，is.numeric 这一行能拦住绝大多数脏数据问题。colnames 重复的问题更隐蔽，常见于合并多批次数据时样本 ID 覆盖，提前检查总比事后排查快。

5.5 iNEXT 在大数据量样本上卡死

现象：iNEXT() 跑一个几百个样本的数据集，半小时不出结果，看起来像死机。

原因：iNEXT 对每个样本都做自举重抽样，nboot=1000 加上样本数量大时计算量呈倍数增长。

解决：先在小样本子集上跑通确认参数，再全量运行；或者把 nboot 降到 100 先看曲线形状。

# 先跑前 6 个样本 small_res <- iNEXT(inext_data[1:6], q = 0, datatype = "abundance", endpoint = 50000, se = TRUE, nboot = 100)

这个坑的本质是调试时用了全量数据。花时间等 iNEXT 跑完是很多项目里浪费时间的重灾区。另外 endpoint 不必设到最大深度的两倍，设到 1.2 倍常常就能看出平台期趋势，曲线后半段的意义远小于前半段。

5.6 phyloseq 的树和 OTU 表对不上，PD 指数全是 NA

现象：用 picante 包或 phyloseq 自带的函数计算 Faith's PD，结果全是 NA，或者干脆报错提示树与 OTU 表不匹配。

原因：OTU 表里的物种名在系统发育树里不存在，或者树里有多余的 tips，匹配函数返回空值。

解决：先取交集，同步修剪 OTU 表和树。

library(phyloseq) common_taxa <- intersect(taxa_names(ps), phy_tree(ps)$tip.label) ps_pruned <- prune_taxa(common_taxa, ps)

细节说明：prune_taxa 只按 taxa_names 过滤 OTU 表，树本身不会自动剪掉多余的 tips。只要 OTU 表里的物种都在树里，多余的 tips 不影响 PD 计算。真正的问题在于 OTU 表里有树外名字时，匹配会静默跳过它们，部分数据被丢弃而你不自知。所以计算前先跑 intersect，再对比长度是否一致。

5.7 分组样本量太小，p 值像随机数

现象：每组只有 3 个样本时，同一个数据集换一次抽平随机种子，p 值在 0.01 和 0.2 之间反复横跳。

原因：非参数检验在样本量极小时效能非常低，秩和检验统计量的离散程度太大，p 值对个别样本的微小变化极其敏感。

解决：要么补充样本，要么在结果里明确标注样本量限制。为了撑起稳健的结论，每组至少要有 5 个样本，理想的微生物组项目每组装载 10 个以上再谈组间比较。样本量不够时，我不会在论文里声称某指数组间显著，只会写成趋势性差异并如实报告 n 值。

6. 进阶技巧：把整套 α 多样性流程封装成一个可复用的函数

真实项目从来不是跑一组数据就结束，而是十几个样本、多个分组、几轮比较，每换一次分组或换一个指数都要重新复制几十行代码。我的做法是把整套流程封装成函数，输入一个 phyloseq 对象和分组列名，输出指数表、检验结果和图对象，一次调用出全套。

run_alpha_analysis <- function(ps, group_col, measures = c("Observed", "Chao1", "Shannon", "Simpson", "InvSimpson")) { alpha <- estimate_richness(ps, measures = measures) alpha$sample_id <- rownames(alpha) meta_df <- data.frame(sample_data(ps)) meta_df$sample_id <- rownames(meta_df) merged <- merge(alpha, meta_df[, c("sample_id", group_col)], by = "sample_id") merged[[group_col]] <- factor(merged[[group_col]]) formula_str <- paste("Shannon ~", group_col) kw <- kruskal.test(as.formula(formula_str), data = merged) p <- ggplot(merged, aes(x = !!sym(group_col), y = Shannon, fill = !!sym(group_col))) + geom_boxplot() + geom_jitter(width = 0.15) + theme_bw() list(alpha_table = merged, test = kw, plot = p) } result <- run_alpha_analysis(ps, "group")

这个函数已经包含前面提到的因子顺序处理、检验公式构造和出图。关键点是 group_col 用字符串传入时，在 ggplot 内部要用 !!sym() 才能让列名被正确识别，这是 R 语言 tidy evaluation 的经典问题，很多人第一次写这类函数就在这里被绕晕。

进阶的下一步是让循环跑遍所有指数，每个指数做一次检验，最后输出一张整合表格。我会在函数里加一个索引参数，让 Shannon 和 Chao1 分别出结果，再用 dplyr 把多张表 bind_rows 汇总。项目里数据更新了、分组调整了、或者审稿人要求换显著性水平，都只需要改一行参数重新跑，不用整个分析从头再来。

我在这件事上吃过一次亏：某次返修意见要求把所有图的字体统一成 Arial，当时所有 ggplot 主题写在主脚本里，我只好改了 20 多张图的 theme 代码。后来我把主题也封装成函数，所有图共用一套主题，再遇到这类全局修改只需要改一处。代码的封装度直接决定你面对修改时的返工量，宁可多花一小时写函数，不要手动画 20 次图。这套流程我已经在多个模拟项目里反复用，每次改数据重跑都只需一两分钟，最省时间的做法是提前把函数写好，希望帮到你。

本文还有配套的精品资源，点击获取