华大员工队列 i99 肠道菌群 α 多样性分布特征分析

摘要

背景与目的：华大员工队列 i99 拥有 10,201 名参与者、21,997 次鸟枪法宏基因组测序（MetaPhlAn 4 物种级相对丰度），但尚未对其肠道菌群 α 多样性的分布特征、时间变异与可重复性做过系统刻画。本分析基于 20% 个体级随机抽样（2,040 人 / 4,390 次 run，4,383 次有效），回答四个问题：分布形态、采样时期差异、个体内 vs 个体间变异来源、技术因素（测序深度）的影响。

主要结果：① α 多样性呈显著左偏非正态分布，Shannon 中位 3.030（IQR 2.55–3.49）、richness 中位 161 SGB（IQR 113–210），落在已发表成人粪便宏基因组研究的常见区间。② 近窗口（2022–2025）α 多样性低于早窗口（2016–2018）：Shannon Δ=−0.148（混合模型，p=7.4e-11），深度匹配子集中 richness Δ=−43.13（p=2.3e-42），配对子队列方向一致；但效应量普遍很小（21 对逐年比较中仅 1 对达大效应），且与测序批次完全共线，只能表述为「批次/时期相关差异」。③ 立项时「个体间差异 >> 个体内波动」的假说被推翻：个体身份仅解释 Shannon 方差的 11.0%，年份仅 3.7%，85.3% 为个体内残余波动；个体内/个体间 SD 比为 1.12，个体内 CV 均值 20.0%。④ 测序深度非主要驱动（深度解释 richness 的 R² 全样本 0.0216、剔除低尾后 0.0010），10M reads 以上 richness 已饱和；但低深度样本集中于 2018 年，使全样本时期效应被低估约 2.4 倍。

方法学要点与自我纠错：全程以个体为独立单位（聚类 bootstrap / 混合模型），避免 run 非独立导致的 p 值偏小；执行中经 Genpilot 对话推翻了原定的 Kruskal-Wallis 方案，并把解析式稀疏化从主方法降级为敏感性分析；交付前经批判性评审，主动推翻了自己「低 ICC 归因于采样窗口」的结论，改用同个体内配对分析（窗口内变化 0.651 反而大于跨窗口 0.589，Wilcoxon p=0.020），相应表述降级为待验证假设。

结论与限制：i99 队列 α 多样性的变异以个体内波动为绝对主导，个体间差异只占约十分之一，因此单次采样难以刻画一个人的「菌群多样性水平」。最大限制是缺少年龄/性别/BMI/用药/饮食等表型与提取批次/文库/机台等技术元数据，使“时期差异”无法与平台换代解耦，也使 85% 的个体内残余方差无法拆分为生物学波动与技术噪声。

引言与背景

背景

肠道菌群的 α 多样性（个体内群落的物种丰富度与均匀度）是微生物组研究中最常用的汇总指标，与宿主代谢、免疫及多项疾病状态相关。在大型人群队列中解读 α 多样性，有三个基础问题必须先回答清楚：它的分布形态是什么（决定后续统计方法）、它在同一个人体内有多稳定（决定单次采样能否代表个体）、以及哪些技术因素会污染它的表观差异（决定关联分析是否需要校正）。

研究对象

i99 队列（华大基因员工健康队列）是本研究的对象。该队列在 DCS 平台上以公共容器数据集形式挂载于 /public/key_project_dataset/i99，包含三块数据：

本研究只使用第一块中的 taxonomic_profile.tsv（物种级相对丰度），因为上游已完成，无需重跑（原始 FASTQ 约 216 TB，重跑成本极高且结果等价）。

已知的队列特征

• 纵向但两段式：采样年份为 2016、2017、2018 与 2022、2023、2024、2025，2019–2021 完全无数据。
• 随访次数高度不平衡：抽样后 819 人仅 1 次采样，645 人 2 次，最多 10 次。
• 无表型元数据：队列不附带年龄、性别、BMI、用药、饮食等信息；已排查平台上的 CIMA 队列表型元数据（428 人），但其样本编号为 CIMA_H###，与 i99 的 HD317... 无法映射。
• 测序深度高度集中：中位 38.9M reads，IQR [38.5M, 39.1M]，但存在低尾（<10M 的 118 个 run，且 7.58% 集中在 2018 年）。

分析设计

鉴于队列规模（21,997 次 run），按用户要求采用 20% 个体级随机抽样：抽样单位为个体而非 run，以保住纵向随访结构；固定随机种子保证可复现；入样概率均匀（0.2）且不按年份补齐，以避免不等权重使总体估计向稀疏年份偏移。抽样代表性以个体级置换检验（1000 次重抽）验证，通过（perm_p=0.226）。

研究方案与方法

数据与抽样

