统计学基础（Statistics Fundamentals for Bioinformatics）

---

一句话说明

统计学是生信分析的"裁判"——你发现的差异菌群、筛选的标志物，到底是真差异还是运气好碰上的？统计方法帮你回答"这个结论靠不靠谱"，面试必考。

---

核心概念（白话版）

1. 描述性统计（Descriptive Statistics）

• 定义：用几个数字概括一堆数据的基本特征，包括集中趋势（数据堆在哪个位置附近）（均值 mean、中位数 median）、离散程度（数据散开的程度）（标准差 std、方差 variance、四分位距 IQR）等
• 白话比方：就像体检报告——身高、体重、血压这几个数字就能大致描述一个人的身体状况。描述性统计就是给你的数据做"体检"
• 关键点：
• 均值（mean）：所有数加起来除以个数，容易被极端值拉偏
• 中位数（median）：排好序取中间那个，不怕极端值
• 标准差（std）：衡量数据离均值有多远，越大说明数据越分散
• 四分位距（IQR）：Q3 - Q1，描述"中间 50% 的数据"分布范围
• 生信应用场景：
• 查看物种丰度表的基本分布：各菌属的平均丰度、标准差
• 样本测序深度的统计：均值、中位数判断测序是否均匀
• 基因表达量的分布概况

---

2. 假设检验（Hypothesis Testing）

假设检验的核心逻辑：先假设"没有差异"（零假设 H0），然后看数据是否能推翻这个假设。

2.1 t 检验（t-test）

• 定义：比较两组数据的均值是否有显著差异。前提假设：数据近似正态分布、方差齐性（equal variance）
• 白话比方：就像比较两个班的考试平均分——A 班 80 分 B 班 75 分，这 5 分差距到底是真的有差异还是偶然波动？t 检验帮你判断
• 关键点：
• 独立样本 t 检验（independent t-test）：两组不同的样本，比如健康人 vs 糖尿病患者
• 配对样本 t 检验（paired t-test）：同一组人治疗前 vs 治疗后
• 单样本 t 检验（one-sample t-test）：一组数据的均值是否等于某个已知值
• 使用前提：数据大致服从正态分布（用 Shapiro-Wilk 检验来判断）
• 生信应用场景：当样本量较大且数据近似正态时，比较健康组 vs T2D 组某个菌属的丰度差异

2.2 Wilcoxon 检验 / Mann-Whitney U 检验

• 定义：非参数检验（non-parametric test，非参数 = 不对数据分布形状做假设的检验方法，更"万金油"），不要求数据服从正态分布，基于数据的排名（rank）来比较两组差异
• 白话比方：不管你考了多少分，只看排名。A 班前 10 名里有 8 个比 B 班前 10 名分高，那 A 班就更强。不关心具体分数，只看排位
• 关键点：
• Wilcoxon 秩和检验（rank-sum test = Mann-Whitney U）：两组独立样本的非参数比较
• Wilcoxon 符号秩检验（signed-rank test）：配对样本的非参数比较
• 适用于：数据偏态分布（skewed，即数据不对称、一头长一头短，像个歪的山丘）、样本量小、有离群值（outlier）
• 生信中用得比 t 检验更多，因为微生物丰度数据通常不是正态分布
• 生信应用场景：T2D 肠道菌群项目中比较两组菌属丰度差异时，首选 Wilcoxon/Mann-Whitney U 检验

2.3 卡方检验（Chi-squared Test）

• 定义：检验两个分类变量（categorical variable）之间是否有关联，比较观察频数和期望频数（如果两件事没关系的话，理论上应该出现的次数）的差异
• 白话比方：统计吃辣和不吃辣的人里得胃病的比例，看"吃辣"和"得胃病"之间有没有关系
• 关键点：
• 输入是频次/计数，不是连续数值
• 要求每个格子的期望频次 >= 5（太少就用 Fisher 精确检验）
• 代码示例：

  # Python：卡方检验
  from scipy.stats import chi2_contingency
  # 构建列联表（2x2为例：吃辣/不吃辣 × 得胃病/没得）
  observed = [[30, 10], [20, 40]]  # 观察频数
  chi2, p, dof, expected = chi2_contingency(observed)
  print(f"卡方值={chi2:.2f}, p={p:.4f}, 自由度={dof}")
  print(f"期望频数:\n{expected}")
  

  # R：卡方检验
  observed <- matrix(c(30, 20, 10, 40), nrow = 2)
  result <- chisq.test(observed)
  print(result)           # 输出卡方值、自由度、p值
  print(result$expected)  # 查看期望频数
  

