GEO 数据挖掘：下载、探针注释、limma 差异分析与常见坑

library(GEOquery)
options(timeout = 600)                        # 默认 60 秒，大文件会中途断开
Sys.setenv(VROOM_CONNECTION_SIZE = 500000L)   # 避免 connection buffer 报错
dir.create("geo_cache", showWarnings = FALSE)

# destdir 固定缓存目录：默认存在 tempdir()，关掉 R 就没了，下次重新下载
gse <- getGEO("GSE16515", destdir = "geo_cache",
              GSEMatrix = TRUE, getGPL = FALSE)  # 不下载 85.6 MB 的 GPL570 注释
length(gse)        # 多平台 Series 返回多个 ExpressionSet，逐个检查
eset <- gse[[1]]
dim(eset)          # 54613 探针 x 52 样本
annotation(eset)   # "GPL570"

# FTP 路径规则：GSE 编号去掉后三位换成 nnn（GSE16515 -> GSE16nnn；GSE1000 以下为 GSEnnn）
# -c 支持断点续传，网络中断后重跑同一命令即可
wget -c https://ftp.ncbi.nlm.nih.gov/geo/series/GSE16nnn/GSE16515/matrix/GSE16515_series_matrix.txt.gz
# 平台注释（GPL570 小于 1000，目录为 GPLnnn）
wget -c https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPLnnn/GPL570/annot/GPL570.annot.gz
# 原始 CEL 文件
wget -c https://ftp.ncbi.nlm.nih.gov/geo/series/GSE16nnn/GSE16515/suppl/GSE16515_RAW.tar

# R 中读本地文件，不再联网
Rscript -e 'library(GEOquery); e <- getGEO(filename = "GSE16515_series_matrix.txt.gz", getGPL = FALSE); print(dim(e))'

library(limma)
ex <- exprs(eset)
# GEO2R 的判断规则
qx <- as.numeric(quantile(ex, c(0, .25, .5, .75, .99, 1), na.rm = TRUE))
LogC <- (qx[5] > 100) || (qx[6] - qx[1] > 50 && qx[2] > 0)
qx; LogC
# GSE16515: 2.28 7.11 18.17 86.87 3283.25 65045 -> TRUE
if (LogC) { ex[ex <= 0] <- NaN; ex <- log2(ex) }

# 看分布：各样本箱线图中位数应大致对齐
boxplot(ex, las = 2, outline = FALSE)
# 明显不齐时再做分位数归一化（GEO2R 的 Force normalization 也是这一步）
# ex <- normalizeBetweenArrays(ex, method = "quantile")

library(hgu133plus2.db)   # GPL570 对应的注释包；其他平台见下表
pk <- AnnotationDbi::select(hgu133plus2.db, keys = rownames(ex),
                            columns = "SYMBOL", keytype = "PROBEID")
multi <- names(which(table(pk$PROBEID) > 1))          # 一个探针对应多个基因：1528 个
pk <- pk[!is.na(pk$SYMBOL) & !pk$PROBEID %in% multi, ] # 去掉无注释和多基因探针
ex <- ex[pk$PROBEID, ]; sym <- pk$SYMBOL               # 剩 43112 个探针、20826 个基因

# 一个基因多个探针：保留平均表达最高的探针（11041 个基因有多个探针）
o <- order(sym, -rowMeans(ex))
keep <- o[!duplicated(sym[o])]
ex_gene <- ex[keep, ]; rownames(ex_gene) <- sym[keep]
# 另一种是取平均：ex_gene <- avereps(ex, ID = sym)

# 没有注释包的平台用 GPL 表（getGEO("GPLxxx") 或 FTP 上的 .annot.gz）
# gpl  <- getGEO("GPL570", destdir = "geo_cache")
# anno <- Table(gpl)[, c("ID", "Gene Symbol")]
# anno <- anno[anno$`Gene Symbol` != "" & !grepl("///", anno$`Gene Symbol`), ]

pd <- pData(eset)
grep(":ch1$", colnames(pd), value = TRUE)   # 分组通常在 xxx:ch1 列
table(pd[["tissue:ch1"]])
# Normal Tissue in Pancreatic Cancer Sample 16 / Tumor Tissue ... 36

group <- factor(ifelse(grepl("Tumor", pd[["tissue:ch1"]]), "Tumor", "Normal"),
                levels = c("Normal", "Tumor"))   # 第一个水平是对照
# 三道核对：样本名一致、交叉表无错配、组内样本数与原文一致
stopifnot(identical(rownames(pd), colnames(ex_gene)))
table(group, pd[["tissue:ch1"]])
patient <- sub("^Pancreatic Sample ?([0-9]*)-.*$", "\\1", pd$title)  # 配对信息，后面用

design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)                       # Normal, Tumor
contr <- makeContrasts(Tumor - Normal, levels = design) # logFC > 0 表示肿瘤中上调
fit  <- lmFit(ex_gene, design)
fit2 <- eBayes(contrasts.fit(fit, contr), trend = TRUE)
tt   <- topTable(fit2, number = Inf)                    # 默认 BH 校正
deg  <- subset(tt, adj.P.Val < 0.05 & abs(logFC) > 1)
table(sign(deg$logFC))                                  # 本页实测：-1 399，1 1136

# 把倍数阈值纳入检验（limma 推荐的做法）
tr <- topTreat(treat(contrasts.fit(fit, contr), lfc = 1, trend = TRUE), number = Inf)
sum(tr$adj.P.Val < 0.05)                                # 本页实测 166