• 输入：/public/key_project_dataset/i99/B2C_COM_META_TSV/{个体}/{年份}/{run}/{run}.taxonomic_profile.tsv
• 抽样：20% 个体级简单随机，seed=20260924，入样概率均匀 0.2，未触发年份兜底；实得 2,040 人 / 4,390 次 run（4,383 次有效，7 次空丰度被剔除）
• 代表性验证：个体级置换检验（1000 次重抽 20% 个体，比较各年 run 数的 χ² 距离），perm_p=0.226，各年 z∈[−1.74, +1.86] 通过
• 精度预估：Shannon 个体级 95% CI 半宽 0.0256，richness 2.60，ICC 0.0339

α 多样性指标

物种级 SGB 定义为 clade 深度 = 7（以 | 分隔）且末级以 s__ 开头，相对丰度在检出的 SGB 内重新归一化。计算：observed richness、Shannon H'=−Σp·ln p、Simpson 1−D、inverse Simpson、Pielou evenness J'、Hill 数 q=0/1/2。同时从文件注释行提取 reads processed 作为测序深度协变量。

统计方法（关键：以个体为独立单位）

由于同一人贡献多次 run，run 不独立。朴素 run 级推断会使置信区间过窄、p 值偏小，并过度代表高频随访者。因此：

执行环境与可复现性

• 在线容器 48 核 / 250 GB；python 3.12 + pandas 2.3 + scipy 1.18 + statsmodels 0.14 + matplotlib 3.11
• 全程只读 /public，不复制原始数据（仅读取约 0.5 GB）；6 个模块的总计算耗时 < 4 分钟
• 所有脚本位于 /work/xuxun/i99_alpha/scripts/，日志位于 /work/xuxun/i99_alpha/logs/；随机种子统一为 seed=20260924

质量把关

• 每个模块脚本在投递前经静态安全审计（0 critical/high）
• 每个模块完成后写 OAA（观察-评估-行动）与自评（dcs_self_review），其中 M3 经两轮、M5 经一轮修订后通过
• 交付前对 5 条核心主张做批判性评审（critical-review-expert 框架），并按评审意见推翻 1 条结论

结果

结果 1｜α 多样性分布显著左偏，人群中心值为 Shannon≈2.97、richness 中位 161 SGB

基于 20% 抽样队列（2,040 人 / 4,383 次有效 run）：Shannon 均值 2.970（个体级 2.944，聚类 bootstrap 95% CI [2.919, 2.968]）、中位 3.030、IQR [2.55, 3.49]、最大 4.71；richness 均值 166.4、中位 161、IQR [113, 210]、范围 [1, 469]；Simpson 1−D 0.857；Pielou evenness 0.589。

Shannon 显著左偏（skew=−0.709、kurtosis=1.069、Shapiro p=4.7e-28、D'Agostino p=9.7e-87），richness 右偏（+0.575）。该中心值落在已发表成人粪便宏基因组研究的常见区间（MetaPhlAn 物种级常检出 100–250 个 SGB、Shannon 多在 2.5–3.5），属中等水平。

需注意：大样本下正态性检验的 p 值极为敏感，且该检验未处理个体内重复相关性；因此「非正态」的证据强度应以外观偏度/峰度与分布图为准，结论仅限于本次抽样样本的描述性统计。

结果 2｜近窗口 α 多样性低于早窗口，方向稳健但幅度小、且与批次共线

个体随机截距混合模型给出 Shannon 的 late(2022–2025) − early(2016–2018) = −0.1482（SE 0.0228，95% CI [−0.193, −0.104]，p=7.4e-11）。在深度匹配子集（reads∈[38M,39M]，2,193 run / 1,417 人）中，richness 的 Δ=−43.13（p=2.3e-42）、Shannon 的 Δ=−0.2639（p=6.4e-18）。配对子队列（两窗口均有采样，n=321）Δ=−0.1306，Wilcoxon p=0.00116（182 人下降 / 139 人上升）。剔除低深度 run 后结论稳定。

必须同时呈现的限定条件：① 效应量普遍很小——21 对逐年比较中 Cliff's delta 为可忽略 9 对、小 9 对、中 2 对、大仅 1 对；② 年份固定效应只解释 Shannon 总方差的 3.7%（见结果 3）；③ 配对子队列中下降者仅占 57%，并非人人同向；④ 2019–2021 无数据，早/邻近为两个采样截面，与平台换代完全共线。因此结论表述为：「2022–2025 批次的 α 多样性低于 2016–2018 批次，方向在深度校正与配对分析中均稳健，但幅度小、异质性强，不能推断为连续时间趋势或个体层面的普遍衰退」。

结果 3｜变异以个体内波动为绝对主导：个体间仅解释约 11% 的方差