• 生信应用场景：比较不同疾病组中某个基因型的频率分布；检验菌群"有/无"在病例/对照中的分布差异

2.4 ANOVA（方差分析，Analysis of Variance）

• 定义：比较 3 组及以上 数据均值是否有显著差异的参数检验方法
• 白话比方：t 检验是两个班比平均分，ANOVA 是全年级 5 个班一起比——看看这些班的平均分有没有真正的差别，还是大家差不多
• 关键点：
• 使用前提：数据近似正态分布 + 各组方差齐性（方差差不多大）。不满足就换 Kruskal-Wallis
• 如果 ANOVA 显著（p < 0.05），只说明"至少有一组不一样"，想知道具体哪两组有差就要做事后检验（post-hoc test，如 Tukey HSD）
• 和 Kruskal-Wallis 的关系：ANOVA 是参数版（假设正态），Kruskal-Wallis 是非参数版（不假设正态）。数据正态用 ANOVA，不正态用 Kruskal-Wallis
• 代码示例：

  # Python：单因素 ANOVA
  from scipy.stats import f_oneway
  group1 = [23, 25, 28, 22, 27]
  group2 = [30, 33, 29, 35, 31]
  group3 = [18, 20, 22, 19, 21]
  f_stat, p_value = f_oneway(group1, group2, group3)
  print(f"F统计量={f_stat:.2f}, p={p_value:.4f}")
  

  # R：单因素 ANOVA
  values <- c(23,25,28,22,27, 30,33,29,35,31, 18,20,22,19,21)
  groups <- factor(rep(c("A","B","C"), each = 5))
  result <- aov(values ~ groups)
  summary(result)
  # 事后检验：看具体哪两组有差
  TukeyHSD(result)
  

• 生信应用场景：当样本分组 >= 3 组且数据近似正态时（如基因表达量经 log 变换后），比较多组间差异

2.5 Kruskal-Wallis 检验

• 定义：Wilcoxon 检验的"多组版"——当有 3 组及以上时，用 Kruskal-Wallis 检验比较它们之间是否有差异。它是 ANOVA 的非参数替代（不要求正态分布）
• 白话比方：Wilcoxon 是两个人比赛，Kruskal-Wallis 是一场锦标赛看所有选手有没有水平差异
• 生信应用场景：比较轻度/中度/重度 T2D 患者三组间的菌群丰度差异（菌群丰度不是正态的，所以用 Kruskal-Wallis 而不是 ANOVA）

---

3. 多重检验校正（Multiple Testing Correction）

• 定义：当你同时对很多个基因/菌属做统计检验时，需要对 p 值进行校正，降低假阳性（false positive）的风险
• 白话比方：你扔硬币 1 次出正面不稀奇（p=0.5），但如果你扔 1000 次，肯定会有连续出 10 次正面的情况——这不是硬币有问题，是你试的次数太多了。多重校正就是解决"试太多次"导致的虚假发现
• 关键点：
• 为什么要校正：检测 1000 个菌属时，即使每个都没差异，按 p < 0.05 的标准也会有约 50 个"假差异"（0.05 × 1000 = 50）
• Bonferroni 校正：最严格，p 值乘以检验次数。比如检测 1000 个菌，原来的 p 要乘以 1000。优点是严格控制假阳性；缺点是太保守，很多真差异也被漏掉了
• FDR / BH 校正（Benjamini-Hochberg）：控制的是"你报出来的显著结果里，假的占多大比例"。比 Bonferroni 宽松合理，是生信中最常用的校正方法
• 校正后的 p 值通常叫 adjusted p-value 或 q-value 或 FDR
• 生信应用场景：
• 差异菌群分析：上百个菌同时做 Wilcoxon 检验，必须做 FDR 校正
• 差异表达基因（DEG）分析：几万个基因同时检验，FDR < 0.05 是标准阈值
• 你的 T2D 项目中 run_de_analysis.py 就用了 scipy.stats.false_discovery_control(method='bh')（需要 scipy >= 1.11.0）

---

