企业官网建设流程全解析

简介：面向微生物生态研究者和R语言初学者的α多样性分析项目代码包，聚焦局域均匀生境下单个样本的物种丰度与多样性评估，涵盖Shannon、Simpson等多样性指数与Sobs、Chao、ACE等丰富度指数的原理及R实现，可应用于土壤、肠道、环境等不同来源的微生物群落数据。压缩包共十个文件，大小仅16KB，以R脚本为核心，配套四个txt文件保存分组信息与分析输出，两个md文档说明分析思路与报告框架，另含gitignore等工程配置；脚本覆盖数据导入、指数计算和箱线图可视化，便于替换数据后快速运行。已有261人学习下载，项目结构清晰，未采用复杂封装，代码直观易改，不仅能生成不同组别Shannon指数分布图以比较组间差异，还可结合md文档理解各指数的生态学含义，适合作为微生物多样性分析的教学示例或论文研究参考。

1. 拿到OTU表先别急着做图：微生物α多样性分析在R里到底要解决什么

拿到16S扩增子测序的OTU表，很多人第一反应是“算几个α多样性指数，画个箱线图就完事”。真正做起来才发现，Shannon、Chao1背后的抽平深度、物种覆盖度、统计检验取舍，每一项都能决定结论站不站得住。微生物α多样性分析在R里做，本质是把OTU计数表变成一组可解释、可复现、敢写进文章的数字和图形，核心工具是vegan、tidyverse和ggplot2三件套。这篇笔记写给两类人：刚拿到测序结果、想自己完成分析的生信新手；已经跑通流程、但被抽平和统计检验反复折腾的研究者。下面按“原理—脚本—可视化—避坑—进阶”递进，每段代码都能直接复制到你的项目里改着跑。

2. 把α多样性的三个底层问题说透：指数语义、样本量效应和R的选型理由

2.1 Shannon、Chao1、Simpson各自在回答什么问题

先从指数语义讲起。Shannon熵来自信息论，公式是 -Σ p_i·ln p_i，p_i为第i个OTU的相对多度。它对稀有种敏感，低丰度物种的数量变化会直接改变熵值，所以它综合反映“丰富度+均匀度”两个维度。Simpson指数常用形式是 1-Σp_i^2，它被优势种主导，群落里如果有一个物种占到90%以上，Simpson会迅速掉下来，所以它回答的是“优势度集中程度”。Chao1用公式 S_est = S_obs + F1^2/(2F2) 估计尚未观测到的物种数，F1是singleton数量、F2是doubleton数量，它回答的是“继续加大测序还可能发现多少物种”。

选择上我的习惯是常规文章同时报告Shannon和Chao1，一个偏综合多样性，一个偏物种储备，审稿人看到这两个数基本就放心。低多样性体系（比如宿主组织样本）里，Simpson比Shannon更容易拉开组间差异，因为优势种占比的变化往往被Shannon里的中丰度物种稀释。还要注意vegan里的命名：index="simpson"返回的是Gini-Simpson形式（1-Σp_i^2），数值越大多样性越高；index="invsimpson"返回的是1/(Σp_i^2)，趋势一致但尺度完全不同。投稿时务必在方法部分写清楚你用哪一种。

2.2 抽平深度为什么决定了你的分析结论是否可信

假设A样本测了5万条reads、B样本只测了1万条reads，即使两者来自完全相同的群落结构，A观测到的物种数也天然比B多。这是测序深度带来的机械偏差，不是生物学差异。α多样性比较的前提是所有样本在同一测序深度下进行，行业里的标准做法是随机抽平（rarefaction）：对每个样本不放回地随机抽取N条reads，N通常取所有样本的最小深度。这样消除深度差异，代价是测序深的样本损失一部分信息。

抽平深度到底取多少，看起来像玄学，其实有迹可循。最小样本只有2000条reads、多数样本有3万条时，硬抽平会让OTU表变得极其稀疏，Chao1估计方差巨大，Shannon被严重低估。我一般先用summary(depth)和直方图看分布，如果存在极端低深度样本，先剔除再抽平，而不是被最小深度绑架。土壤这类高多样性样本建议抽平深度不低于10000，宿主组织这类低多样性样本可以降到1000左右，具体阈值按实际建库质量定。