总体 ICC（个体级聚类 bootstrap 1000 次）= 0.1152（95% CI [0.0738, 0.1569]；MixedLM 复核 0.1101），重复性 R=0.2755。统一方差分解：个体间 11.0%、年份固定效应 3.7%、个体内残余 85.3%。个体内 SD 均值 0.5589、个体内 CV 均值 20.0%（95% CI [19.1, 21.0]，中位 16.6%）、个体内极差均值 1.004；个体均值 SD=0.5009（IQR 宽 0.67）；个体内/个体间 SD 比 = 1.12。

方法学修正（本项目最重要的自我纠错）：初步分析曾依据「段内早期 ICC=0.032 vs 段内近期 ICC=0.222」推断低 ICC 主要来自采样窗口/批次差异。批判性评审指出该归因方法学上不成立（不同子集的 ICC 在人数、重复次数、间隔、批次上均不可比，且早期 ICC 的置信区间跨 0）。补充的同个体内配对分析给出了相反方向的证据：在 285 名两窗口均有窗口内重复的人中，窗口内平均 |ΔShannon| = 0.651 反而大于跨窗口 |Δ| = 0.589（Wilcoxon p=0.020，仅 121/285 人跨窗口变化更大）。因此正确结论是：低 ICC 既不由长间隔解释、也不由窗口解释；85% 的方差属个体内残余波动，而在缺少提取批次/文库/机台元数据时，其中的生物学波动与技术噪声无法分离。

结果 4｜测序深度不是主要驱动因素，但低深度混杂使时期效应被低估约 2.4 倍

深度解释 richness 的 R²：全样本 0.0216 [0.0113, 0.0311]、剔除 <10M 后 0.0010、38–39M 子集内 0.0094；Spearman(深度, richness) 全样本 rho=−0.002（p=0.9）。解析式稀疏化到 10M reads：观测 169.0 → 稀疏化 169.0，收缩 0.0%（10M 以上 richness 已饱和）。window×depth 交互 LRT p=5.1e-64。

各方法方向完全一致（全部为负），但幅度差异显著：全样本未校正 Δrichness=−18.07，深度匹配子集为 −43.13（约 2.4 倍）。机理是 2018 年低深度样本（该年 7.58% 的 run <10M，其他年 0–1.2%）压低了早窗口的表观 richness。这不与稀疏化结论矛盾：稀疏化针对「把已测深的样本降到 10M」（饱和区无影响），而混杂来自「本来就低于 10M 的样本」（未饱和区 richness 被低估）。

重要保留：测序深度是后测变量而非随机分配变量。若低深度与低微生物载量、样本保存或个体用药状态相关，则按深度匹配/剔除将构成选择偏倚入口而非单纯的测量误差校正。故 −43.13 应表述为「深度匹配子集内的效应」，不外推为全队列真实效应量。

结果 5｜菌群组成整体稳定

属水平平均相对丰度 Top：Phocaeicola 13.95%、Bacteroides 9.54%、Segatella 7.24%、Prevotella 6.76%、Neisseria 2.69%、Faecalibacterium 1.86%、Roseburia 1.82%、Alistipes 1.70%、Lachnospira 1.60%。

逐年堆叠组成未见单一菌属的系统性替换或塌陷，提示结果 2 的群体级小幅下降更可能来自多个低丰度 SGB 的普遍窄幅变化，而非某个主导菌属失衡。需声明本图为 rel_ab 归一化值、非绝对定量，只能比较组成比例。

讨论

1. 主要发现的意义

「个体内波动主导」是本研究最稳健也最具方法学含义的结论。 个体身份只解释 Shannon 方差的 11.0%，而 85.3% 属个体内残余波动；个体内 CV 达 20%、个体内极差均值 1.004（Shannon 单位）。这意味着：用单次粪便样本推断某个体的「菌群多样性水平」在统计上是很弱的。对队列研究的直接启示是：若要用 α 多样性做暴露-结局关联，单次测量会造成严重的回归稀释（regression dilution），需要重复采样或误差校正。

值得强调的是，这个低可重复性不能归因于随访间隔长。同一批人内部的对照显示，窗口内相邻采样的变化幅度（0.651）反而大于跨 4 年窗口的变化（0.589，Wilcoxon p=0.020）；间隔与 |ΔShannon| 的秩相关仅 0.045。这对「长期菌群稳定性差」的直觉是一个反例。

2. 时期差异的解读边界

两个采样窗口间存在统计上非常显著的下降（Shannon p=7.4e-11，深度匹配后 richness p=2.3e-42），但显著性来自精度而非幅度：效应量多为小或可忽略，年份只解释 3.7% 的方差，配对子队列中下降者仅占 57%。更关键的是，2019–2021 年完全无数据，使「早期/近期」与「测序平台与流程换代」在结构上完全共线。因此本研究不主张存在个体层面的菌群多样性衰退趋势，只主张「两个采样批次之间存在一致性方向的差异」。要把它升级为生物学结论，至少需要：批次元数据、同批次内的纵向重复、以及表型协变量。

