统计学习路线（生信方向）¶

# === 均值、中位数、众数 ===
expression <- c(2.3, 5.1, 3.8, 45.2, 4.5, 3.9, 5.0, 4.2)  # 基因表达量

mean(expression)     # 均值：9.25（被离群值45.2拉高了！）
median(expression)   # 中位数：4.35（不受离群值影响，更"靠谱"）

# 白话解释：
# 均值 = 所有数加起来除以个数（容易被极端值带偏）
# 中位数 = 排好序后最中间的数（对离群值免疫）
# 生信中：表达量数据通常右偏（少数基因极高表达），所以常看中位数

# === 方差、标准差、四分位距 ===
var(expression)      # 方差：样本方差（偏差的平方的平均）
sd(expression)       # 标准差：方差开根号（和原数据同单位）
IQR(expression)      # 四分位距：Q3-Q1（中间50%数据的范围）

# 白话解释：
# 方差/标准差 = 数据有多"散"（越大越散）
# 标准差的好处 = 单位和原数据一样（方差的单位是"平方"，不直观）
# IQR = 不受离群值影响的"散度"指标

# === CV（变异系数）—— 比较不同量纲数据的离散度 ===
cv <- sd(expression) / mean(expression) * 100  # CV = 标准差/均值 × 100%
# 白话：相对于均值，数据波动了百分之多少
# 生信中：CV常用来评估技术重复之间的一致性

# === 正态分布（高斯分布）===
# 白话：大部分数据集中在中间，越往两边越少，像一口钟
# 性质：均值±1SD包含68%数据，±2SD包含95%，±3SD包含99.7%

# 检验数据是否正态分布
shapiro.test(expression)    # Shapiro-Wilk检验（p>0.05表示不拒绝正态假设）

# === log转换 ===
# 生信中基因表达数据通常是右偏的（少数基因表达量极高）
# log转换后更接近正态分布，更适合做统计检验
log2_expr <- log2(expression + 1)  # +1避免log(0)

# === 可视化分布 ===
par(mfrow = c(1, 2))               # 一行两列画图
hist(expression, main = "原始数据", breaks = 20)      # 直方图
hist(log2_expr, main = "log2转换后", breaks = 20)     # log2后的直方图

# QQ图（判断正态性的可视化方法）
qqnorm(log2_expr)       # 画QQ图
qqline(log2_expr)       # 加参考线（点在线上 = 接近正态）

# === 核心概念 ===
# H0（零假设）：两组没差异（"无罪推定"）
# H1（备择假设）：两组有差异（"有罪"）
# p值：在H0成立的情况下，观察到当前结果（或更极端）的概率
# p < 0.05：拒绝H0，认为差异"统计显著"
# p >= 0.05：不能拒绝H0（注意：不是"证明没差异"！）

# 白话类比：
# 你怀疑硬币被动了手脚（H1），H0是"硬币公平"
# 你扔了100次，正面90次，这个概率极低（p非常小）
# 所以你拒绝"硬币公平"的假设，认为硬币有问题

# === 两类错误 ===
# I类错误（假阳性）：实际没差异，你说有差异（冤枉好人）
# II类错误（假阴性）：实际有差异，你说没差异（放过坏人）
# alpha（显著性水平）：允许的I类错误率，通常0.05
# beta：II类错误率
# power（统计功效）= 1 - beta：检测到真实差异的能力

# === t检验（比较两组均值）===
# 前提：数据近似正态分布，方差齐性
wt_expr <- c(12.3, 15.6, 14.1, 13.8, 15.2)    # WT组表达量
ko_expr <- c(3.2, 2.8, 4.1, 3.5, 2.9)          # KO组表达量

# 独立样本t检验（两组独立样本）
t.test(wt_expr, ko_expr)
# 白话：比较WT和KO的平均表达量是否有显著差异

# 配对t检验（同一样本处理前后）
t.test(before, after, paired = TRUE)

# === Wilcoxon检验（非参数，不要求正态分布）===
# 白话：数据不是正态分布时，用这个代替t检验
wilcox.test(wt_expr, ko_expr)
# 生信中：样本量小或数据明显偏斜时，优先用非参数检验

# === 卡方检验（比较比例/频数）===
# 例：某突变在病人和正常人中的出现频率是否不同？
observed <- matrix(c(30, 70, 10, 90), nrow = 2)  # 2x2列联表
# 病人组：30有突变，70无突变
# 正常组：10有突变，90无突变
chisq.test(observed)

# === Fisher精确检验（小样本时代替卡方检验）===
fisher.test(observed)
# 白话：样本少的时候卡方检验不准，用Fisher更靠谱