# 配对样本：全部配对时把 patient 放进设计矩阵
# design <- model.matrix(~ patient + group)
# 部分配对（GSE16515 有 16 对加 20 个单独肿瘤）：
corfit <- duplicateCorrelation(ex_gene, design, block = patient)
fitb <- lmFit(ex_gene, design, block = patient, correlation = corfit$consensus)
ttb  <- topTable(eBayes(contrasts.fit(fitb, contr), trend = TRUE), number = Inf)

library(GEOquery); library(edgeR); library(limma); library(org.Hs.eg.db)
acc <- "GSE164073"
hasRNASeqQuantifications(acc)   # TRUE 表示有 NCBI 生成的 counts
# 直接下载 raw counts（这个链接格式在 2026-09-02 后不被 reCAPTCHA 拦截）
url <- paste0("https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=", acc,
              "&format=file&file=", acc, "_raw_counts_GRCh38.p13_NCBI.tsv.gz")
f <- file.path("geo_cache", paste0(acc, "_raw_counts.tsv.gz"))
if (!file.exists(f)) download.file(url, f, mode = "wb")
cnt <- as.matrix(read.delim(f, row.names = 1, check.names = FALSE))  # 行名是 Entrez Gene ID
sym <- mapIds(org.Hs.eg.db, rownames(cnt), "SYMBOL", "ENTREZID")     # 代替被拦截的注释文件

pd  <- pData(getGEO(acc, destdir = "geo_cache", getGPL = FALSE)[[1]])[colnames(cnt), ]
inf <- factor(ifelse(grepl("SARS", pd[["infection:ch1"]]), "CoV2", "mock"), levels = c("mock", "CoV2"))
tis <- factor(pd[["tissue:ch1"]])
design <- model.matrix(~ tis + inf)

y <- DGEList(cnt)
keep <- filterByExpr(y, design = design)          # 39376 -> 16933 个基因
y <- normLibSizes(y[keep, , keep.lib.sizes = FALSE])
v <- voom(y, design)
tt <- topTable(eBayes(lmFit(v, design)), coef = "infCoV2", number = Inf)
tt$SYMBOL <- sym[rownames(tt)]

library(sva)
# m：两个同平台 GSE 合并后的基因矩阵（各自已确认 log2）；batch、group 为因子
table(batch, group)          # 先看设计：每个批次是否都有两组

# 差异分析：batch 放进设计矩阵
design <- model.matrix(~ group + batch)
fit <- eBayes(lmFit(m, design), trend = TRUE)
tt  <- topTable(fit, coef = "groupTumor", number = Inf)

# 画 PCA、热图用的校正矩阵（不再拿去做 limma）
m_plot <- removeBatchEffect(m, batch = batch, design = model.matrix(~ group))
# 一定要用 ComBat 时，带上 mod 保护分组效应
m_cb <- ComBat(m, batch = batch, mod = model.matrix(~ group))

library(ggplot2); library(ggrepel); library(pheatmap)
v <- tt; v$gene <- rownames(v)
v$change <- with(v, ifelse(adj.P.Val < 0.05 & logFC > 1, "Up",
                    ifelse(adj.P.Val < 0.05 & logFC < -1, "Down", "NS")))
lab <- head(v[v$change != "NS", ][order(v$adj.P.Val[v$change != "NS"]), ], 10)
p <- ggplot(v, aes(logFC, -log10(adj.P.Val), colour = change)) +
  geom_point(size = 0.8, alpha = 0.6) +
  scale_colour_manual(values = c(Down = "#2b6cb0", NS = "grey75", Up = "#c53030")) +
  geom_vline(xintercept = c(-1, 1), linetype = 2) +
  geom_hline(yintercept = -log10(0.05), linetype = 2) +   # 纵轴与筛选用同一个 P
  geom_text_repel(data = lab, aes(label = gene), size = 3, colour = "black") +
  labs(x = "log2 fold change (Tumor vs Normal)", y = "-log10 adjusted P") + theme_bw()
ggsave("volcano.png", p, width = 6, height = 5, dpi = 300)

# 热图：上调、下调各取 adj.P 最小的 25 个
up <- head(rownames(v)[v$change == "Up"][order(v$adj.P.Val[v$change == "Up"])], 25)
dn <- head(rownames(v)[v$change == "Down"][order(v$adj.P.Val[v$change == "Down"])], 25)
o  <- order(group)
ann <- data.frame(group = group[o], row.names = colnames(ex_gene)[o])
pheatmap(ex_gene[c(up, dn), o], scale = "row", cluster_cols = FALSE,
         annotation_col = ann, show_colnames = FALSE, fontsize_row = 6,
         breaks = seq(-2, 2, length.out = 101),   # 截断 z 值，少数极端值不压缩色阶
         filename = "heatmap.png", width = 7, height = 8)

library(TCGAbiolinks)
q <- GDCquery(project = "TCGA-PAAD",
              data.category = "Transcriptome Profiling",
              data.type = "Gene Expression Quantification",   # 不写这一行会多出一倍文件
              workflow.type = "STAR - Counts")                # HTSeq - Counts 已不存在
GDCdownload(q, files.per.chunk = 20)   # 分块下载，失败后重跑会跳过已下载的文件
se  <- GDCprepare(q)
cnt <- SummarizedExperiment::assay(se, "unstranded")         # 做 DESeq2/edgeR 用这一层
table(se$sample_type)                  # TCGA-PAAD：Primary Tumor 178，Solid Tissue Normal 4