3. 技术因素的两种角色

测序深度在本队列呈现「高度集中（中位 38.9M）+ 低尾」的分布。分析显示深度总体上不是多样性的驱动因素：全样本秩相关近乎为零（rho=−0.002），剔除低尾后深度解释的 richness 方差仅 0.10%，且 10M reads 以上 richness 已饱和（稀疏化收缩 0.0%）。

但深度在本队列中扮演了两种不同角色：① 在低尾（2018 年 7.58% 的 run <10M）它是混杂因素，压低早窗口表观 richness，使时期效应被低估约 2.4 倍；② 同时它可能是选择偏倚的入口——低深度可能与低微生物载量、样本保存或个体用药状态相关，那么按深度匹配/剔除得到的效应只适用于「可成功测序的样本」。区分这两种角色需要额外的技术元数据与样本处理记录，本研究无法完成。

4. 与文献的关系

本队列的多样性水平（Shannon≈3.0、richness 中位 161 SGB）落在已发表成人粪便宏基因组研究的常见区间，属中等水平，未见异常。可重复性（ICC=0.115）低于文献常报的 0.2–0.5，一个重要的结构性解释是队列高度同质（华大员工在年龄、职业、地域、生活方式上接近），个体间方差被压缩会直接降低 ICC，这不等于个体的绝对变异更大。此外已有纵向研究报告两时间点间 Shannon 与检出物种数显著下降（Cell Reports 202501496-2)），本研究方向一致，但幅度小且与批次共线，不构成独立证据。

5. 本研究的自我纠错

值得记录的是，本研究在交付前推翻了自己的一个结论。初步分析曾依据「段内早期 ICC=0.032 vs 段内近期 ICC=0.222」推断低 ICC 主要来自采样窗口/批次差异；批判性评审指出该归因方法学不成立（不同子集的 ICC 在人数、重复次数、间隔、批次上均不可比，且早期 ICC 的置信区间跨 0）。补充的同个体内配对分析不仅未能支持该归因，反而给出相反方向的证据。相应结论已从「结论」降级为「待验证假设」。

6. 局限

1. 无表型元数据（年龄/性别/BMI/用药/饮食）——无法调整关键混杂，也无法解释时期差异的驱动机制；已核实 CIMA 元数据（428 人）因编号体系不同不可映射。

2. 无技术批次元数据（DNA 提取批次、文库、测序机台）——85.3% 的个体内残余方差无法拆分；时期效应与平台换代完全共线。

3. 两段式采样——2019–2021 无数据，只能做两窗口截面比较。

4. 相对丰度非绝对定量——无法判断菌群总载量是否变化。

5. 20% 抽样——代表性置换检验通过（p=0.226），但配对子队列（n=321）与早期段内 ICC（n=605）等分层结论的外推性受限。

6. 人群边界——i99 为健康、同质的华大员工队列，结论不可直接外推到一般人群或患者。

7. 下一步建议

按性价比排序：① 争取补齐表型与技术批次元数据（解除最大限制）；② 用全队列 21,997 run 复跑（脚本已保留 --full 开关）以消除抽样质疑；③ 复用 species_relab.parquet 做 β 多样性与差异丰度，检验表观下降是否伴随组成位移；④ 数据集中已具备同队列 HUMAnN3 功能谱（pathabundance、KO/GO/Pfam/eggnog），可做功能层面的 α 多样性与差异分析——这是本项目尚未使用的现成资源；⑤ 建立宏基因组与宿主 WGS（5,339 样本）的样本号映射后，可开展菌群-宿主遗传关联。

参考文献

数据来源

1. DCS Cloud 容器公共数据集：/public/key_project_dataset/i99/B2C_COM_META_TSV（i99 / 华大员工队列宏基因组 MetaPhlAn 4 + HUMAnN3 结果表）。

2. 方法学参考：MetaPhlAn 4 / mpa_vJun23_CHOCOPhlAnSGB_202403 参考库（SGB 级相对丰度）。

方法学与对照文献（本次交付中实际检索到的来源）

3. 测序深度对宏基因组多样性估计的影响（5–20M reads 评估，建议最低深度）：ScienceDirect

4. 纵向人群肠道菌群动态（瑞典一年期人群研究，个体内功能波动）：Cell Host & Microbe

5. 粪便标志物个体内变异（CV%）评估：PMC11490648

6. 两时间点间 Shannon 多样性与检出物种数显著下降的报道：Cell Reports 202501496-2)

7. 肠道微生物变异性与重复采样设计：Am J Clin Nutr00065-0)

> 说明：本报告的定量结论均来自本次对 i99 数据的实机计算，上列文献仅用于水平对标与方法学参照；未检索到的文献不做引用，亦不编造 PMID/DOI。