# === ANOVA（方差分析，比较3组及以上的均值）===
# 白话：t检验只能比两组，ANOVA能比多组
group <- factor(c(rep("A", 5), rep("B", 5), rep("C", 5)))  # 3组
values <- c(rnorm(5, 10), rnorm(5, 12), rnorm(5, 15))      # 模拟数据
anova_result <- aov(values ~ group)
summary(anova_result)

# 事后检验（ANOVA显著后，看具体哪两组之间有差异）
TukeyHSD(anova_result)

你要比较什么？
├── 两组均值
│   ├── 数据正态？ → 是 → t检验（独立/配对）
│   └── 数据不正态？ → Wilcoxon检验
├── 三组及以上均值
│   ├── 数据正态？ → 是 → ANOVA + 事后检验
│   └── 数据不正态？ → Kruskal-Wallis检验
├── 比例/频数
│   ├── 样本量大？ → 卡方检验
│   └── 样本量小？ → Fisher精确检验
└── 两变量相关性
    ├── 线性关系？ → Pearson相关
    └── 非线性/有序？ → Spearman相关

# === 为什么需要校正？ ===
# 白话类比：
# 你扔1次硬币，出现正面的概率是50%
# 你扔20次硬币，至少有1次正面的概率是 1 - 0.5^20 ≈ 99.9999%

# 生信中的问题：
# 基因组有 20000 个基因，每个都做一次假设检验
# 即使所有基因都没差异，alpha=0.05 也会有 20000 × 0.05 = 1000 个假阳性！
# 这就叫"多重检验问题"

# === Bonferroni 校正（最严格）===
p_values <- c(0.001, 0.01, 0.03, 0.04, 0.06, 0.10)  # 6个p值
p_bonferroni <- p.adjust(p_values, method = "bonferroni")
# 原理：每个p值乘以检验次数（n=6）
# 0.001 × 6 = 0.006 < 0.05 → 显著
# 0.01 × 6 = 0.06 > 0.05 → 不显著了！
# 缺点：太严格了，很多真正的差异也会被"冤枉"

# === BH校正（Benjamini-Hochberg，生信最常用！）===
p_bh <- p.adjust(p_values, method = "BH")
# 原理：控制的是FDR（False Discovery Rate，假发现率）
# FDR = 在你报告的"显著"结果中，预期有多少比例是假阳性
# padj < 0.05 意味着：你报告的显著基因中，最多5%是假的

# === 对比 ===
comparison <- data.frame(
  p_value = p_values,
  Bonferroni = p_bonferroni,
  BH_FDR = p_bh
)
print(comparison)
# 可以看到：BH校正比Bonferroni宽松，保留更多真阳性

# === 白话版 FDR ===
# 场景：你做了差异分析，发现200个"显著"基因
# p_adj < 0.05（FDR < 5%）意味着：
# 这200个基因中，大约 200 × 0.05 = 10 个可能是假阳性
# 其他 190 个大概率是真正的差异基因

# === 面试常问：p-value 和 adjusted p-value 的区别？===
# p-value：单个检验的显著性（只考虑这一个基因）
# adjusted p-value（padj/FDR）：考虑了你同时做了几万次检验后的校正值
# 生信差异分析中，必须看 padj，不能只看 p-value！

# === q-value（另一种FDR衡量方式）===
# q-value ≈ BH adjusted p-value
# 生信论文中 q-value 和 padj 经常互换使用
# DESeq2 输出的 padj 就是 BH 校正后的 p 值

# === 简单线性回归 ===
# 白话：找一条直线，最好地描述x和y的关系
# 例：基因表达量 vs 药物浓度
drug_dose <- c(0, 0.1, 0.5, 1, 2, 5, 10)        # 药物浓度
gene_expr <- c(1.2, 1.5, 2.3, 3.1, 4.8, 8.2, 12.5)  # 基因表达量

model <- lm(gene_expr ~ drug_dose)   # lm = linear model
summary(model)
# 关键输出：
# Coefficients: 截距(Intercept)和斜率(drug_dose)
# R-squared: 决定系数，0-1之间，越接近1拟合越好
# p-value: 斜率是否显著不为0

# === 多元线性回归 ===
# 白话：多个因素同时影响结果
# 例：表达量 = β0 + β1×年龄 + β2×BMI + β3×性别
model_multi <- lm(expression ~ age + BMI + sex, data = clinical_df)
summary(model_multi)

# === 生信应用：批次效应校正 ===
# 用线性模型把批次效应从表达数据中"减掉"
library(limma)
# 设计矩阵：包含感兴趣的变量（condition）和要去掉的变量（batch）
design <- model.matrix(~ condition + batch, data = sample_info)
# removeBatchEffect 去除批次效应
corrected_expr <- removeBatchEffect(expr_matrix, batch = sample_info$batch)