这里要特意区分一个常见误用：把每列除以总数得到相对丰度表，以为这样也能消除深度差异。归一化没有消除“观测物种数”的差异，反而破坏了Chao1的计数语义——singleton和doubleton的判定必须基于绝对整数计数，归一化之后F1、F2的数值彻底失真。rrarefy抽平是保留整数计数结构的常规路径，别用相对丰度替代它。

2.3 为什么这个场景选R而不是Python或QIIME2

生态学领域的多样性指数计算，R的vegan包是事实标准。diversity()、specnumber()、estimateR()、rrarefy()这些函数经过大量论文验证，计算语义明确，边界情况处理得比自写脚本稳。QIIME2能输出α多样性结果，但组间统计检验和自定义作图最终还是要回到R；Python的scikit-bio也不是不能算，只是生态学社区的数据格式、中间产物和审稿人习惯都更偏向R生态。

R的另一个优势在数据处理链路本身。biom或TSV格式的OTU表读入后，需要转置、过滤、合并分组信息，tidyverse的read_delim、left_join、group_by这套管道操作非常顺手，ggplot2出图质量也足够发表。如果你的上下游已经跑在QIIME2里，最稳妥的组合是上游出OTU表、R做α多样性指数和统计验证，这个分工是很多项目里沟通成本最低的组合。有一个细节必须记住：vegan::rrarefy是随机抽平，不固定随机种子的话，同一份数据每次运行得到的OTU表和指数都会有轻微波动。后面脚本里直接在抽平前set.seed，这是可复现分析的第一条底线。

3. 从OTU表到多样性指数：一份可直接跑的R脚本与参数说明

3.1 读入并清洗OTU表：行与列的方向别搞反

常见输入格式是上游流程导出的TSV：第一列是OTU ID，后面每列是一个样本，值为reads计数。vegan的diversity()和rowSums()默认把“行”当作样本维度，所以第一步就是把它转置成“行=样本、列=OTU”。这一步方向搞反的话，代码不会报错，但算出来的所有指数全部是错的。

# otu_table.txt: 第一列为OTU ID，其余列为样本计数 library(tidyverse) otu_raw <- read_delim("otu_table.txt", delim = "\t") # 把第一列重命名为OTU_ID，避免各流程输出列名不一致 names(otu_raw)[1] <- "OTU_ID" otu_raw <- column_to_rownames(otu_raw, "OTU_ID") # 转置为 行=样本、列=OTU otu_t <- as.data.frame(t(otu_raw)) # 过滤总计数<=2的OTU，降低spurious OTU干扰 otu_filt <- otu_t[, colSums(otu_t) > 2] # 转为整数矩阵，vegan的计数接口要求整数 otu_mat <- as.matrix(otu_filt) storage.mode(otu_mat) <- "integer" # 防御性检查：确认行名是样本ID head(rownames(otu_mat)) head(colnames(otu_mat))

这段代码里read_delim按制表符读取，column_to_rownames把第一列升为行名，t()完成转置。colSums(otu_t)统计每个OTU在所有样本中的总reads，大于2的保留，这一步是低丰度过滤。storage.mode改成integer是为了让vegan内部计数逻辑不出现浮点边界问题，如果输入数据其实是0到1的比例数，这行代码会把大量值截断成0，反而能帮你暴露数据问题。

参数上的选择：delim="\t"对应制表符分隔文件，如果上游导出的是CSV就改成delim=","。过滤阈值>2是多数16S项目的默认值，测序深度超过5万reads时可以提高到>5或>10。head(rownames(otu_mat))这一步是人工确认方向，看到的是样本ID才算正常，如果看到的还是OTU编号，说明转置没生效。

3.2 抽平深度怎么定：先看分布再动手

