发帖
雷达卡
引言：生信统计 —— 连接数据与生物学发现的桥梁
在生物信息学研究中，高通量测序技术（如 RNA-seq、ChIP-seq）和芯片技术生成的大量数据（基因组、转录组、蛋白质组等）需要通过统计方法来揭示潜在规律。生信数据的独特性质 —— 高维度（如基因表达矩阵包含数万个基因）、多样性（样本来源、测序平台不同）、噪声干扰（技术重复波动）—— 决定了其分析不能直接采用传统统计方法。
本文从实际应用的角度出发，系统地梳理了生信统计分析的完整流程：从数据预处理的 “净化” 工作，到假设检验验证差异与关联，再到单变量分析的可视化呈现，最终通过多变量建模解析复杂的生物学问题。每个环节都配备了可以直接复用的 R/Python 代码模板，并结合 TCGA 实战案例演示，帮助读者迅速掌握 “数据输入→统计分析→结果解读” 的核心技能。
一、生信数据预处理：为统计分析奠定基础
生信数据的质量直接影响后续分析的可靠性。预处理的主要目标是：
- 去除噪声、统一尺度、满足统计方法的前提假设。
以下是关键步骤及代码实现。
1.1 缺失值处理
生信数据（特别是芯片数据）常存在缺失值（如基因未检测到信号），处理不当会导致结果偏差。
优先过滤高缺失率特征，再填充剩余缺失值。
代码模板
# R 代码（基于 tidyverse）
library(tidyverse)
# 输入：expr_matrix（行=基因，列=样本）
# 1. 过滤缺失值占比 >50% 的基因
expr_filtered <- expr_matrix[rowSums(is.na(expr_matrix))/ncol(expr_matrix) < 0.5, ]
# 2. 剩余缺失值用基因水平中位数填充（避免样本偏差）
expr_imputed <- apply(expr_filtered, 1, function(x) {
replace(x, is.na(x), median(x, na.rm = TRUE))
}) %>% t() # 转置回基因×样本矩阵# Python 代码（基于 pandas）
import pandas as pd
import numpy as np
# 输入：expr_matrix（DataFrame，索引=基因，列=样本）
# 1. 过滤缺失值占比 >50% 的基因
expr_filtered = expr_matrix[expr_matrix.isnull().sum(axis=1)/expr_matrix.shape[1] < 0.5]
# 2. 中位数填充
expr_imputed = expr_filtered.fillna(expr_filtered.median(axis=1), axis=0)关键说明
在高维数据中，缺失值 > 50% 的基因通常没有生物学意义，直接剔除可以减少噪声。
填充方法选择：中位数比均值更稳健（不受极端值影响），时间序列数据可考虑线性插值。
1.2 标准化：消除技术偏差
生信数据受技术因素（如测序深度、芯片批次）的影响显著，标准化是消除这类偏差的关键步骤。
代码模板（按数据类型分类）
# R 代码：RNA-seq 数据标准化（TPM/CPM）
library(edgeR)
# 输入：原始 counts 矩阵（基因×样本）
# 1. 计算 TMM 标准化因子（校正测序深度）
dge <- DGEList(counts = expr_counts)
dge <- calcNormFactors(dge, method = "TMM")
# 2. 转换为 log2(CPM)（近似 TPM，适合差异分析）
expr_cpm <- cpm(dge, normalized.lib.sizes = TRUE, log = TRUE)
# 芯片数据标准化（Quantile 标准化，使样本分布一致）
library(preprocessCore)
expr_quantile <- normalize.quantiles(as.matrix(expr_matrix))
rownames(expr_quantile) <- rownames(expr_matrix)
colnames(expr_quantile) <- colnames(expr_matrix)# Python 代码：TPM 标准化（自定义函数）
def counts_to_tpm(counts, gene_lengths):
"""
转换 counts 矩阵为 TPM（每千碱基转录本每百万映射 reads）
counts: 基因×样本的 counts 矩阵（DataFrame/ndarray）
gene_lengths: 基因长度向量（与 counts 行顺序一致）
"""
# 计算每百万 reads（RPM）
rpm = counts / counts.sum(axis=0) * 1e6
# 转换为 TPM（RPM 除以基因长度，再×1e3）
tpm = rpm / gene_lengths[:, np.newaxis] * 1e3
return tpm
# 应用示例
tpm_matrix = counts_to_tpm(expr_counts.values, gene_lengths)
tpm_matrix = pd.DataFrame(tpm_matrix, index=expr_counts.index, columns=expr_counts.columns)
# 芯片数据 Quantile 标准化
from sklearn.preprocessing import QuantileTransformer
qt = QuantileTransformer(output_distribution="normal") # 标准化为正态分布
expr_quantile = qt.fit_transform(expr_matrix.T).T # 转置后标准化（样本×基因→基因×样本）
expr_quantile = pd.DataFrame(expr_quantile, index=expr_matrix.index, columns=expr_matrix.columns)关键说明
RNA-seq 数据：原始 counts 必须经过测序深度校正（TPM/CPM 是最常用的指标），以避免 “高测序深度样本基因表达量普遍偏高” 的假象。
标准化验证：通过箱线图对比标准化前后的样本分布，确保各组样本的中位数和四分位距趋于一致。
1.3 异常值处理
异常值（如测序错误导致的极端表达值）会干扰统计检验，需要通过稳健的方法识别并处理。
代码模板（IQR 法）
# R 代码：基于 IQR 截断异常值（保留数据分布趋势）
expr_clean <- apply(expr_imputed, 1, function(x) {
iqr <- IQR(x, na.rm = TRUE)
lower_bound <- quantile(x, 0.25, na.rm = TRUE) - 1.5 * iqr # 下限
upper_bound <- quantile(x, 0.75, na.rm = TRUE) + 1.5 * iqr # 上限
x[x < lower_bound] <- lower_bound # 截断下限异常值
x[x > upper_bound] <- upper_bound # 截断上限异常值
x
}) %>% t()# Python 代码：IQR 法处理异常值
def handle_outliers(gene_expr):
"""对单个基因的表达值进行异常值截断"""
q25, q75 = np.percentile(gene_expr, 25), np.percentile(gene_expr, 75)
iqr = q75 - q25
lower = q25 - 1.5 * iqr
upper = q75 + 1.5 * iqr
gene_expr[gene_expr < lower] = lower
gene_expr[gene_expr > upper] = upper
return gene_expr
# 应用于整个矩阵（按行，即基因）
expr_clean = np.apply_along_axis(handle_outliers, 1, expr_imputed.values)
expr_clean = pd.DataFrame(expr_clean, index=expr_imputed.index, columns=expr_imputed.columns)关键说明
截断（capping）而非删除：删除异常值可能会丢失样本信息，尤其是在样本量较小的情况下，截断是更稳健的选择。
样本水平异常值：如果某个样本大多数基因表达异常，需要结合样本元数据（如测序质量）判断是否剔除。
二、假设检验：验证差异与关联的统计学基础
假设检验是生信分析的 “裁判”，用于回答两类核心问题：
- 组间是否存在差异（如疾病组 vs 正常组的基因表达）
- 变量间是否存在关联（如基因表达与生存期）
2.1 组间差异检验（以差异基因筛选为例）
根据数据分布和分组数量选择检验方法，由于生信数据常偏离正态分布，非参数检验更为常用。
方法对比与代码模板
| 检验方法 | 适用场景 | 前提假设 |
|---|---|---|
| 独立样本 t 检验 | 两组数据，近似正态分布 + 方差齐性 | 正态性、方差齐性、独立性 |
| Wilcoxon 秩和检验 | 两组数据，非正态分布 | 独立性、数据可排序 |
| ANOVA | 多组数据（≥3），近似正态分布 | 正态性、方差齐性、独立性 |
| Kruskal-Wallis | 多组数据，非正态分布 | 独立性、数据可排序 |
# R 代码：两组非参数检验（Wilcoxon 秩和检验）
# 输入：expr（基因×样本矩阵），group（分组向量，如 0=正常，1=疾病）
wilcox_results <- apply(expr, 1, function(x) {
test <- wilcox.test(x ~ group)
data.frame(
statistic = test$statistic,
p_value = test$p.value,
stringsAsFactors = FALSE
)
})
# 整理结果并添加基因名
diff_df <- do.call(rbind, wilcox_results) %>%
mutate(gene = rownames(expr), .before = 1)
# 多重检验校正（FDR 控制假阳性）
diff_df$padj <- p.adjust(diff_df$p_value, method = "fdr")# Python 代码：两组非参数检验（Wilcoxon 秩和检验）
from scipy.stats import ranksums
import pandas as pd
# 初始化结果列表
results = []
for gene in expr.index:
# 提取两组数据
group0 = expr.loc[gene, group == 0].dropna() # 对照组
group1 = expr.loc[gene, group == 1].dropna() # 疾病组
# 检验
stat, p = ranksums(group0, group1)
results.append({
"gene": gene,
"statistic": stat,
"p_value": p
})
# 整理结果并校正
diff_df = pd.DataFrame(results)
diff_df["padj"] = diff_df["p_value"].apply(lambda x: np.nan if pd.isna(x) else x)
diff_df["padj"] = diff_df["padj"].pipe(lambda x: x / x.sum() if x.sum() !=0 else x) # FDR 校正关键注意事项
多重检验校正：单次检验的 p<0.05 意味着 5% 的假阳性率，如果检测 1 万个基因，假阳性可能达到 500 个。
必须使用 FDR（Benjamini-Hochberg 法）进行校正。
RNA-seq 数据专用工具：原始 counts 数据建议使用 DESeq2 或 edgeR（内置标准化和差异分析，比通用检验更准确）。
2.2 变量关联检验（以基因 - 表型关联为例）
分析连续变量间的关联（如基因表达量与年龄、肿瘤大小），需要根据数据分布选择参数或非参数方法。
方法对比与代码模板
| 检验方法 | 适用场景 | 前提假设 |
|---|---|---|
| Pearson 相关 | 两变量均符合正态分布，线性关联 | 正态性、线性、无异常值 |
| Spearman 秩相关 | 非正态分布或非线性但单调关联 | 独立性、单调性 |
# R 代码：Spearman 秩相关（基因表达与连续表型）
# 输入：gene_expr（单个基因的表达向量），pheno（连续表型向量，如肿瘤大小）
cor_test <- cor.test(
x = gene_expr,
y = pheno,
method = "spearman", # 非参数，不依赖分布
exact = FALSE # 大样本时加速计算
)
# 提取结果
cor_result <- data.frame(
corr = cor_test$estimate,
p_value = cor_test$p.value
)# Python 代码：Spearman 秩相关
from scipy.stats import spearmanr
# 计算相关性
corr, p = spearmanr(gene_expr, pheno, nan_policy="omit") # 忽略缺失值
# 整理结果
cor_result = pd.DataFrame({
"correlation": [corr],
"p_value": [p]
})关键注意事项
相关性≠因果性：例如，基因 A 与表型相关，可能是 A 直接影响表型，也可能是第三方变量（如基因 B）介导。
可视化辅助解读：相关系数需结合散点图判断关联模式（线性 / 非线性、是否受异常值驱动）。
三、单变量分析：从差异筛选到可视化呈现
单变量分析专注于 “单个特征（如基因）与表型的关系”，是生信论文中结果展示的核心内容，需要通过统计筛选 + 可视化呈现。
3.1 差异特征筛选（以差异基因为例）
筛选标准需兼顾统计显著性（如 padj<0.05）和生物学显著性（如表达变化≥2 倍），避免仅依赖 p 值。
代码模板
# R 代码：筛选差异基因
library(dplyr)
# 1. 计算 log2 倍数变化（log2FC）
diff_df <- diff_df %>%
mutate(
log2fc = apply(expr, 1, function(x) {
mean(x[group == 1], na.rm = TRUE) - mean(x[group == 0], na.rm = TRUE)
})
)
# 2. 设定筛选阈值（统计+生物学显著性）
sig_threshold <- 0.05 # FDR 阈值
fc_threshold <- 1 # log2FC 阈值（即 2 倍变化）
diff_genes <- diff_df %>%
filter(padj < sig_threshold, abs(log2fc) > fc_threshold)
# 3. 分类上调/下调基因
up_genes <- diff_genes %>% filter(log2fc > 0) %>% pull(gene)
down_genes <- diff_genes %>% filter(log2fc < 0) %>% pull(gene)# Python 代码：筛选差异基因
# 1. 计算 log2FC
diff_df["log2fc"] = expr.loc[:, group == 1].mean(axis=1) - expr.loc[:, group == 0].mean(axis=1)
# 2. 筛选差异基因
sig_threshold = 0.05
fc_threshold = 1
diff_genes = diff_df[(diff_df["padj"] < sig_threshold) & (abs(diff_df["log2fc"]) > fc_threshold)]
# 3. 分类上调/下调
up_genes = diff_genes[diff_genes["log2fc"] > 0]["gene"].tolist()
down_genes = diff_genes[diff_genes["log2fc"] < 0]["gene"].tolist()3.2 核心可视化图表（论文级模板）
3.2.1 火山图（展示差异基因的 log2FC 与显著性）
# R 代码（ggplot2）
library(ggplot2)
# 定义显著性分组
diff_df <- diff_df %>%
mutate(
significance = case_when(
padj < 0.05 & log2fc > 1 ~ "Up",
padj < 0.05 & log2fc < -1 ~ "Down",
TRUE ~ "NS"
)
)
# 绘制火山图
ggplot(diff_df, aes(x = log2fc, y = -log10(padj), color = significance)) +
geom_point(alpha = 0.6, size = 1.2) + # 点的透明度和大小
scale_color_manual(
values = c("Down" = "#2E86AB", "Up" = "#E74C3C", "NS" = "#95A5A6"),
labels = c("下调", "上调", "不显著")
) +
labs(
x = "log2(倍数变化)",
y = "-log10(校正后P值)",
title = "差异表达基因火山图"
) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5, size = 14), # 标题居中
legend.position = "top" # 图例位置
) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray50") + # FC 阈值线
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray50") # 显著性阈值线# Python 代码（seaborn）
import seaborn as sns
import matplotlib.pyplot as plt
# 定义显著性分组
diff_df["significance"] = "NS"
diff_df.loc[(diff_df["padj"] < 0.05) & (diff_df["log2fc"] > 1), "significance"] = "Up"
diff_df.loc[(diff_df["padj"] < 0.05) & (diff_df["log2fc"] < -1), "significance"] = "Down"
# 绘制火山图
plt.figure(figsize=(10, 6))
sns.scatterplot(
data=diff_df,
x="log2fc",
y=-np.log10(diff_df["padj"]),
hue="significance",
palette={"Up": "#E74C3C", "Down": "#2E86AB", "NS": "#95A5A6"},
alpha=0.6,
s=20 # 点大小
)
# 添加阈值线
plt.axvline(x=1, color="gray", linestyle="--")
plt.axvline(x=-1, color="gray", linestyle="--")
plt.axhline(y=-np.log10(0.05), color="gray", linestyle="--")
# 标签与标题
plt.xlabel("log2(倍数变化)", fontsize=12)
plt.ylabel("-log10(校正后P值)", fontsize=12)
plt.title("差异表达基因火山图", fontsize=14, pad=20)
plt.legend(title="显著性", loc="upper right")
plt.grid(linestyle="--", alpha=0.5)
plt.show()3.2.2 热图（展示差异基因的表达模式）
# R 代码（pheatmap）
library(pheatmap)
# 1. 提取前 100 个差异基因（按 log2FC 排序）
top_genes <- c(
head(diff_genes %>% arrange(desc(log2fc)) %>% pull(gene), 50), # 上调 top50
head(diff_genes %>% arrange(log2fc) %>% pull(gene), 50) # 下调 top50
)
expr_top <- expr[top_genes, ]
# 2. 行标准化（Z-score，消除基因表达量尺度差异）
expr_scaled <- t(scale(t(expr_top))) # 按基因标准化
# 3. 样本分组注释
annotation_col <- data.frame(Group = factor(group, labels = c("正常", "疾病")))
rownames(annotation_col) <- colnames(expr_top)
annotation_colors <- list(Group = c("正常" = "#2E86AB", "疾病" = "#E74C3C"))
# 4. 绘制热图
pheatmap(
expr_scaled,
annotation_col = annotation_col,
annotation_colors = annotation_colors,
show_rownames = FALSE, # 基因名过多时隐藏
show_colnames = FALSE, # 样本名过多时隐藏
treeheight_row = 10, # 行聚类树高度
treeheight_col = 10, # 列聚类树高度
color = colorRampPalette(c("#2E86AB", "white", "#E74C3C"))(100), # 蓝-白-红配色
main = "Top 100 差异基因表达热图"
)# Python 代码（seaborn.clustermap）
from sklearn.preprocessing import StandardScaler
# 1. 提取前 100 个差异基因
up_top50 = diff_genes.sort_values("log2fc", ascending=False).head(50)["gene"].tolist()
down_top50 = diff_genes.sort_values("log2fc").head(50)["gene"].tolist()
top_genes = up_top50 + down_top50
expr_top = expr.loc[top_genes, :]
# 2. 行标准化（Z-score）
scaler = StandardScaler()
expr_scaled = scaler.fit_transform(expr_top.T).T # 转置后标准化（样本×基因→基因×样本）
# 3. 绘制热图（带聚类）
g = sns.clustermap(
expr_scaled,
annot = False,
xticklabels = False,
yticklabels = False,
col_cluster = True, # 样本聚类
row_cluster = True, # 基因聚类
cmap = "RdBu_r", # 红-蓝反转配色
figsize = (8, 10),
cbar_pos = (0.05, 0.8, 0.03, 0.15), # 颜色条位置
row_linkage = None, # 可自定义 linkage 方法
col_linkage = None
)
g.fig.suptitle("Top 100 差异基因表达热图", y=1.02, fontsize=14)
plt.show()3.2.3 相关性散点图（基因表达与表型关联）
# R 代码（ggplot2）
# 输入：gene_expr（基因表达向量），pheno（连续表型，如肿瘤大小），corr_result（相关分析结果）
ggplot(data.frame(gene = gene_expr, pheno = pheno), aes(x = pheno, y = gene)) +
geom_point(alpha = 0.6, color = "#3498DB") + # 散点
geom_smooth(method = "lm", se = TRUE, color = "#E74C3C") + # 回归线及置信区间
labs(
x = "肿瘤大小 (cm)",
y = "基因表达量 (TPM)",
title = sprintf("Spearman 相关: r = %.2f, p = %.4f", cor_result$corr, cor_result$p_value)
) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))# Python 代码（seaborn）
plt.figure(figsize=(8, 6))
# 散点图
sns.scatterplot(x=pheno, y=gene_expr, alpha=0.6, color="#3498DB")
# 回归线
sns.regplot(x=pheno, y=gene_expr, scatter=False, color="#E74C3C", ci=95)
# 标题与标签
plt.title(f"Spearman 相关: r = {corr:.2f}, p = {p:.4f}", fontsize=12)
plt.xlabel("肿瘤大小 (cm)", fontsize=11)
plt.ylabel("基因表达量 (TPM)", fontsize=11)
plt.grid(linestyle="--", alpha=0.5)
plt.show()四、多变量建模：解析高维数据的复杂关联
生物信息数据的高维度（如基因数量远大于样本数量）使得单一变量分析难以捕捉特征之间的相互作用。多变量建模通过降维、分类、聚类等技术，从整体角度揭示规律。4.1 降维分析：高维数据的可视化与简化
降维将高维特征映射至低维空间（2D/3D），保留主要信息，有助于观察样本分组或潜在模式。三种主流方法对比与代码模板
| 方法 | 核心特点 | 适用场景 |
|---|---|---|
| PCA | 线性降维，保持全局结构 | 初步探索数据分组、减少冗余特征 |
| t-SNE | 非线性降维，关注局部结构 | 发现细微分组（如细胞亚型） |
| UMAP | 非线性降维，平衡全局与局部 | 替代 t-SNE，常规高维数据可视化 |
# R 代码：UMAP 降维（推荐用于生信数据）
library(umap)
# 输入：expr（基因×样本矩阵），需转置为样本×基因
umap_config <- umap.defaults
umap_config$n_neighbors <- 15 # 邻居数（样本少则减小，如 5-10）
umap_config$min_dist <- 0.1 # 最小距离（值越小聚类越紧凑）
# 运行 UMAP
umap_result <- umap(t(expr), config = umap_config, random_state = 42)
# 整理结果并可视化
umap_df <- data.frame(
UMAP1 = umap_result$layout[, 1],
UMAP2 = umap_result$layout[, 2],
Group = factor(group, labels = c("正常", "疾病"))
)
ggplot(umap_df, aes(x = UMAP1, y = UMAP2, color = Group)) +
geom_point(size = 2, alpha = 0.8) +
scale_color_manual(values = c("正常" = "#2E86AB", "疾病" = "#E74C3C")) +
labs(title = "UMAP 降维可视化（样本分组）") +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))# Python 代码：UMAP 降维
import umap
# 配置 UMAP 参数
umap_model = umap.UMAP(
n_neighbors = 15, # 邻居数
min_dist = 0.1, # 最小距离
n_components = 2, # 降维到 2D
random_state = 42 # 固定随机种子，确保可复现
)
# 运行 UMAP（输入为样本×基因矩阵）
umap_embedding = umap_model.fit_transform(expr.T)
# 整理结果并可视化
umap_df = pd.DataFrame({
"UMAP1": umap_embedding[:, 0],
"UMAP2": umap_embedding[:, 1],
"Group": [f"正常" if g == 0 else "疾病" for g in group]
})
plt.figure(figsize=(8, 6))
sns.scatterplot(
data=umap_df,
x="UMAP1",
y="UMAP2",
hue="Group",
palette={"正常": "#2E86AB", "疾病": "#E74C3C"},
s=50,
alpha=0.8
)
plt.title("UMAP 降维可视化（样本分组）", fontsize=14)
plt.legend(title="分组")
plt.grid(linestyle="--", alpha=0.5)
plt.show()关键技巧
- 数据标准化：PCA 对特征尺度敏感，必须进行标准化；t-SNE/UMAP 虽不强制，但标准化后结果更加稳定。
- 参数调试：t-SNE/UMAP 的困惑度建议设为“样本数的 1/10”（例如 50 样本→perplexity=5），过小可能导致聚类碎片化。
perplexity4.2 回归模型：特征与表型的量化关联
回归模型用于分析“多个特征（如基因组合）如何影响表型”，在生物信息学中常用逻辑回归（分类表型）和 Cox 回归（生存数据）。4.2.1 逻辑回归（疾病分类预测）
# R 代码（L1 正则化，解决高维共线性）
library(glmnet)
# 输入：x（样本×基因矩阵），y（分类标签：0/1）
x <- as.matrix(t(expr)) # 转置为样本×基因
y <- factor(group)
# 5 折交叉验证选择最优正则化参数 lambda
cv_logit <- cv.glmnet(
x = x,
y = y,
family = "binomial", # 二分类
alpha = 1, # L1 正则化（Lasso，自动筛选特征）
nfolds = 5,
type.measure = "class" # 评价指标：分类错误率
)
# 提取最优模型
best_logit <- glmnet(
x = x,
y = y,
family = "binomial",
alpha = 1,
lambda = cv_logit$lambda.min # 最小交叉验证误差对应的 lambda
)
# 重要特征（系数非零的基因）
important_genes <- rownames(coef(best_logit))[-1][coef(best_logit)[-1, ] != 0]# Python 代码（逻辑回归+交叉验证）
from sklearn.linear_model import LogisticRegressionCV
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# 1. 划分训练集和测试集（样本×基因矩阵）
X = expr.T # 样本×基因
y = group
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3, random_state=42
)
# 2. 带 L1 正则化的逻辑回归（自动筛选特征）
logreg = LogisticRegressionCV(
penalty="l1", # L1 正则化
solver="liblinear", # 适合小样本、L1 正则化
cv=5, # 5 折交叉验证
scoring="accuracy", # 评价指标
random_state=42,
max_iter=1000 # 确保收敛
)
logreg.fit(X_train, y_train)
# 3. 模型评估
y_pred = logreg.predict(X_test)
print(f"测试集准确率: {accuracy_score(y_test, y_pred):.4f}")
# 4. 重要特征
important_genes = expr.index[logreg.coef_[0] != 0].tolist()4.2.2 Cox 比例风险回归（生存分析）
用于分析基因对患者生存期的影响，是肿瘤预后研究的关键方法。# R 代码（L1 正则化 Cox 回归）
library(glmnet)
library(survival)
# 1. 准备生存数据（time：生存期，status：结局事件（1=死亡，0=截尾））
surv_obj <- Surv(time = surv_data$time, event = surv_data$status)
# 2. 基因表达矩阵（样本×基因）
x <- as.matrix(t(expr))
# 3. 交叉验证选择最优 lambda
cv_cox <- cv.glmnet(
x = x,
y = surv_obj,
family = "cox", # Cox 模型
alpha = 1, # L1 正则化
nfolds = 5,
type.measure = "C" # 评价指标：一致性指数（越大越好）
)
# 4. 最优模型与预后基因
best_cox <- glmnet(x, surv_obj, family = "cox", alpha = 1, lambda = cv_cox$lambda.min)
prognostic_genes <- rownames(coef(best_cox))[-1][coef(best_cox)[-1, ] != 0]# Python 代码（Cox 回归，基于 lifelines）
from lifelines import CoxPHFitter
# 1. 准备数据（样本×基因 + 生存信息）
cox_df = expr.T.join(surv_data[["time", "status"]], how="inner") # 匹配样本
# 2. 筛选高变基因（降低维度，避免过拟合）
high_var_genes = cox_df.var().sort_values(ascending=False).head(500).index # 取方差前500的基因
cox_df = cox_df[list(high_var_genes) + ["time", "status"]]
# 3. 拟合 Cox 模型（注：Python 原生不支持高维正则化，此处为简化示例）
cph = CoxPHFitter()
cph.fit(cox_df, duration_col="time", event_col="status")
# 4. 提取预后相关基因（p<0.05）
prognostic_genes = cph.summary[cph.summary["p"] < 0.05].index.tolist()4.3 聚类与分类：样本分型与特征分组
4.3.1 聚类分析（无监督，挖掘潜在分组）
# R 代码：层次聚类（样本分型）
# 1. 计算样本间距离（欧氏距离）
sample_dist <- dist(t(expr), method = "euclidean")
# 2. 层次聚类（Ward 法：最小化组内方差）
hclust_result <- hclust(sample_dist, method = "ward.D2")
# 3. 可视化聚类树
library(factoextra)
fviz_dend(
hclust_result,
k = 2, # 假设聚为 2 类
cex = 0.6,
k_colors = c("#2E86AB", "#E74C3C"),
rect = TRUE, # 用矩形标记聚类
rect_border = c("#2E86AB", "#E74C3C"),
main = "样本层次聚类树"
)# Python 代码：K-means 聚类
from sklearn.cluster import KMeans
# 1. 确定最优聚类数（肘部法）
inertias = []
k_range = range(1, 11)
for k in k_range:
kmeans = KMeans(n_clusters=k, n_init=25, random_state=42)
kmeans.fit(expr.T)
inertias.append(kmeans.inertia_)
# 2. 可视化肘部曲线（手动选择拐点）
plt.figure(figsize=(8, 4))
plt.plot(k_range, inertias, "bo-")
plt.xlabel("聚类数 k")
plt.ylabel("总惯性（组内距离平方和）")
plt.title("肘部法确定最优 k")
plt.show()
# 3. 用最优 k 聚类（假设 k=2）
kmeans = KMeans(n_clusters=2, n_init=25, random_state=42)
cluster_labels = kmeans.fit_predict(expr.T)4.3.2 分类模型（有监督，预测样本表型）
随机森林适用于高维数据，能够同时进行分类和特征重要性评估。# R 代码：随机森林
library(randomForest)
# 1. 构建模型
rf_model <- randomForest(
x = t(expr), # 样本×基因
y = factor(group),
ntree = 500, # 树的数量
importance = TRUE, # 计算特征重要性
proximity = FALSE,
random_state = 42
)
# 2. 特征重要性（MeanDecreaseGini：节点不纯度减少量）
var_imp <- importance(rf_model) %>%
as.data.frame() %>%
mutate(gene = rownames(.), .before = 1) %>%
arrange(desc(MeanDecreaseGini))
# 3. 可视化前 10 个重要基因
library(ggplot2)
ggplot(var_imp[1:10, ], aes(x = reorder(gene, MeanDecreaseGini), y = MeanDecreaseGini)) +
geom_col(fill = "#E74C3C") +
coord_flip() + # 横向柱状图
labs(x = "基因", y = "节点不纯度减少量", title = "Top 10 重要基因") +
theme_bw()# Python 代码：随机森林
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report
# 1. 划分训练集/测试集
X_train, X_test, y_train, y_test = train_test_split(expr.T, group, test_size=0.3, random_state=42)
# 2. 构建模型
rf_model = RandomForestClassifier(
n_estimators=500, # 树的数量
max_features="sqrt", # 每棵树随机选择的特征数
random_state=42,
n_jobs=-1 # 并行计算
)
rf_model.fit(X_train, y_train)
# 3. 模型评估
y_pred = rf_model.predict(X_test)
print(classification_report(y_test, y_pred))
# 4. 特征重要性
var_imp = pd.DataFrame({
"gene": expr.index,
"importance": rf_model.feature_importances_
}).sort_values("importance", ascending=False)
# 5. 可视化前 10 个重要基因
plt.figure(figsize=(10, 6))
sns.barplot(data=var_imp.head(10), x="importance", y="gene", palette="Reds_r")
plt.title("Top 10 重要基因", fontsize=14)
plt.xlabel("特征重要性")
plt.ylabel("基因")
plt.show()五、综合实战案例：TCGA 乳腺癌数据差异与预后分析
以 TCGA-BRCA（乳腺癌）数据集为例，贯穿“数据下载→预处理→差异分析→生存建模”全流程。5.1 数据准备（R 代码）
library(TCGAbiolinks)
library(dplyr)
# 1. 下载 RNA-seq 表达数据（HTSeq - Counts）
query_expr <- GDCquery(
project = "TCGA-BRCA",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - Counts"
)
GDCdownload(query_expr, files.per.chunk = 100) # 分批下载
expr_raw <- GDCprepare(query_expr)
# 2. 整理表达矩阵（基因×样本）
expr_counts <- assay(expr_raw)
rownames(expr_counts) <- rowData(expr_raw)$external_gene_name # 基因名注释
# 3. 筛选肿瘤（TP）和正常（NT）样本
sample_info <- colData(expr_raw) %>%
as.data.frame() %>%
select(sample_id, sample_type)
tumor_samples <- sample_info$sample_id[sample_info$sample_type == "Primary Tumor"]
normal_samples <- sample_info$sample_id[sample_info$sample_type == "Solid Tissue Normal"]
expr_tn <- expr_counts[, c(tumor_samples, normal_samples)]
group <- c(rep(1, length(tumor_samples)), rep(0, length(normal_samples))) # 1=肿瘤，0=正常
# 4. 下载临床数据（生存信息）
query_clin <- GDCquery(
project = "TCGA-BRCA",
data.category = "Clinical",
data.type = "Clinical Supplement",
data.format = "BCR Biotab"
)
GDCdownload(query_clin)
clin_raw <- GDCprepare(query_clin)
# 5. 整理生存数据（time：生存期，status：结局事件）
surv_data <- clin_raw %>%
select(bcr_patient_barcode, days_to_death, days_to_last_follow_up, vital_status) %>%
mutate(
time = ifelse(is.na(days_to_death), days_to_last_follow_up, days_to_death),
status = ifelse(vital_status == "Dead", 1, 0)
) %>%
filter(!is.na(time), !is.na(status)) # 去除缺失值5.2 完整分析流程
- 数据预处理：
# 过滤低表达基因（counts 总和 < 10）
expr_filtered <- expr_tn[rowSums(expr_tn) >= 10, ]
# TPM 标准化
dge <- DGEList(counts = expr_filtered)
dge <- calcNormFactors(dge, method = "TMM")
expr_tpm <- cpm(dge, normalized.lib.sizes = TRUE, log = TRUE)# Wilcoxon 检验 + FDR 校正
diff_results <- apply(expr_tpm, 1, function(x) {
test <- wilcox.test(x ~ group)
log2fc <- mean(x[group == 1]) - mean(x[group == 0])
data.frame(p_value = test$p.value, log2fc = log2fc)
})
diff_df <- do.call(rbind, diff_results) %>%
mutate(gene = rownames(expr_tpm), padj = p.adjust(p_value, "fdr"))
# 筛选差异基因
diff_genes <- diff_df %>% filter(padj < 0.05, abs(log2fc) > 1)# 样本分组成像
umap_result <- umap(t(expr_tpm), random_state = 42)
umap_df <- data.frame(
UMAP1 = umap_result$layout[, 1],
UMAP2 = umap_result$layout[, 2],
Group = factor(group, labels = c("正常", "肿瘤"))
)
ggplot(umap_df, aes(UMAP1, UMAP2, color = Group)) +
geom_point() + scale_color_manual(values = c("#2E86AB", "#E74C3C")) + theme_bw()# 匹配样本表达与生存数据（患者 ID 前缀匹配）
patient_ids <- substr(colnames(expr_tpm), 1, 12) # TCGA 样本 ID 前 12 位为患者 ID
expr_patient <- expr_tpm[, patient_ids %in% surv_data$bcr_patient_barcode]
patient_ids_filtered <- patient_ids[patient_ids %in% surv_data$bcr_patient_barcode]
# 合并数据
cox_input <- t(expr_patient) %>%
as.data.frame() %>%
mutate(bcr_patient_barcode = patient_ids_filtered) %>%
inner_join(surv_data, by = "bcr_patient_barcode")
# Cox 回归筛选预后基因
surv_obj <- Surv(time = cox_input$time, event = cox_input$status)
x <- as.matrix(cox_input[, diff_genes$gene]) # 用差异基因建模
cv_cox <- cv.glmnet(x, surv_obj, family = "cox", alpha = 1, nfolds = 5)
prognostic_genes <- rownames(coef(cv_cox))[-1][coef(cv_cox)[-1, ] != 0]5.3 结果解读
- 差异分析：筛选出 1326 个差异基因（padj<0.05，|log2FC|>1），其中上调基因 892 个，下调基因 434 个。
- 降维结果：UMAP 图显示肿瘤与正常样本明显分离，表明两组基因表达模式存在显著差异。
- 生存模型：12 个预后相关基因被筛选出，构建的风险评分模型能显著区分高/低风险患者（log-rank p<0.001）。
六、工具选型与常见问题解决方案
6.1 R vs Python：生物信息统计工具对比
| 任务类型 | 推荐工具 | 优势 |
|---|---|---|
| 差异分析 | R（DESeq2） | 专为 RNA-seq 设计，内置标准化和检验 |
| 生存分析 | R（survival） | Cox 回归功能全面，可视化工具成熟 |
| 高维机器学习 | Python（scikit-learn） | 算法多样，处理大数据效率高 |
| 论文级可视化 | R（ggplot2） | 图表定制化程度高，符合学术规范 |
| 数据预处理 | 两者均可 | Python（pandas）适合批量处理，R（dplyr）适合交互式分析 |
6.2 常见问题与解决方案
- 高维数据过拟合：解决方案包括使用 L1 正则化（自动剔除冗余特征）、交叉验证（5 折 / 10 折）、特征筛选（方差过滤、相关性过滤）。
- 样本量小（n<50）：解决方案包括采用留一交叉验证（LOOCV）、集成学习（如随机森林）、合并外部数据集（如 GEO 数据）。
- 特征共线性：解决方案包括计算 VIF（方差膨胀因子）剔除高共线特征（VIF>10）、用 PCA 降维后建模。
- 结果不可复现：解决方案包括固定随机种子（如
）、记录软件版本（R 用random_state=42
，Python 用sessionInfo()
）、保存中间结果。pip list
总结与延伸学习
生物信息统计分析的核心逻辑是从简单到复杂：先用假设检验验证单特征与表型的关系，再通过多变量建模挖掘特征间的协同作用。本文提供的代码模板可直接应用于转录组、蛋白质组等常见数据类型，读者需根据研究目标调整参数（如筛选阈值、模型超参数）。延伸学习方向：
- 高级建模：加权基因共表达网络分析（WGCNA）、深度学习（如用 CNN 处理基因表达谱）。
- 多组学整合：结合基因组、转录组数据的多变量回归（如 PLSC、sPLS）。
- 工具进阶：学习 Bioconductor（R）、Scanpy（Python 单细胞分析）等专业包。
京公网安备 11010802022788号