4. 相关性分析（Correlation Analysis）

4.1 Pearson 相关系数

• 定义：衡量两个变量之间的线性相关程度，取值 -1 到 1
• 白话比方：身高和体重一般一起增长，这就是正相关（r 接近 1）；温度和穿衣厚度反着来，就是负相关（r 接近 -1）；学号和成绩没关系，r 接近 0
• 关键点：
• r = 1 完全正相关，r = -1 完全负相关，r = 0 无线性关系
• 前提：两个变量都近似正态分布，关系是线性的
• 对离群值敏感

4.2 Spearman 相关系数

• 定义：基于排名的相关系数，衡量两个变量的单调关系（不一定是直线，只要一个增另一个也增就行）
• 白话比方：Pearson 只认直线关系，Spearman 更灵活——只要趋势一致就算相关，不管是直线还是曲线
• 关键点：
• 不要求正态分布，不受离群值影响
• 生信中比 Pearson 更常用，因为微生物丰度数据通常是非正态的
• 生信应用场景：
• 分析菌属丰度与临床指标（如 HbA1c 糖化血红蛋白、BMI）的相关性
• 菌属间的共现关系（co-occurrence）分析

---

5. 降维分析（Dimensionality Reduction）

5.1 PCA（Principal Component Analysis，主成分分析）

• 定义：将高维数据（比如几百个菌属的丰度）压缩到 2-3 个维度来可视化，同时尽量保留数据的主要差异信息
• 白话比方：你的简历有几十项指标（成绩、项目经历、技能等），HR 面试时只看"综合实力"和"发展潜力"两个维度打分——这就是降维，把复杂信息压缩成少数关键指标
• 关键点：
• 基于原始数据的方差做数学变换（把原来的多个指标重新组合成新的"综合指标"）
• PC1（第一主成分）捕获最大的变异方向，PC2 捕获第二大的
• 适用于欧氏距离（Euclidean distance）适用的场景（比如基因表达量）
• 图上离得近的样本 = 性质相似

5.2 PCoA（Principal Coordinates Analysis，主坐标分析）

• 定义：和 PCA 类似的降维可视化，但输入是距离矩阵而不是原始数据，可以使用多种距离度量
• 白话比方：PCA 是拿原始数据直接降维（像看全班的原始成绩单），PCoA 是先算好任意两个人"有多不同"（距离矩阵），再画图。因为距离的算法可以换，所以更灵活
• 关键点：
• 微生物生态学中常用 Bray-Curtis 距离（基于物种组成的不相似性指数，non-metric dissimilarity，不满足三角不等式所以不是严格的"距离"而是"不相似度"）。白话说：只看共有物种的相似程度，不像欧氏距离那样把"都没有"也算进去
• PCA vs PCoA 核心区别：PCA 直接用原始数据 + 欧氏距离；PCoA 可以选择任意距离度量（Bray-Curtis、UniFrac 等）
• 微生物组研究几乎都用 PCoA（因为 Bray-Curtis 距离比欧氏距离更适合物种丰度数据）
• 生信应用场景：
• 你的 T2D 项目中可视化健康组 vs T2D 组的菌群结构差异
• PERMANOVA 检验（用 vegan::adonis2）配合 PCoA 图展示组间差异是否显著

---

6. 正态分布与非正态分布

• 定义：
• 正态分布（Normal Distribution）：钟形曲线，数据对称地分布在均值两侧，大部分数据集中在均值附近。也叫高斯分布（Gaussian Distribution）
• 非正态分布：不满足上述特征的分布，如右偏分布（right-skewed）、双峰分布等
• 白话比方：全班考试成绩如果大多数人在 70-80 分之间，特别高和特别低的很少，画出来像一个钟——这就是正态分布。但微生物丰度通常是少数菌占大量、多数菌很少，像一个歪的山——这就是非正态
• 关键点：
• 为什么重要：很多经典统计方法（t 检验、Pearson 相关、ANOVA）都假设数据是正态分布的
• 怎么判断：
  • Shapiro-Wilk 检验（scipy.stats.shapiro），p > 0.05 = 可以认为是正态的。最常用，适合小样本（n < 50）
  • Q-Q 图（Quantile-Quantile Plot）：把数据的分位数和理论正态分布的分位数画在一张散点图上，如果数据是正态的，点会排成一条直线。偏离直线越多 = 越不正态。Python 用 scipy.stats.probplot(data, plot=plt)，R 用 qqnorm(data); qqline(data)
  • K-S 检验（Kolmogorov-Smirnov）（scipy.stats.kstest(data, 'norm', args=(mean, std))）：比较数据分布和理论正态分布的差距，p > 0.05 = 可认为正态。适合大样本（n > 50），小样本不敏感