# 统计每个样本的测序深度 depth <- rowSums(otu_mat) summary(depth) hist(depth, breaks = 30, main = "Sequencing depth per sample") # 剔除低于5000 reads的样本，阈值按项目实际情况调整 keep <- depth >= 5000 otu_sub <- otu_mat[keep, ] # 用剩余样本的最小深度作为抽平深度 rare_depth <- min(rowSums(otu_sub))

summary给出深度的六个分位数，hist画出分布形状，这一步能直接看出有没有极端低深度样本。如果min(depth)和第三四分位数差一个数量级，优先剔除低深度样本而不是硬抽平。我这里设的阈值5000是经验值：土壤、肠道内容物这类高多样性样本建议不低于10000；宿主组织活检这类低生物量样本可以降到1000。把rare_depth单独保存，下一步抽平直接引用，避免在两个地方改同一个数值。

3.3 用vegan计算Shannon、Simpson、Chao1与观测物种数

library(vegan) # 随机抽平，必须固定随机种子 set.seed(2024) otu_rare <- rrarefy(otu_sub, sample = rare_depth) # 基于抽平后的表计算指数 shannon <- diversity(otu_rare, index = "shannon") simpson <- diversity(otu_rare, index = "simpson") observed <- specnumber(otu_rare) chao1_est <- estimateR(otu_rare)["S.chao1", ] # 组装成数据框 alpha_df <- tibble( sample = names(shannon), shannon = shannon, simpson = simpson, observed = observed, chao1 = chao1_est )

rrarefy是vegan的随机抽平函数，sample参数传抽平深度。set.seed(2024)保证同一份数据在任何机器上抽平结果一致，这个数字本身不重要，但必须固定。diversity计算Shannon和Simpson，specnumber统计每个样本的非零OTU数，estimateR返回多行估计结果，["S.chao1", ]取出Chao1那一行。

提示：抽平后可以用rowSums(otu_rare)验证每一行是否完全相等。如果出现不一致，说明输入矩阵混入了负值或非整数，回去检查storage.mode那一步。

参数上，index="shannon"和"simpson"是vegan接受的字符串，拼写错误会直接报错。estimateR在样本里没有doubleton时可能输出Inf，这种情况很常见，具体处理放第5章讲。

3.4 合并分组信息并导出整洁表格

# group_info.txt: 两列，sample和group group_df <- read_delim("group_info.txt", delim = "\t") alpha_df <- alpha_df %>% left_join(group_df, by = "sample") # 固定分组顺序，避免ggplot2按字母序乱排 alpha_df$group <- factor(alpha_df$group, levels = c("CK", "Treatment")) # 导出CSV，后续统计和作图都基于这个文件 write_csv(alpha_df, "alpha_diversity.csv")

left_join按sample列匹配分组信息，factor通过levels显式设定分组顺序，这样箱线图x轴顺序和后续统计检验的组顺序都受控。分组文件里的样本ID必须与OTU表完全一致，短横线和下划线混用是匹配失败最常见的原因。导出的CSV既是后续统计的输入，也可以直接作为投稿补充材料。

4. 组间比较与可视化：箱线图、小提琴图和统计检验的连招

4.1 用ggplot2画分组箱线图，样本点必须露出来

library(ggplot2) p <- ggplot(alpha_df, aes(x = group, y = shannon, fill = group)) + geom_boxplot(outlier.shape = NA, width = 0.5) + geom_jitter(width = 0.15, size = 1.8, alpha = 0.6) + scale_fill_manual(values = c("#4DBBD5", "#E64B35")) + theme_classic(base_size = 14) + labs(x = NULL, y = "Shannon Index") ggsave("shannon_group.pdf", p, width = 5, height = 4, useDingbats = FALSE)

outlier.shape=NA把箱线图的离群点隐藏，由geom_jitter展示所有原始数据点，避免同一个数值被画两次。width=0.15控制抖动横向范围，size=1.8控制点大小，alpha=0.6让重叠点更透明。theme_classic是期刊最常见的白底细边框主题，比默认灰色网格更适合直接进文章。

如果组内样本量充足（每组≥8个），可以考虑用小提琴图展示分布密度：

