GSE19697：术后复发 vs 无复发乳腺癌的基因表达差异分析
新辅助化疗前肿瘤活检 · Affymetrix HG-U133 Plus 2.0 (GPL570) · limma
2026-10-16
1 研究背景与数据
临床问题：新辅助化疗可以缩小肿瘤、利于手术完全切除，并清除微小转移以降低复发风险。但不同患者对化疗的敏感性差异很大，临床上需要能够预测疗效/复发的分子标志物。
本数据 GSE19697（Lin et al., Clin Cancer Res 2009, PMID 19967557）是华盛顿大学一项新辅助化疗试验的基线（化疗前）肿瘤活检表达谱：
| 项目 | 内容 |
|---|---|
| 平台 | Affymetrix HG-U133 Plus 2.0（GPL570，54 675 个探针集） |
| 样本 | 24 例化疗前 core biopsy（来自 70 例中的 basal-like 亚组） |
| 分组（本分析） | recur = 1 术后复发 vs
recur = 0 无复发 |
| 其他临床变量 | ER / PR / HER2 状态、组织学分级、阳性淋巴结数、T 分期、种族等 |
| 数据来源 | GEO: GSE19697 |
原始数据有两个层次，本流程两者都会用到：
GSE19697_series_matrix.txt.gz：作者用 MAS 5.0 处理好的表达值（线性尺度、scaled to target intensity 1500），同时携带全部样本临床注释。GSE19697_RAW.tar：24 个原始.CEL文件，用于从原始荧光强度出发做标准化的 RMA 处理，是更规范、可复现的做法。
分析策略：以 CEL → RMA 为主流程；若原始 CEL 不存在，则退化为 series matrix 的 MAS 5.0 值 +
log2+ 分位数标准化。两条路径的结果会在最后一节做一致性比较。
2 环境准备
library(limma) # 差异表达分析（线性模型 + 经验贝叶斯）
library(affy) # 读取 CEL、RMA 标准化
library(hgu133plus2.db) # HG-U133 Plus 2.0 探针注释
library(AnnotationDbi)
library(Biobase)
library(ggplot2)
library(ggrepel)
library(pheatmap)
library(RColorBrewer)
library(matrixStats)
set.seed(2026)3 数据下载与解压
data_dir <- "data" # 数据目录（与 Rmd 同级）
cel_dir <- file.path(data_dir, "CEL") # CEL 文件目录
res_dir <- "results" # 结果输出目录
dir.create(cel_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(res_dir, recursive = TRUE, showWarnings = FALSE)
sm_file <- file.path(data_dir, "GSE19697_series_matrix.txt.gz")
raw_file <- file.path(data_dir, "GSE19697_RAW.tar")
sm_url <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE19nnn/GSE19697/matrix/GSE19697_series_matrix.txt.gz"
raw_url <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE19nnn/GSE19697/suppl/GSE19697_RAW.tar"
# 已下载过就跳过，避免重复下载（series matrix 约 3 MB，RAW tar 约 112 MB）
download_if_missing <- function(url, dest) {
if (!file.exists(dest) || file.size(dest) == 0) {
message("正在下载 ", basename(dest), " ...")
utils::download.file(url, dest, mode = "wb", quiet = TRUE)
} else {
message("已存在，跳过下载：", basename(dest))
}
}
download_if_missing(sm_url, sm_file)
download_if_missing(raw_url, raw_file)
# base R 解压 gzip，避免引入额外依赖
gunzip_file <- function(f) {
out <- sub("\\.gz$", "", f)
con_in <- gzfile(f, "rb")
con_out <- file(out, "wb")
repeat {
chunk <- readBin(con_in, "raw", n = 1e6)
if (length(chunk) == 0) break
writeBin(chunk, con_out)
}
close(con_in); close(con_out)
unlink(f)
out
}
# RAW tar -> CEL/*.CEL
cel_files <- list.files(cel_dir, pattern = "\\.CEL$", full.names = TRUE)
if (length(cel_files) == 0) {
utils::untar(raw_file, exdir = cel_dir)
for (f in list.files(cel_dir, pattern = "\\.CEL\\.gz$", full.names = TRUE)) gunzip_file(f)
cel_files <- list.files(cel_dir, pattern = "\\.CEL$", full.names = TRUE)
}
cel_files <- cel_files[order(basename(cel_files))] # 按 GSM 编号排序，与临床表一致
knitr::kable(
data.frame(
文件 = c(basename(sm_file), basename(raw_file), paste0("CEL 文件（", length(cel_files), " 个）")),
大小 = c(paste(round(file.size(sm_file) / 1024^2, 1), "MB"),
paste(round(file.size(raw_file) / 1024^2, 1), "MB"),
paste(round(sum(file.size(cel_files)) / 1024^2, 1), "MB")),
路径 = c(data_dir, data_dir, cel_dir)
),
caption = "已就绪的数据文件"
)| 文件 | 大小 | 路径 |
|---|---|---|
| GSE19697_series_matrix.txt.gz | 3.2 MB | data |
| GSE19697_RAW.tar | 106.8 MB | data |
| CEL 文件（24 个） | 310.2 MB | data/CEL |
4 临床注释的提取与格式化
series matrix 的注释信息全部写在文件头部的 !Sample_*
行里（每行一个变量，制表符分隔各样本）。下面用 base R 直接解析，便于看清
SOFT / GEOSoft 文件的结构。
sm_lines <- readLines(gzfile(sm_file), warn = FALSE)
# 取出某个 "!Sample_xxx" 行里所有样本的取值
parse_row <- function(lines, tag) {
hit <- grep(paste0("^!", tag, "\t"), lines, value = TRUE)
if (length(hit) == 0) return(NULL)
gsub('^"|"$', "", strsplit(hit[1], "\t")[[1]][-1])
}
# 每个 !Sample_characteristics_ch1 行是一种临床变量（key: value 形式），本数据共 12 行
char_lines <- grep("^!Sample_characteristics_ch1\t", sm_lines, value = TRUE)
char_list <- lapply(char_lines, function(l) {
v <- gsub('^"|"$', "", strsplit(l, "\t")[[1]][-1])
sub("^[^:]*:\\s*", "", v) # 去掉 "key: " 前缀，只留取值
})
names(char_list) <- vapply(char_lines, function(l) {
sub(":.*$", "", gsub('^"', "", strsplit(l, "\t")[[1]][2])) # 变量名取第一条记录的 key
}, character(1))
pheno <- data.frame(gsm = parse_row(sm_lines, "Sample_geo_accession"),
title = parse_row(sm_lines, "Sample_title"),
as.data.frame(char_list, check.names = FALSE, stringsAsFactors = FALSE),
check.names = FALSE)
# 变量重命名 + 类型转换（recur：1 = 术后复发，0 = 无复发）
pheno <- within(pheno, {
recurrence <- factor(ifelse(recur == "1", "Recurrence", "NoRecurrence"),
levels = c("NoRecurrence", "Recurrence"))
# 显式给出 levels，某些变量在本队列中只有一个水平（如 PR 全为阴性）
er <- factor(er, levels = c("0", "1"), labels = c("negative", "positive"))
pr <- factor(pr, levels = c("0", "1"), labels = c("negative", "positive"))
her2 <- factor(her2, levels = c("0", "1"), labels = c("negative", "positive"))
grade <- factor(grade)
t_stage <- factor(`t stage`)
pos_nodes <- as.numeric(`positive nodes`)
nodes_exam <- as.numeric(`nodes examined`)
race <- factor(race)
})
cat("样本数：", nrow(pheno),
" 无复发：", sum(pheno$recurrence == "NoRecurrence"),
" 复发：", sum(pheno$recurrence == "Recurrence"), "\n")样本数： 24 无复发： 18 复发： 6 knitr::kable(
data.frame(GSM = pheno$gsm,
样本 = sub("Pre-Treatment breast tumor ", "tumor ", pheno$title),
复发 = as.character(pheno$recurrence),
ER = pheno$er, PR = pheno$pr, HER2 = pheno$her2,
分级 = pheno$grade, 阳性淋巴结 = pheno$pos_nodes, T分期 = pheno$t_stage),
caption = "24 例样本的临床注释（解析自 series matrix）"
)| GSM | 样本 | 复发 | ER | PR | HER2 | 分级 | 阳性淋巴结 | T分期 |
|---|---|---|---|---|---|---|---|---|
| GSM491610 | tumor 7411 | NoRecurrence | negative | negative | positive | 3 | 0 | 2 |
| GSM491611 | tumor 7412 | Recurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491612 | tumor 7415 | Recurrence | negative | negative | negative | 3 | 3 | 2 |
| GSM491613 | tumor 7421 | NoRecurrence | negative | negative | negative | 2 | 0 | 2 |
| GSM491614 | tumor 7422 | NoRecurrence | negative | negative | negative | 3 | 0 | 3 |
| GSM491615 | tumor 8159 | NoRecurrence | negative | negative | negative | 3 | 5 | 2 |
| GSM491616 | tumor 8170 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491617 | tumor 8189 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491618 | tumor 8190 | NoRecurrence | negative | negative | negative | 3 | 1 | 2 |
| GSM491619 | tumor 8191 | Recurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491620 | tumor 8199 | NoRecurrence | negative | negative | negative | 3 | 0 | 1 |
| GSM491621 | tumor 8200 | NoRecurrence | negative | negative | negative | 3 | 1 | 2 |
| GSM491622 | tumor 8209 | Recurrence | negative | negative | positive | 3 | 1 | 2 |
| GSM491623 | tumor 8210 | NoRecurrence | positive | negative | negative | 3 | 0 | 2 |
| GSM491624 | tumor 8413 | Recurrence | negative | negative | negative | 2 | 10 | 3 |
| GSM491625 | tumor 11519 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491626 | tumor 11524 | Recurrence | negative | negative | negative | 3 | 14 | 3 |
| GSM491627 | tumor 11526 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491628 | tumor 11529 | NoRecurrence | negative | negative | negative | 3 | 1 | 2 |
| GSM491629 | tumor 11530 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491630 | tumor 11532 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491631 | tumor 11571 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491632 | tumor 11582 | NoRecurrence | negative | negative | negative | 3 | 0 | 2 |
| GSM491633 | tumor 11585 | NoRecurrence | negative | negative | negative | 3 | 2 | 2 |
# 两组临床特征的分布（小样本，仅作描述性统计）
clin_vars <- c("er", "pr", "her2", "grade", "t_stage", "race")
tab <- do.call(rbind, lapply(clin_vars, function(v) {
t <- table(pheno[[v]], pheno$recurrence)
data.frame(变量 = v, 水平 = rownames(t),
无复发 = t[, "NoRecurrence"], 复发 = t[, "Recurrence"], row.names = NULL)
}))
knitr::kable(tab, caption = "两组临床特征分布")| 变量 | 水平 | 无复发 | 复发 |
|---|---|---|---|
| er | negative | 17 | 6 |
| er | positive | 1 | 0 |
| pr | negative | 18 | 6 |
| pr | positive | 0 | 0 |
| her2 | negative | 17 | 5 |
| her2 | positive | 1 | 1 |
| grade | 2 | 1 | 1 |
| grade | 3 | 17 | 5 |
| t_stage | 1 | 1 | 0 |
| t_stage | 2 | 16 | 4 |
| t_stage | 3 | 1 | 2 |
| race | African American | 8 | 2 |
| race | Caucasian | 10 | 4 |
cat("\n阳性淋巴结数（均值 ± 标准差）：无复发 ",
sprintf("%.2f ± %.2f", mean(pheno$pos_nodes[pheno$recurrence == "NoRecurrence"]),
sd(pheno$pos_nodes[pheno$recurrence == "NoRecurrence"])),
"，复发 ",
sprintf("%.2f ± %.2f", mean(pheno$pos_nodes[pheno$recurrence == "Recurrence"]),
sd(pheno$pos_nodes[pheno$recurrence == "Recurrence"])), "\n")
阳性淋巴结数（均值 ± 标准差）：无复发 0.56 ± 1.25 ，复发 4.67 ± 5.92 5 表达数据与标准化
5.1 读取 series matrix（备用路径 + 稳健性比较）
series matrix 的表体以 ID_REF 行开头，最后一行是
!series_matrix_table_end
标记，读入后必须把它剔除，否则会多出一行全为缺失值的”探针”。本文件表体含
54 674 个探针集，比芯片设计的 54 675 个少 1
个（作者提交数据时的差异，不影响分析）。
skip_n <- grep("^\"?ID_REF", sm_lines)[1]
expr_mas5 <- as.matrix(read.delim(gzfile(sm_file), skip = skip_n, header = TRUE,
row.names = 1, check.names = FALSE))
# 剔除表尾的 "!series_matrix_table_end" 行（不含任何表达值）
expr_mas5 <- expr_mas5[!startsWith(rownames(expr_mas5), "!"), , drop = FALSE]
storage.mode(expr_mas5) <- "double"
stopifnot(!anyNA(expr_mas5))
cat(sprintf("MAS 5.0 表达矩阵：%d 个探针集 × %d 个样本\n", nrow(expr_mas5), ncol(expr_mas5)))MAS 5.0 表达矩阵：54674 个探针集 × 24 个样本cat(sprintf("取值范围：%.1f ~ %.1f（线性尺度，尚未做 log 转换）\n", min(expr_mas5), max(expr_mas5)))取值范围：0.6 ~ 531159.8（线性尺度，尚未做 log 转换）5.2 CEL → RMA（主流程）
RMA 一个函数就完成了四步：背景校正 → log2 转换 → 分位数标准化 → 探针集汇总（median polish），因此 RMA 结果无需再做额外标准化。
use_cel <- TRUE # TRUE：用 RAW tar 里的 CEL 走 RMA；FALSE：退化为 MAS5 + 分位数标准化
cel_ok <- use_cel && length(cel_files) == nrow(pheno) &&
requireNamespace("hgu133plus2cdf", quietly = TRUE)
if (cel_ok) {
abatch <- ReadAffy(filenames = cel_files) # 原始 CEL：未做任何标准化
raw_pm <- log2(pmax(exprs(abatch), 1)) # 标准化前的 PM 强度（log2 仅为可视化）
eset <- rma(abatch) # RMA 标准化
expr_norm <- exprs(eset)
} else {
message("未使用 CEL：退化为 series matrix 的 MAS 5.0 值 + log2 + 分位数标准化")
raw_pm <- log2(pmax(expr_mas5, 1))
expr_norm <- normalizeBetweenArrays(raw_pm, method = "quantile")
}Background correcting
Normalizing
Calculating Expression# ReadAffy 保留 ".CEL" 扩展名，统一成 GSM 编号后再与临床表对应
colnames(raw_pm) <- sub("[.]CEL$", "", basename(colnames(raw_pm)))
colnames(expr_norm) <- sub("[.]CEL$", "", basename(colnames(expr_norm)))
# 让表达矩阵的列顺序与临床表严格对应
pheno <- pheno[match(colnames(expr_norm), pheno$gsm), , drop = FALSE]
stopifnot(!anyNA(pheno$gsm), identical(colnames(expr_norm), pheno$gsm))
cat(sprintf("标准化后表达矩阵：%d 个探针集 × %d 个样本（log2 尺度）\n",
nrow(expr_norm), ncol(expr_norm)))标准化后表达矩阵：54675 个探针集 × 24 个样本（log2 尺度）5.3 标准化效果可视化
par(mfrow = c(1, 2), mar = c(7, 4, 3, 1))
boxplot(raw_pm, las = 2, outline = FALSE, cex.axis = 0.8,
main = "标准化前（原始 PM 强度，log2）", ylab = expression(log[2]~"intensity"))
boxplot(expr_norm, las = 2, outline = FALSE, cex.axis = 0.8,
main = if (cel_ok) "RMA 标准化后" else "分位数标准化后",
ylab = expression(log[2]~"expression"))
par(mfrow = c(1, 1))# 密度图：标准化后各样本的分布应基本重合
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plotDensities(raw_pm, legend = FALSE, main = "标准化前")
plotDensities(expr_norm, legend = FALSE,
main = if (cel_ok) "RMA 标准化后" else "分位数标准化后")
par(mfrow = c(1, 1))# 用前 5000 个高变探针集做 PCA，观察样本整体分布是否与分组有关
probe_var <- rowVars(expr_norm)
top_var <- order(probe_var, decreasing = TRUE)[seq_len(5000)]
pca <- prcomp(t(expr_norm[top_var, ]), scale. = TRUE)
pv <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)
pca_df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2],
recurrence = pheno$recurrence)
ggplot(pca_df, aes(PC1, PC2, color = recurrence)) +
geom_point(size = 3.5, alpha = 0.9) +
stat_ellipse(level = 0.68, linetype = 2) +
labs(title = "PCA（前 5000 个高变探针集）",
x = paste0("PC1 (", pv[1], "%)"), y = paste0("PC2 (", pv[2], "%)"))
6 探针注释与基因水平汇总
anno <- AnnotationDbi::select(hgu133plus2.db, keys = rownames(expr_norm),
keytype = "PROBEID",
columns = c("SYMBOL", "ENTREZID", "GENENAME"))
anno <- anno[!is.na(anno$SYMBOL), ]
anno <- anno[!duplicated(anno$PROBEID), ]
cat(sprintf("探针集 %d 个，其中 %d 个可注释到基因符号，涉及 %d 个基因\n",
nrow(expr_norm), nrow(anno), length(unique(anno$SYMBOL))))探针集 54675 个，其中 44647 个可注释到基因符号，涉及 21355 个基因# 一个基因对应多个探针时，保留平均表达最高的探针（基因水平汇总常用做法）
anno$probe_mean <- rowMeans(expr_norm[anno$PROBEID, , drop = FALSE])
anno_sorted <- anno[order(anno$SYMBOL, -anno$probe_mean), ]
keep_probe <- anno_sorted$PROBEID[!duplicated(anno_sorted$SYMBOL)]
expr_gene <- expr_norm[keep_probe, , drop = FALSE]
rownames(expr_gene) <- anno_sorted$SYMBOL[match(keep_probe, anno_sorted$PROBEID)]
expr_gene <- expr_gene[order(rownames(expr_gene)), , drop = FALSE]
cat(sprintf("基因水平表达矩阵：%d 个基因 × %d 个样本\n", nrow(expr_gene), ncol(expr_gene)))基因水平表达矩阵：21355 个基因 × 24 个样本7 差异表达分析（limma）
group <- pheno$recurrence
design <- model.matrix(~ group) # 截距 + groupRecurrence（复发相对于无复发的 log2 差异）
colnames(design) <- c("Intercept", "Recurrence_vs_NoRecurrence")
fit <- lmFit(expr_norm, design) # 探针集水平拟合
fit <- eBayes(fit) # 经验贝叶斯收缩，提高小样本下的稳定性
tt <- topTable(fit, coef = "Recurrence_vs_NoRecurrence", number = Inf, sort.by = "P")
tt$probe_id <- rownames(tt)
tt$symbol <- anno$SYMBOL[match(tt$probe_id, anno$PROBEID)]
tt <- tt[, c("probe_id", "symbol", "logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "B")]
# 显著性判定：FDR < 0.05 且 |log2FC| > 1
sig <- decideTests(fit, p.value = 0.05, lfc = 1, adjust.method = "BH")
cat("差异表达探针集统计（FDR < 0.05, |log2FC| > 1）：\n")差异表达探针集统计（FDR < 0.05, |log2FC| > 1）：print(summary(sig)) Intercept Recurrence_vs_NoRecurrence
Down 0 0
NotSig 0 54675
Up 54675 0de_fdr <- subset(tt, adj.P.Val < 0.05 & abs(logFC) > 1)
cat(sprintf("\n严格标准共筛出 %d 个探针集，最小 FDR = %.3f\n", nrow(de_fdr), min(tt$adj.P.Val)))
严格标准共筛出 0 个探针集，最小 FDR = 1.000cat("复发组仅 6 例，经 54 675 次检验的 BH 校正后没有探针集达到 FDR < 0.05，\n")复发组仅 6 例，经 54 675 次检验的 BH 校正后没有探针集达到 FDR < 0.05，cat("这是小样本芯片研究的常见结果，因此下面另用一个探索性标准（名义 P 值）筛选候选基因。\n\n")这是小样本芯片研究的常见结果，因此下面另用一个探索性标准（名义 P 值）筛选候选基因。# 探索性判定：名义 P < 0.05 且 |log2FC| > 1 —— 仅用于候选基因筛选与可视化，必须外部验证
cand <- subset(tt, P.Value < 0.05 & abs(logFC) > 1)
cand <- cand[order(cand$P.Value), ]
cat(sprintf("探索性标准（P < 0.05 且 |logFC| > 1）：%d 个探针集，上调 %d 个，下调 %d 个\n",
nrow(cand), sum(cand$logFC > 0), sum(cand$logFC < 0)))探索性标准（P < 0.05 且 |logFC| > 1）：39 个探针集，上调 27 个，下调 12 个knitr::kable(head(tt[!is.na(tt$symbol), ], 20), digits = 3,
caption = "按 P 值排序的前 20 个探针集（Recurrence vs NoRecurrence）")| probe_id | symbol | logFC | AveExpr | t | P.Value | adj.P.Val | B | |
|---|---|---|---|---|---|---|---|---|
| 210501_x_at | 210501_x_at | EIF3K | -0.811 | 10.287 | -4.594 | 0.000 | 1 | -2.517 |
| 201441_at | 201441_at | COX6B1 | -0.725 | 10.366 | -4.525 | 0.000 | 1 | -2.558 |
| 225969_at | 225969_at | ALKBH6 | -0.854 | 6.242 | -4.485 | 0.000 | 1 | -2.581 |
| 221494_x_at | 221494_x_at | EIF3K | -0.775 | 10.339 | -4.301 | 0.000 | 1 | -2.693 |
| 203262_s_at | 203262_s_at | FAM50A | -0.675 | 8.013 | -4.109 | 0.000 | 1 | -2.811 |
| 227138_at | 227138_at | CRTAP | 0.566 | 5.773 | 4.043 | 0.000 | 1 | -2.853 |
| 202053_s_at | 202053_s_at | ALDH3A2 | -1.122 | 8.133 | -4.032 | 0.000 | 1 | -2.859 |
| 212716_s_at | 212716_s_at | EIF3K | -0.707 | 10.654 | -4.002 | 0.001 | 1 | -2.878 |
| 224913_s_at | 224913_s_at | TIMM50 | -0.748 | 8.070 | -3.901 | 0.001 | 1 | -2.943 |
| 202475_at | 202475_at | TMEM147 | -0.732 | 9.806 | -3.864 | 0.001 | 1 | -2.966 |
| 212955_s_at | 212955_s_at | POLR2I | -0.718 | 8.915 | -3.856 | 0.001 | 1 | -2.971 |
| 226437_at | 226437_at | YIF1B | -0.814 | 7.607 | -3.844 | 0.001 | 1 | -2.979 |
| 1552414_at | 1552414_at | WFDC9 | 0.382 | 3.463 | 3.803 | 0.001 | 1 | -3.005 |
| 1555235_s_at | 1555235_s_at | IQCF3 | 0.392 | 3.526 | 3.732 | 0.001 | 1 | -3.051 |
| 224467_s_at | 224467_s_at | PDCD2L | -0.813 | 7.225 | -3.726 | 0.001 | 1 | -3.055 |
| 203343_at | 203343_at | UGDH | 0.718 | 7.099 | 3.696 | 0.001 | 1 | -3.075 |
| 215154_at | 215154_at | ULK2 | 0.381 | 4.507 | 3.687 | 0.001 | 1 | -3.080 |
| 218702_at | 218702_at | SARS2 | -0.556 | 5.251 | -3.677 | 0.001 | 1 | -3.087 |
| 200820_at | 200820_at | PSMD8 | -0.692 | 9.705 | -3.669 | 0.001 | 1 | -3.092 |
| 211047_x_at | 211047_x_at | AP2S1 | -0.524 | 9.801 | -3.653 | 0.001 | 1 | -3.102 |
8 关键可视化
8.1 火山图
# 按探索性标准着色（P < 0.05 且 |log2FC| > 1）
tt$sig_class <- factor(
ifelse(tt$P.Value < 0.05 & tt$logFC > 1, "Up in recurrence",
ifelse(tt$P.Value < 0.05 & tt$logFC < -1, "Down in recurrence", "Not significant")),
levels = c("Up in recurrence", "Down in recurrence", "Not significant"))
# 标注最显著的候选基因（前 15 个有基因符号的）
lab <- head(tt[!is.na(tt$symbol) & tt$sig_class != "Not significant", ], 15)
ggplot(tt, aes(logFC, -log10(P.Value), color = sig_class)) +
geom_point(alpha = 0.6, size = 1.3) +
scale_color_manual(values = c("Up in recurrence" = "#FF6E6E",
"Down in recurrence" = "#4FC3F7",
"Not significant" = "grey55")) +
geom_vline(xintercept = c(-1, 1), linetype = 2, color = "grey60") +
geom_hline(yintercept = -log10(0.05), linetype = 2, color = "grey60") +
geom_text_repel(data = lab, aes(label = symbol), size = 3, max.overlaps = 30,
box.padding = 0.4, segment.color = "grey70") +
labs(title = "火山图：术后复发 vs 无复发",
subtitle = sprintf("虚线为 |log2FC| = 1 与 P = 0.05（探索性筛选）；FDR 校正后最小 FDR = %.3f，无显著探针集",
min(tt$adj.P.Val)),
x = expression(log[2]~"fold change (Recurrence / NoRecurrence)"),
y = expression(-log[10]~"(P value)"), color = NULL)
8.2 热图（候选差异基因）
# 取探索性候选基因（|log2FC| > 1 中 P 值最小的若干探针集），按行（基因）做 z-score 标准化
top_n <- 50
hm_probe <- head(cand$probe_id, top_n)
if (length(hm_probe) < 10) hm_probe <- head(subset(tt, abs(logFC) > 1)$probe_id, top_n)
hm_mat <- expr_norm[hm_probe, , drop = FALSE]
hm_lab <- tt$symbol[match(hm_probe, tt$probe_id)]
rownames(hm_mat) <- ifelse(is.na(hm_lab), hm_probe, hm_lab)
hm_mat <- t(scale(t(hm_mat))) # 按基因标准化，突出组间相对高低
ann_col <- data.frame(Recurrence = pheno$recurrence,
ER = pheno$er,
PR = pheno$pr,
HER2 = pheno$her2,
Grade = pheno$grade,
row.names = pheno$gsm)
pheatmap(hm_mat,
color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100),
annotation_col = ann_col,
cluster_rows = TRUE, cluster_cols = TRUE,
show_colnames = FALSE, fontsize_row = 6.5, border_color = NA,
treeheight_row = 25, treeheight_col = 25,
main = sprintf("探索性候选差异基因（%d 个探针集；行为基因的 z-score）", length(hm_probe)))
8.3 代表基因的表达箱线图
top_genes <- head(unique(na.omit(cand$symbol)), 6)
if (length(top_genes) < 2) top_genes <- head(unique(na.omit(tt$symbol)), 6)
gb <- do.call(rbind, lapply(top_genes, function(g) {
if (!g %in% rownames(expr_gene)) return(NULL)
data.frame(symbol = g, expression = expr_gene[g, ], recurrence = pheno$recurrence)
}))
ggplot(gb, aes(recurrence, expression, fill = recurrence)) +
geom_boxplot(outlier.shape = NA, alpha = 0.7, width = 0.6) +
geom_jitter(width = 0.15, size = 1.6, alpha = 0.9) +
facet_wrap(~ symbol, scales = "free_y", nrow = 2) +
scale_fill_manual(values = c(NoRecurrence = "#4FC3F7", Recurrence = "#FF6E6E")) +
labs(title = "代表性差异基因的表达水平", x = NULL, fill = NULL,
y = expression(log[2]~"expression")) +
theme(legend.position = "bottom")
9 结果输出
out_file <- file.path(res_dir, "GSE19697_limma_recurrence_vs_norecurrence.csv")
write.csv(tt, out_file, row.names = FALSE)
cat("差异表达结果已保存：", out_file, "\n")差异表达结果已保存： results/GSE19697_limma_recurrence_vs_norecurrence.csv out_hm <- file.path(res_dir, "GSE19697_top_diff_genes_expression.csv")
hm_expr <- expr_norm[hm_probe, , drop = FALSE]
rownames(hm_expr) <- rownames(hm_mat) # 行名用基因符号（无法注释的探针保留探针 ID）
write.csv(hm_expr, out_hm)
cat("热图所用基因的表达矩阵已保存：", out_hm, "\n")热图所用基因的表达矩阵已保存： results/GSE19697_top_diff_genes_expression.csv 10 稳健性检查与小结
# 对照：MAS 5.0 值 + 分位数标准化 是否给出相似的 logFC？
if (cel_ok) {
expr_alt <- normalizeBetweenArrays(log2(pmax(expr_mas5, 1)), method = "quantile")
tt_alt <- topTable(eBayes(lmFit(expr_alt, design)),
coef = "Recurrence_vs_NoRecurrence", number = Inf, sort.by = "none")
common <- intersect(rownames(tt), rownames(tt_alt))
cmp <- data.frame(rma = tt$logFC[match(common, rownames(tt))],
mas5 = tt_alt$logFC[match(common, rownames(tt_alt))])
rho <- cor(cmp$rma, cmp$mas5, method = "spearman", use = "complete.obs")
cat(sprintf("RMA 与 MAS5+分位数标准化两种策略的 logFC（Spearman）相关系数：%.3f\n", rho))
cat(sprintf(" 两种策略下的候选基因（P<0.05 且 |logFC|>1）重叠比例：%.1f%%\n",
100 * mean(rownames(tt)[tt$P.Value < 0.05 & abs(tt$logFC) > 1] %in%
rownames(tt_alt)[tt_alt$P.Value < 0.05 & abs(tt_alt$logFC) > 1])))
print(ggplot(cmp, aes(rma, mas5)) +
geom_point(alpha = 0.15, size = 0.6) +
geom_abline(slope = 1, intercept = 0, color = "red", linetype = 2) +
labs(title = sprintf("两种标准化策略的 logFC 一致性（rho = %.3f）", rho),
x = "logFC (CEL + RMA)", y = "logFC (MAS5 + quantile)"))
} else {
cat("本次运行未使用 CEL 原始数据，跳过稳健性比较。\n")
}RMA 与 MAS5+分位数标准化两种策略的 logFC（Spearman）相关系数：0.587
两种策略下的候选基因（P<0.05 且 |logFC|>1）重叠比例：64.1%
结果解读要点
- 分组方向：
recur = 1（术后复发）为实验组、recur = 0（无复发）为对照组，logFC > 0表示该基因在复发组中表达更高。 - 标准化是必需的：RMA
已完成背景校正、
log2转换、分位数标准化与探针集汇总，可直接送入 limma；而 series matrix 的 MAS 5.0 值是线性尺度，必须先log2再做分位数标准化，否则会出现大量假阳性。 - 统计模型：
limma的线性模型 + 经验贝叶斯方差收缩，能在 n = 24 的小样本下稳定估计方差，是芯片差异表达分析的标准做法。当前模型只含复发分组，未校正 ER/HER2/分级等混杂（该队列以 basal-like 三阴性为主，异质性相对有限）。 - 结果怎么读：本队列复发组仅 6 例，BH 校正（54 675
次检验）后没有任何探针集达到 FDR < 0.05（最小 FDR ≈
1.0），因此火山图与热图展示的是名义 P
值下的探索性候选基因（P < 0.05 且
|log2FC| > 1）。这类结果只能作为假设生成，需要在更大、更均衡的队列（如 TCGA、METABRIC）中独立验证，不能直接当作临床结论。 - 提升效能的思路：加入 ER/HER2/分级等协变量、改用
eBayes(trend = TRUE)、或针对特定通路做基因集检验（如limma::camera、roast），都比降低阈值更可靠。
sessionInfo()R version 4.6.1 (2026-06-24)
Platform: x86_64-conda-linux-gnu
Running under: Ubuntu 24.04.4 LTS
Matrix products: default
BLAS/LAPACK: /home/wangzg/miniconda3/envs/r-4.6/lib/libopenblasp-r0.3.34.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Asia/Shanghai
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] hgu133plus2cdf_2.18.0 matrixStats_1.5.0 RColorBrewer_1.1-3
[4] pheatmap_1.0.13 ggrepel_0.9.8 ggplot2_4.0.3
[7] hgu133plus2.db_3.13.0 org.Hs.eg.db_3.23.1 AnnotationDbi_1.74.0
[10] IRanges_2.46.0 S4Vectors_0.50.1 affy_1.90.0
[13] Biobase_2.72.0 BiocGenerics_0.58.1 generics_0.1.4
[16] limma_3.68.5
loaded via a namespace (and not attached):
[1] KEGGREST_1.52.2 gtable_0.3.6 xfun_0.60
[4] bslib_0.12.0 lattice_0.23-1 vctrs_0.7.3
[7] tools_4.6.1 tibble_3.3.1 RSQLite_3.53.3
[10] blob_1.3.0 pkgconfig_2.0.3 S7_0.2.2
[13] lifecycle_1.0.5 compiler_4.6.1 farver_2.1.2
[16] textshaping_1.0.5 Biostrings_2.80.1 statmod_1.5.2
[19] Seqinfo_1.2.0 htmltools_0.5.9 sass_0.4.10
[22] yaml_2.3.12 preprocessCore_1.74.0 pillar_1.11.1
[25] crayon_1.5.3 jquerylib_0.1.4 MASS_7.3-66
[28] thematic_0.1.8 cachem_1.1.0 tidyselect_1.2.1
[31] digest_0.6.39 dplyr_1.2.1 labeling_0.4.3
[34] fastmap_1.2.0 grid_4.6.1 cli_3.6.6
[37] magrittr_2.0.5 withr_3.0.3 rappdirs_0.3.4
[40] scales_1.4.0 bit64_4.8.4 rmarkdown_2.31
[43] XVector_0.52.0 httr_1.4.8 affyio_1.82.0
[46] bit_4.6.0 otel_0.2.0 ragg_1.5.2
[49] png_0.1-9 memoise_2.0.1 evaluate_1.0.5
[52] knitr_1.51 rlang_1.3.0 Rcpp_1.1.2
[55] glue_1.8.1 DBI_1.3.0 BiocManager_1.30.27
[58] rstudioapi_0.19.0 jsonlite_2.0.0 R6_2.6.1
[61] systemfonts_1.3.2