• 不满足正态怎么办：用非参数检验（Wilcoxon、Spearman、Kruskal-Wallis），或者对数据做变换（log 变换、CLR 变换）
• 微生物丰度数据几乎都是非正态的，所以生信中非参数方法是主流
• 生信应用场景：选择统计方法前，先检查数据分布决定用参数检验还是非参数检验

---

7. p 值与效应量（p-value & Effect Size）

• 定义：
• p 值：假设"没有差异"（H0）是真的情况下，观察到当前数据（或更极端情况）的概率。p 值越小，说明数据越不支持"没有差异"
• 效应量（Effect Size）：衡量差异的"实际大小"，比如 Cohen's d（把两组均值差除以标准差，得到一个"标准化后的差距大小"，0.2=小差异、0.5=中等、0.8=大差异）、log2FC
• 白话比方：p 值告诉你"这个差异是不是真的"，效应量告诉你"这个差异有多大"。就像法官判案——p 值是定罪（有没有犯罪），效应量是量刑（罪有多重）。一个 p = 0.001 但 log2FC = 0.01 的基因，虽然"确实有差异"但差异极其微小，生物学上可能没意义
• 关键点：
• p < 0.05 是常用阈值，但不是绝对标准
• p 值 ≠ 效应量：大样本量下，极小的差异也能得到很小的 p 值
• 差异基因/差异菌群的筛选通常要求双重标准：p 值（或 FDR）< 0.05 且 |log2FC| > 1
• 你的项目中用的就是 |log2FC| > 1 且 FDR < 0.05 的标准
• 生信应用场景：火山图（Volcano Plot）的 x 轴是 log2FC（效应量），y 轴是 -log10(p-value)，同时展示统计显著性和生物学效应

---

常用工具/命令

工具/库	语言	用途	安装方式
scipy.stats	Python	t 检验、Wilcoxon、Shapiro-Wilk、Spearman、FDR 校正等	pip install scipy
statsmodels	Python	回归分析、多重检验校正（multipletests）	pip install statsmodels
scikit-learn	Python	PCA 降维、随机森林、交叉验证	pip install scikit-learn
numpy	Python	数组运算、基础数学计算	pip install numpy
pandas	Python	数据读取、整理、分组操作	pip install pandas
scikit-bio	Python	Bray-Curtis 距离、PCoA、Shannon 指数	pip install scikit-bio
stats（内置）	R	t.test、wilcox.test、chisq.test、cor.test、p.adjust	R 内置，无需安装
vegan	R	多样性分析、PERMANOVA（adonis2）、PCoA	install.packages("vegan")
ggplot2	R	数据可视化（箱线图、PCA 图等）	install.packages("ggplot2")

---

实操代码

Python 实操

#!/usr/bin/env python3
"""
统计学基础实操 —— Python 版
包含：描述性统计、正态性检验、t检验、Wilcoxon检验、
      Spearman相关、FDR校正、PCA降维
"""

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# ============================================================
# 1. 模拟数据：健康组 vs T2D组 的 5 个菌属丰度
# ============================================================
np.random.seed(42)  # 固定随机种子，保证每次运行结果一致

# 模拟 10 个健康人的菌属丰度（单位：相对丰度 %）
healthy = pd.DataFrame({
    'Bacteroides':    np.random.lognormal(2.0, 0.5, 10),  # 拟杆菌属
    'Lactobacillus':  np.random.lognormal(1.5, 0.6, 10),  # 乳杆菌属
    'Faecalibacterium': np.random.lognormal(1.8, 0.4, 10),  # 粪杆菌属（产短链脂肪酸）
    'Escherichia':    np.random.lognormal(0.5, 0.8, 10),  # 大肠杆菌属
    'Prevotella':     np.random.lognormal(1.2, 0.7, 10),  # 普雷沃氏菌属
})