# === 逻辑回归（用于二分类问题）===
# 白话：预测"是/否"类结果的概率
# 例：根据微生物组特征预测是否患有T2D

model_logistic <- glm(
  disease ~ Bacteroides + Prevotella + Faecalibacterium,
  data = microbiome_df,
  family = binomial    # binomial = 二分类（逻辑回归）
)
summary(model_logistic)
# 关键输出：
# Coefficients: 每个特征的系数（正=增加患病风险，负=降低）
# Odds Ratio = exp(coefficient)：优势比
# 例：exp(0.5) = 1.65，表示该菌每增加1单位，患病几率增加65%

# 预测
probabilities <- predict(model_logistic, newdata = test_df, type = "response")
predictions <- ifelse(probabilities > 0.5, "T2D", "Control")

# === 面试：线性回归 vs 逻辑回归 ===
# 线性回归：预测连续值（表达量、浓度）
# 逻辑回归：预测分类结果（患病/健康、响应/不响应）

library(survival)     # 生存分析核心包
library(survminer)    # 好看的生存曲线

# === KM生存分析 ===
# 白话：比较两组病人"谁活得更久"
# time: 随访时间（月）
# status: 是否发生事件（1=死亡, 0=存活/失访）
# group: 分组（如基因高表达 vs 低表达）

surv_obj <- Surv(time = clinical$time, event = clinical$status)  # 创建生存对象

# KM拟合
fit <- survfit(surv_obj ~ group, data = clinical)

# 画KM曲线
ggsurvplot(
  fit,
  data = clinical,
  pval = TRUE,                    # 显示log-rank检验p值
  conf.int = TRUE,                # 显示置信区间
  risk.table = TRUE,              # 显示风险表
  palette = c("#E41A1C", "#377EB8"),  # 自定义颜色
  xlab = "Time (months)",
  ylab = "Survival probability",
  title = "Overall Survival by Gene Expression"
)

# === 中位生存时间 ===
print(fit)  # 输出中包含 median survival time

# === Cox比例风险模型 ===
# 白话：分析哪些因素影响生存时间，每个因素的"风险比"是多少
cox_model <- coxph(
  Surv(time, status) ~ age + stage + gene_expr_high,
  data = clinical
)
summary(cox_model)
# 关键输出：
# HR (Hazard Ratio)：风险比
# HR > 1: 增加死亡风险（危险因素）
# HR < 1: 降低死亡风险（保护因素）
# HR = 1: 无影响
# 例：HR = 2.5 表示该因素使死亡风险增加150%

# 森林图（展示多个变量的HR）
ggforest(cox_model, data = clinical)

# === 贝叶斯定理（白话版）===
# P(A|B) = P(B|A) × P(A) / P(B)

# 白话类比：
# 你做了一个基因检测，结果显示"阳性"
# 这个检测的灵敏度（真阳性率）= 95%
# 特异度（真阴性率）= 90%
# 人群中这个病的发病率 = 1%

# 问：检测阳性，你真的有病的概率是多少？

sensitivity <- 0.95     # 灵敏度：有病检出阳性的概率
specificity <- 0.90     # 特异度：没病检出阴性的概率
prevalence <- 0.01      # 发病率

# 贝叶斯计算
p_positive <- sensitivity * prevalence +     # 真阳性
              (1 - specificity) * (1 - prevalence)  # 假阳性
ppv <- (sensitivity * prevalence) / p_positive  # 阳性预测值

cat(sprintf("检测阳性后真的有病的概率: %.1f%%\n", ppv * 100))
# 结果约：8.8%！（远低于直觉的95%！）
# 原因：发病率太低，大部分阳性都是假阳性

# === 生信中的贝叶斯应用 ===
# 1. DESeq2 内部：用贝叶斯收缩(shrinkage)估计log2FC
#    - 低表达/高方差的基因，FC被"拉"向0（更保守）
# 2. 系统发育分析：MrBayes 用贝叶斯方法推断进化树
# 3. 变异检测：GATK 用贝叶斯模型判断真突变 vs 测序错误
# 4. 宏基因组：Kraken2 用贝叶斯分类器对reads分类

# === 白话理解 ===
# 先验（Prior）：你在看到数据之前的"预判"
# 似然（Likelihood）：数据支持各个假设的程度
# 后验（Posterior）：综合先验和数据后的"更新判断"

# 后验 ∝ 先验 × 似然
# 白话：你原来的想法 + 新看到的证据 → 更新后的想法

# === 例子：DESeq2 的 shrinkage ===
# 先验：大部分基因不会有巨大的表达变化
# 数据：某个低表达基因看起来FC=10（但可能只是测了3个reads变成30个）
# 后验：考虑到低表达高方差，把FC收缩为更保守的估计

# 这就是为什么 DESeq2 的 lfcShrink() 比原始 log2FC 更可靠