p_violin <- ggplot(alpha_df, aes(x = group, y = shannon, fill = group)) + geom_violin(trim = FALSE, alpha = 0.4) + geom_jitter(width = 0.1, size = 1.5) + theme_classic(base_size = 14)

小提琴图能展示分布形态，但组内样本少于8个时形状极不稳定，容易给审稿人留下过度解读的印象。样本少时箱线图加抖动点更稳妥。PDF导出记得加useDingbats=FALSE，否则部分期刊排版系统会把Dingbats字体显示成乱码。

4.2 先做正态性检验再选统计方法：非参数检验是默认值

Shannon这类指数在每组3到5个生物学重复时经常不是正态分布，直接套t检验的前提不成立。先做Shapiro-Wilk检验再决定统计方法。

# 按组做正态性检验 alpha_df %>% group_by(group) %>% summarise(p_shapiro = shapiro.test(shannon)$p.value) # 两组比较：默认非配对的Wilcoxon秩和检验 wilcox.test(shannon ~ group, data = alpha_df) # 三组及以上：先做Kruskal-Wallis总检验 kruskal.test(shannon ~ group, data = alpha_df) # 两两比较并做FDR校正 library(rstatix) pairwise <- alpha_df %>% wilcox_test(shannon ~ group, p.adjust.method = "fdr")

shapiro.test的p值小于0.05说明该组数据显著偏离正态，这时Wilcoxon是更稳妥的默认选择。wilcox.test的公式接口shannon ~ group会自动按group拆两组。kruskal.test是多样本非参数总检验，先看总体是否存在组间差异，再决定是否值得做两两比较。rstatix包的wilcox_test输出包含p.adj列，比基础R的pairwise.wilcox.test更适合放进管道流。

参数上的选择：p.adjust.method="fdr"对应Benjamini-Hochberg校正，适合2到4组的两两比较，比bonferroni保留更多统计功效。每组样本量小于5时，Wilcoxon的精确p值会呈现离散化，比如0.1或0.066，这时不要只依赖p值，建议同时报告效应量，用rstatix::wilcox_effsize计算，或者直接看箱线图的中位数差异。p值只回答“有没有差异”，不回答“差异多大”，α多样性分析里组间分布重叠程度同样重要。

5. α多样性分析避坑指南：五个常见翻车现场与排查路径

这一节记录五个最容易让分析结论翻车的点，每一条都是“现象—原因—解决”的完整链路。前两个坑会直接毁掉结果，后三个坑则是代码不报错但结论失真，属于更难排查的类型。

5.1 抽平深度被极端样本拖垮，稀有种被成批丢弃

现象：某个样本只有几百条reads，其余样本有两三万，抽平后所有样本Shannon暴跌，原本显著的组间差异消失。

原因：抽平深度取的是所有样本的最小值，一个离群低深度样本把全体样本的深度上限拉低。大量稀有OTU在抽平中被随机丢弃，多样性信息损失过大。

解决：先看summary(depth)和直方图，剔除严重低深度的样本后再抽平。如果低深度样本恰好是某个处理组的全部样本，说明该组测序本身失败，需要补测而不是强行抽平比较。

5.2 Chao1算出Inf，低频OTU太少时估计发散

现象：estimateR输出的S.chao1列为Inf，或附带“estimation is very poor”警告。

原因：Chao1公式分母是2*F2。当样本里不存在doubleton（出现两次的OTU）时，分母为0，估计值发散到无穷。这在低深度或低多样性样本里很常见。

解决：这种情况直接在结果里改用Observed species，并在方法部分写清楚“该样本因低频OTU过少，Chao1估计不稳定，采用观测物种数代替”。不要把一个Inf值硬保留进图表，审稿人看到Inf会直接质疑整个流程。

5.3 行列方向搞反，代码不报错但结果全是错的

现象：算出来的Shannon全部在0到0.1之间，或者specnumber结果全是1，数据看起来异常整齐。