# 模拟 10 个 T2D 患者的菌属丰度
# T2D 特征：Faecalibacterium 下降，Escherichia 上升
t2d = pd.DataFrame({
    'Bacteroides':    np.random.lognormal(2.1, 0.5, 10),
    'Lactobacillus':  np.random.lognormal(1.3, 0.6, 10),
    'Faecalibacterium': np.random.lognormal(1.0, 0.4, 10),  # 明显降低
    'Escherichia':    np.random.lognormal(1.5, 0.8, 10),    # 明显升高
    'Prevotella':     np.random.lognormal(1.1, 0.7, 10),
})

print("=" * 60)
print("1. 描述性统计")
print("=" * 60)
print("\n健康组：")
print(healthy.describe().round(2))  # describe() 一键输出均值、标准差、四分位数等
print("\nT2D 组：")
print(t2d.describe().round(2))

# ============================================================
# 2. 正态性检验（Shapiro-Wilk）
# ============================================================
print("\n" + "=" * 60)
print("2. 正态性检验（Shapiro-Wilk）")
print("  p > 0.05 → 可以认为是正态分布")
print("  p < 0.05 → 非正态分布，要用非参数检验")
print("=" * 60)

for col in healthy.columns:
    stat_h, p_h = stats.shapiro(healthy[col])  # 检验健康组
    stat_t, p_t = stats.shapiro(t2d[col])      # 检验 T2D 组
    normal_h = "正态" if p_h > 0.05 else "非正态"
    normal_t = "正态" if p_t > 0.05 else "非正态"
    print(f"  {col:20s}  健康组 p={p_h:.4f}({normal_h})  T2D组 p={p_t:.4f}({normal_t})")

# ============================================================
# 3. 假设检验：t 检验 vs Wilcoxon 检验
# ============================================================
print("\n" + "=" * 60)
print("3. 两组比较：t 检验 & Wilcoxon 检验")
print("=" * 60)

pvalues_t = []       # 存储 t 检验的 p 值
pvalues_wilcox = []  # 存储 Wilcoxon 检验的 p 值

for col in healthy.columns:
    # t 检验（假设正态分布）
    t_stat, p_t = stats.ttest_ind(healthy[col], t2d[col])

    # Wilcoxon 秩和检验 = Mann-Whitney U 检验（不假设正态）
    u_stat, p_w = stats.mannwhitneyu(healthy[col], t2d[col], alternative='two-sided')

    pvalues_t.append(p_t)
    pvalues_wilcox.append(p_w)

    print(f"  {col:20s}  t检验 p={p_t:.4f}  Wilcoxon p={p_w:.4f}")

# ============================================================
# 4. 多重检验校正（FDR / Benjamini-Hochberg）
# ============================================================
print("\n" + "=" * 60)
print("4. FDR 校正（Benjamini-Hochberg）")
print("  校正后 FDR < 0.05 = 显著差异菌属")
print("=" * 60)

# 方法一：用 scipy（需要 scipy >= 1.11.0 才支持 false_discovery_control）
fdr_scipy = stats.false_discovery_control(pvalues_wilcox, method='bh')

# 方法二：用 statsmodels（更经典的方式）
from statsmodels.stats.multitest import multipletests
reject, fdr_sm, _, _ = multipletests(pvalues_wilcox, alpha=0.05, method='fdr_bh')
# reject: 是否拒绝 H0（True/False）
# fdr_sm: 校正后的 p 值

for i, col in enumerate(healthy.columns):
    sig = "← 显著" if fdr_sm[i] < 0.05 else ""
    print(f"  {col:20s}  原始p={pvalues_wilcox[i]:.4f}  FDR={fdr_sm[i]:.4f} {sig}")

# ============================================================
# 5. 相关性分析（Spearman）
# ============================================================
print("\n" + "=" * 60)
print("5. Spearman 相关性（以健康组为例）")
print("  |rho| > 0.6 = 强相关")
print("=" * 60)

# 计算 Bacteroides 和 Faecalibacterium 的相关性
rho, p_corr = stats.spearmanr(healthy['Bacteroides'], healthy['Faecalibacterium'])
print(f"  Bacteroides vs Faecalibacterium: rho={rho:.3f}, p={p_corr:.4f}")

# 计算整个矩阵的 Spearman 相关系数
corr_matrix = healthy.corr(method='spearman')  # pandas 内置方法
print("\n  Spearman 相关矩阵：")
print(corr_matrix.round(2).to_string())

# ============================================================
# 6. PCA 降维可视化
# ============================================================
print("\n" + "=" * 60)
print("6. PCA 降维")
print("=" * 60)

# 合并两组数据
all_data = pd.concat([healthy, t2d], ignore_index=True)
labels = ['Healthy'] * 10 + ['T2D'] * 10  # 分组标签

# 标准化（PCA 前必须做，让每个特征的均值=0，方差=1）
scaler = StandardScaler()
data_scaled = scaler.fit_transform(all_data)

# 执行 PCA，降到 2 维
pca = PCA(n_components=2)  # 保留前 2 个主成分
pc = pca.fit_transform(data_scaled)

# 查看每个主成分解释了多少方差
print(f"  PC1 解释方差比例：{pca.explained_variance_ratio_[0]:.1%}")
print(f"  PC2 解释方差比例：{pca.explained_variance_ratio_[1]:.1%}")
print(f"  合计解释：{sum(pca.explained_variance_ratio_):.1%}")

# 输出 PCA 坐标（可用于画图）
pca_df = pd.DataFrame({
    'PC1': pc[:, 0],
    'PC2': pc[:, 1],
    'Group': labels
})
print("\n  PCA 坐标（前5行）：")
print(pca_df.head().to_string(index=False))

---

R 实操

#!/usr/bin/env Rscript
# 统计学基础实操 —— R 语言版
# 包含：描述性统计、正态性检验、t检验、Wilcoxon检验、
#       Spearman相关、FDR校正、PCoA降维

# ============================================================
# 1. 模拟数据
# ============================================================
set.seed(42)  # 固定随机种子

# 健康组 10 个样本，5 个菌属的丰度
healthy <- data.frame(
  Bacteroides       = rlnorm(10, meanlog = 2.0, sdlog = 0.5),
  Lactobacillus     = rlnorm(10, meanlog = 1.5, sdlog = 0.6),
  Faecalibacterium  = rlnorm(10, meanlog = 1.8, sdlog = 0.4),
  Escherichia       = rlnorm(10, meanlog = 0.5, sdlog = 0.8),
  Prevotella        = rlnorm(10, meanlog = 1.2, sdlog = 0.7)
)

# T2D 组 10 个样本
t2d <- data.frame(
  Bacteroides       = rlnorm(10, meanlog = 2.1, sdlog = 0.5),
  Lactobacillus     = rlnorm(10, meanlog = 1.3, sdlog = 0.6),
  Faecalibacterium  = rlnorm(10, meanlog = 1.0, sdlog = 0.4),  # 明显降低
  Escherichia       = rlnorm(10, meanlog = 1.5, sdlog = 0.8),   # 明显升高
  Prevotella        = rlnorm(10, meanlog = 1.1, sdlog = 0.7)
)

cat("========== 1. 描述性统计 ==========\n")
cat("健康组：\n")
print(summary(healthy))   # summary() 输出最小值、四分位数、均值、最大值
cat("\nT2D 组：\n")
print(summary(t2d))

# ============================================================
# 2. 正态性检验（Shapiro-Wilk）
# ============================================================
cat("\n========== 2. 正态性检验 ==========\n")
for (col in colnames(healthy)) {
  p_h <- shapiro.test(healthy[[col]])$p.value  # 健康组
  p_t <- shapiro.test(t2d[[col]])$p.value      # T2D 组
  cat(sprintf("  %20s  健康组 p=%.4f  T2D组 p=%.4f\n", col, p_h, p_t))
}

# ============================================================
# 3. 假设检验
# ============================================================
cat("\n========== 3. 两组比较 ==========\n")
p_values <- c()  # 存放原始 p 值

for (col in colnames(healthy)) {
  # t 检验
  p_t <- t.test(healthy[[col]], t2d[[col]])$p.value

  # Wilcoxon 秩和检验（非参数）
  p_w <- wilcox.test(healthy[[col]], t2d[[col]])$p.value
  p_values <- c(p_values, p_w)

  cat(sprintf("  %20s  t检验 p=%.4f  Wilcoxon p=%.4f\n", col, p_t, p_w))
}