原因：OTU表没有转置就直接传给vegan。vegan默认对每一行计算指数，此时每行是某个OTU、每列是样本，计算语义完全错乱，而代码本身不报错。

解决：在清洗阶段统一转置，并用head(rownames(otu_mat))确认行名是样本ID。这种错误是典型的黑匣子陷阱，数字看起来正常但含义全错，防御性检查必须做。

5.4 相对丰度表被当成计数表，Chao1完全失真

现象：用公司返回的门水平相对丰度表计算，Shannon值能精确到小数点后四位，Chao1和Observed几乎相等。

原因：相对丰度表已经归一化到100%，丢失了测序深度信息。Chao1依赖singleton和doubleton的绝对计数，比例数据下F1、F2统计完全失真。

解决：α多样性指数只用原始OTU计数表计算。如果分析起点只有相对丰度表，最多计算Observed species和Simpson这类基于比例的指数，不能计算Chao1，也不应该做抽平。遇到这种数据，回上游找原始counts表才是正路。

5.5 每组只有两三个重复，Wilcoxon检验失去意义

现象：wilcox.test不报错，但p值要么是离散值0.1，要么提示“cannot compute exact p-value with ties”。

原因：非参数检验在小样本下统计功效极低，组内出现重复值就做不了精确检验，p值分辨率受限。

解决：报告效应量而不是只依赖p值，用rstatix::wilcox_effsize计算秩效应大小。如果所有组都没有4个以上生物学重复，α多样性组间差异在统计上本身就很难下结论，这是实验设计问题，不是代码能弥补的。开始实验前就按每组至少5个重复设计，能省掉后面所有检验相关的麻烦。

6. 进阶习惯：把α多样性分析脚本收敛成带参数的函数

等所有指数和图表都稳定之后，我会把散落的脚本收敛成一个函数，输入只留路径和两个参数。这样换数据集时不用改代码，也避免“上周跑的结果这周复现不出来”的尴尬。

compute_alpha <- function(otu_path, group_path, output_dir = "results", depth_min = 5000, seed = 2024) { dir.create(output_dir, showWarnings = FALSE) otu <- read_delim(otu_path, delim = "\t") names(otu)[1] <- "OTU_ID" otu <- column_to_rownames(otu, "OTU_ID") otu <- as.data.frame(t(otu)) otu <- as.matrix(otu) storage.mode(otu) <- "integer" depth <- rowSums(otu) cat("Removed", sum(depth < depth_min), "samples below depth", depth_min, "\n") otu <- otu[depth >= depth_min, ] set.seed(seed) rare_depth <- min(rowSums(otu)) cat("Rarefaction depth:", rare_depth, "\n") otu_rare <- rrarefy(otu, sample = rare_depth) alpha <- data.frame( sample = rownames(otu_rare), shannon = diversity(otu_rare, "shannon"), simpson = diversity(otu_rare, "simpson"), chao1 = estimateR(otu_rare)["S.chao1", ], observed = specnumber(otu_rare) ) groups <- read_delim(group_path, delim = "\t") alpha <- left_join(alpha, groups, by = "sample") write_csv(alpha, file.path(output_dir, "alpha_diversity.csv")) saveRDS(sessionInfo(), file.path(output_dir, "session_info.rds")) alpha } # 调用 result <- compute_alpha("otu_table.txt", "group_info.txt", output_dir = "results_2025")

这里有几个我之前踩过的坑。函数内部一定要把剔除样本数和抽平深度打印出来，否则换数据集时根本不知道有多少样本被静默丢掉。saveRDS(sessionInfo())是最近养成的习惯，把R版本和包版本随结果一起存档，半年后复现时至少能还原当时的运行环境。每个项目单独用日期命名输出目录，避免覆盖旧结果。

这个函数不算优雅，但它替我省下了大量重复劳动。现在拿到任何一份OTU表，我的工作流不是重新复制粘贴脚本，而是改两个路径、跑一遍、看日志输出，确认没有样本被异常剔除，再去看图形和统计结果。代码组织到这一步，α多样性分析才算真正跑通了。希望帮到你。

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