# ============================================================
# 4. FDR 校正
# ============================================================
cat("\n========== 4. FDR 校正 ==========\n")
# p.adjust() 是 R 内置函数，method = "BH" 就是 Benjamini-Hochberg
fdr <- p.adjust(p_values, method = "BH")

for (i in seq_along(colnames(healthy))) {
  sig <- ifelse(fdr[i] < 0.05, "<- 显著", "")
  cat(sprintf("  %20s  原始p=%.4f  FDR=%.4f %s\n",
              colnames(healthy)[i], p_values[i], fdr[i], sig))
}

# ============================================================
# 5. Spearman 相关
# ============================================================
cat("\n========== 5. Spearman 相关 ==========\n")
cor_result <- cor.test(healthy$Bacteroides, healthy$Faecalibacterium,
                       method = "spearman")
cat(sprintf("  Bacteroides vs Faecalibacterium: rho=%.3f, p=%.4f\n",
            cor_result$estimate, cor_result$p.value))

# ============================================================
# 6. PCoA（主坐标分析）
# ============================================================
cat("\n========== 6. PCoA ==========\n")

# 如果没装 vegan 包，先安装：install.packages("vegan")
library(vegan)

# 合并数据
all_data <- rbind(healthy, t2d)
group <- c(rep("Healthy", 10), rep("T2D", 10))

# 计算 Bray-Curtis 距离矩阵
bc_dist <- vegdist(all_data, method = "bray")

# 执行 PCoA
pcoa_result <- cmdscale(bc_dist, k = 2, eig = TRUE)
# k=2 表示降到 2 维；eig=TRUE 返回特征值用于计算方差解释比

# 计算每个轴解释的方差比例
eig_vals <- pcoa_result$eig
var_explained <- eig_vals / sum(abs(eig_vals)) * 100

cat(sprintf("  PCoA1 解释方差：%.1f%%\n", var_explained[1]))
cat(sprintf("  PCoA2 解释方差：%.1f%%\n", var_explained[2]))

# PERMANOVA 检验：组间差异是否显著
perm <- adonis2(bc_dist ~ group, permutations = 999)
cat(sprintf("  PERMANOVA p-value: %.4f\n", perm$`Pr(>F)`[1]))
cat("  （p < 0.05 说明两组菌群结构有显著差异）\n")

---

和你项目的关联

你的简历项目涉及 2 型糖尿病（T2D）肠道菌群宏基因组分析 + 随机森林分类，以下是具体用到统计学的地方：

项目环节	用到的统计方法	具体做了什么
差异菌群筛选	Wilcoxon 检验 + FDR 校正	逐个菌属比较 T2D 组和健康组的丰度，用 BH 方法校正 p 值，FDR < 0.05 的菌属才算显著差异
效应量评估	log2 Fold Change	计算每个菌属在两组间的丰度变化倍数，
Alpha 多样性	Shannon 指数	计算每个样本的物种多样性，T2D 患者通常 Shannon 指数更低（多样性下降）
Beta 多样性	Bray-Curtis + PCoA	用 Bray-Curtis 距离计算样本间的菌群组成差异，PCoA 可视化，PERMANOVA 检验组间差异显著性
相关性分析	Spearman 相关	分析核心菌属（如 Faecalibacterium）与临床指标（HbA1c、BMI）的相关性
随机森林建模前	描述性统计 + 特征筛选	查看各菌属丰度分布，剔除低丰度菌（prevalence filter，出现率 < 5% 的菌属不纳入模型）
模型评估	交叉验证 + 置换重要性	5 折分层交叉验证评估模型稳定性，permutation importance 筛选核心菌属

你的代码中的实际应用： - project/scripts/run_de_analysis.py：用 scipy.stats.ttest_ind 做 t 检验，用 scipy.stats.false_discovery_control(method='bh') 做 FDR 校正 - project/scripts/calc_diversity.py：手动实现 Shannon 指数和 Bray-Curtis 距离计算 - t2d_ai_project/scripts/train_rf_prjna686835_thesis.py：用 sklearn 的 StratifiedKFold 做 5 折交叉验证，permutation_importance 做特征重要性评估

---

面试怎么答

Q1: 什么是 p 值？怎么理解？

p 值就是"假设没有差异的情况下，出现当前结果的概率"。比如我比较健康组和 T2D 组的某个菌属丰度，如果 p = 0.01，意思是"如果这两组真的没差异，那出现现在这么大差距的概率只有 1%"——概率这么小，我就有理由认为它们确实有差异。通常用 p < 0.05 作为显著性的门槛。

但 p 值有个容易误解的点：p 值不是"差异存在的概率"，它只是一个证据强度的指标。而且 p 值不能告诉你差异有多大，所以在生信分析中我们通常还要结合 log2FC 这样的效应量指标一起看。

Q2: 什么时候用 t 检验，什么时候用 Wilcoxon 检验？

关键看数据分布。如果数据近似正态分布，用 t 检验；如果数据是偏态的、样本量又小、或者有离群值，用 Wilcoxon 检验更稳妥。

实际在生信中，微生物丰度数据基本都是非正态的——少数菌很多、多数菌很少，分布很偏。所以我在做 T2D 菌群差异分析时主要用 Wilcoxon 秩和检验。判断方法也很简单：先用 Shapiro-Wilk 检验看一下正态性，p < 0.05 就说明不是正态的，那就用 Wilcoxon。

Q3: 为什么要做多重检验校正？FDR 和 Bonferroni 有什么区别？

因为同时测很多个基因或菌属，假阳性会累积。打个比方，如果你只测 1 个菌有没有差异，按 p < 0.05 的标准，犯错概率是 5%。但如果你同时测 1000 个菌，就算它们全都没差异，也会有大约 50 个被"误判"成有差异的（1000 × 0.05 = 50）。多重校正就是解决这个问题。

Bonferroni 校正最简单粗暴：把 p 值乘以检验次数。1000 个检验的话，原来 p = 0.001 的校正后变成 1。这样假阳性控制得很严，但太保守，很多真差异也被漏掉了。

FDR（BH 方法）更聪明：它控制的是"你报出来的显著结果里，假的占多大比例"。比如你筛出 100 个差异菌属，FDR = 0.05 意味着其中大概 5 个是假阳性，95 个是真的。这个比 Bonferroni 合理得多，所以在生信中 FDR 是标准做法。

Q4: PCA 和 PCoA 有什么区别？

两个都是降维可视化的方法，目的一样——把高维数据画到二维图上看样本分布。

PCA 是直接拿原始数据做线性变换，默认用的是欧氏距离。适合基因表达量这种连续的、近似正态的数据。

PCoA 是先算好样本间的距离矩阵，然后基于距离矩阵做降维。好处是距离的算法可以自己选——微生物组研究通常用 Bray-Curtis 距离，因为它比欧氏距离更适合处理物种丰度数据（考虑了"你有我没有"的情况）。

所以在我的 T2D 肠道菌群项目中用的是 PCoA + Bray-Curtis 距离，不是 PCA。这也是微生物组领域的标准做法。

Q5: 你的项目中用了哪些统计方法？

我的项目是 2 型糖尿病肠道菌群的宏基因组分析，统计方法贯穿了整个流程：

首先，做差异菌群分析时，因为菌属丰度是非正态的，我用了 Wilcoxon 秩和检验 逐个比较健康组和 T2D 组，然后用 BH 方法做 FDR 校正，同时要求 |log2FC| > 1，双重标准筛选差异菌属。

然后，多样性分析用了 Shannon 指数 看 alpha 多样性，用 Bray-Curtis + PCoA 看 beta 多样性，还用 PERMANOVA 检验两组菌群结构的整体差异是否显著。

相关性分析方面，我用 Spearman 相关 分析核心菌属和临床指标的关系，因为数据不是正态的。

最后建模阶段，我用了 随机森林，通过 5 折分层交叉验证 评估模型性能，用 permutation importance 排序特征重要性来确定核心菌属。在特征筛选之前还做了 prevalence filter，出现率低于 5% 的菌属不纳入模型。

---

延伸阅读

1. 《生物统计学》（李春喜 等） —— 中文教材，适合入门。讲解常用统计方法和生物学应用，语言友好
2. StatQuest YouTube 频道（Josh Starmer） —— 英文但有字幕，用动画讲统计和机器学习，通俗易懂。搜索 "StatQuest t-test"、"StatQuest PCA"、"StatQuest FDR" 即可（https://www.youtube.com/@statquest）
3. SciPy 官方统计文档 —— Python 统计函数的权威参考，包含所有检验方法的用法和示例（https://docs.scipy.org/doc/scipy/reference/stats.html）