GEO Data Mining: Download, Probe Annotation, limma Differential Expression and Common Pitfalls

Choose the data

GSE, GSM, GPL, GDS and which file to download

A GSE is a study (Series), a GSM is one sample in it, a GPL is the array or sequencing platform, and a GDS is a dataset GEO curated by hand in earlier years. For differential expression the unit you download is the GSE; platform details are in the GPL.

Pipeline for NCBI-generated counts: HISAT2 alignment to GRCh38 (GCA_000001405.15), featureCounts, gene annotation Annotation Release 109.20190905. Runs with alignment rate below 50% and single-cell data are skipped; counts from several runs of one GSM are summed. GEO2R uses these counts with DESeq2 for RNA-seq.

Species: the official page (revised 2026-07-08) still says mouse data are "expected to become available in 2025". An E-utilities search for "rnaseq counts"[Filter] on 2026-10-11 returned 0 mouse-only Series; about 27,700 records carry the filter, all human.

GDS status: counting by year with E-utilities gives 279 new GDS in 2014, 64 in 2015, 2 in 2016 and 0 from 2018 on. Following old tutorials with getGEO("GDSxxxx") only reaches old data.

Download

Downloading with GEOquery, and what to do when it is slow or fails

The code below downloaded the GSE16515 Series Matrix (6.4 MB) in 8.3 s on the test machine. With getGPL = TRUE it also fetches the 85.6 MB GPL570 annotation; on poor networks most failures happen at that step.

A mirror widely used in China is geoChina() from the AnnoProbe package (Biotrainee), which downloads the submitter-processed ExpressionSet (.Rdata) from a Tencent Cloud server. The server was still reachable when tested on 2026-10-11. Its GSE list is built into the package, whose repository was last updated in November 2022, so datasets released after that are not mirrored.

GEOquery versions: this page used GEOquery 2.74.0 (Bioconductor 3.20). The development version on GitHub (2.81.x) plans to change getGEO's return type from ExpressionSet to SummarizedExperiment; code using exprs() and pData() will then need assay() and colData().

library(GEOquery)
options(timeout = 600)                        # default is 60 s; large files get cut off
Sys.setenv(VROOM_CONNECTION_SIZE = 500000L)   # avoids the connection buffer error
dir.create("geo_cache", showWarnings = FALSE)

# Fixed cache directory: the default is tempdir(), which is wiped when R exits
gse <- getGEO("GSE16515", destdir = "geo_cache",
              GSEMatrix = TRUE, getGPL = FALSE)  # skip the 85.6 MB GPL570 annotation
length(gse)        # multi-platform Series return several ExpressionSets
eset <- gse[[1]]
dim(eset)          # 54613 probes x 52 samples
annotation(eset)   # "GPL570"

# FTP path rule: replace the last three digits of the GSE number with nnn
# (GSE16515 -> GSE16nnn; accessions below 1000 use GSEnnn)
# -c resumes partial downloads; rerun the same command after an interruption
wget -c https://ftp.ncbi.nlm.nih.gov/geo/series/GSE16nnn/GSE16515/matrix/GSE16515_series_matrix.txt.gz
# Platform annotation (GPL570 is below 1000, so the folder is GPLnnn)
wget -c https://ftp.ncbi.nlm.nih.gov/geo/platforms/GPLnnn/GPL570/annot/GPL570.annot.gz
# Raw CEL files
wget -c https://ftp.ncbi.nlm.nih.gov/geo/series/GSE16nnn/GSE16515/suppl/GSE16515_RAW.tar

# Read the local file in R without going online
Rscript -e 'library(GEOquery); e <- getGEO(filename = "GSE16515_series_matrix.txt.gz", getGPL = FALSE); print(dim(e))'

Preprocessing

Log2 or not: decide from the values, not the description

limma expects log-scale input. Series Matrix values may already be log2 or may be linear, and the sample description is not reliable: the data_processing field of GSE16515 says GC-RMA (which normally outputs log2), yet the matrix is on a linear scale with a maximum of 65045.

GEO2R's auto-detection only looks at samples assigned to groups and transforms the whole matrix or none of it. The log2(dat + 1) line common in shared templates runs unconditionally; on data already in log2 it compresses fold changes to nearly zero.

If a dataset contains negative values or many zeros (common with MAS5 background subtraction or unlogged Illumina data), a direct log2 produces NaN. Decide after the check whether to reprocess from the raw files.

library(limma)
ex <- exprs(eset)
# The rule GEO2R uses
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) }

# Check distributions: sample medians in the boxplot should roughly line up
boxplot(ex, las = 2, outline = FALSE)
# Only if they clearly do not: quantile normalization (what GEO2R's Force normalization does)
# ex <- normalizeBetweenArrays(ex, method = "quantile")

Annotation

Probe annotation and genes with several probes

The Bioconductor annotation package for a platform is more current than the GPL table. On GPL570 the two disagree for 3874 of 41,039 single-gene probes, mostly because genes were renamed: FAM122C is now PABIR3 and CCDC11 is now CFAP53.

library(hgu133plus2.db)   # annotation package for GPL570; see the table for other platforms
pk <- AnnotationDbi::select(hgu133plus2.db, keys = rownames(ex),
                            columns = "SYMBOL", keytype = "PROBEID")
multi <- names(which(table(pk$PROBEID) > 1))          # probes mapping to several genes: 1528
pk <- pk[!is.na(pk$SYMBOL) & !pk$PROBEID %in% multi, ] # drop unannotated and multi-gene probes
ex <- ex[pk$PROBEID, ]; sym <- pk$SYMBOL               # 43112 probes, 20826 genes left

# Several probes per gene: keep the probe with the highest mean (11041 genes have >1 probe)
o <- order(sym, -rowMeans(ex))
keep <- o[!duplicated(sym[o])]
ex_gene <- ex[keep, ]; rownames(ex_gene) <- sym[keep]
# Alternative, averaging: ex_gene <- avereps(ex, ID = sym)

# Platforms without a package: use the GPL table (getGEO("GPLxxx") or .annot.gz on FTP)
# gpl  <- getGEO("GPL570", destdir = "geo_cache")
# anno <- Table(gpl)[, c("ID", "Gene Symbol")]
# anno <- anno[anno$`Gene Symbol` != "" & !grepl("///", anno$`Gene Symbol`), ]

Groups

Taking groups from pData: the step most likely to fail silently

In the ExpressionSet returned by getGEO, pData rows are in the same order as the expression matrix columns. Errors arise after sorting pData, merging it with an external clinical table or typing a group vector by hand. R raises no error, yet the results change completely.

If sample counts do not match the paper, check whether the GSE spans several platforms (length(gse) > 1) or mixes in cell lines or technical replicates. The 78 samples of GSE15471 come from 36 patients: the 6 samples with _rep in the title are replicate arrays of normal and tumor tissue from 3 patients. Treating replicate arrays as independent samples overstates sample size; merge them by patient and tissue with avereps, or use duplicateCorrelation with patient as the block.

pd <- pData(eset)
grep(":ch1$", colnames(pd), value = TRUE)   # group labels usually live in xxx:ch1 columns
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"))   # first level is the reference
# Three checks: sample names match, no mismatches in the cross-table, group sizes match the paper
stopifnot(identical(rownames(pd), colnames(ex_gene)))
table(group, pd[["tissue:ch1"]])
patient <- sub("^Pancreatic Sample ?([0-9]*)-.*$", "\\1", pd$title)  # pairing, used later

Differential expression

limma: design matrix, contrast matrix and cutoffs

A two-group comparison can be written two equivalent ways: ~0 + group with makeContrasts, or ~group and coef = 2. On GSE15471 the largest logFC difference between the two was 8.8 × 10⁻¹⁵, so use whichever you can read. For more than two groups, use ~0 + group and spell out each contrast.

Basis for the cutoffs: adj.P controls the false discovery rate with the BH method, and 0.05 is GEO2R's default. |logFC| > 1 is a convention without statistical grounding. The limma topTable help states that when p-values and fold changes are not highly correlated, filtering by p-value and then fold change can push the false discovery rate above the nominal level, and recommends treat when a fold-change threshold is wanted. If you use the conventional cutoff, state it and provide results for all genes in the supplement.

DEG counts depend strongly on sample size: for the same tissue, GSE15471 (78 samples) gives 2306 up-regulated genes and GSE16515 (52 samples) gives 1136. When up and down counts are unbalanced, treat up- and down-regulated genes separately in heatmaps and enrichment.

design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)                       # Normal, Tumor
contr <- makeContrasts(Tumor - Normal, levels = design) # logFC > 0 means up in tumor
fit  <- lmFit(ex_gene, design)
fit2 <- eBayes(contrasts.fit(fit, contr), trend = TRUE)
tt   <- topTable(fit2, number = Inf)                    # BH adjustment by default
deg  <- subset(tt, adj.P.Val < 0.05 & abs(logFC) > 1)
table(sign(deg$logFC))                                  # measured here: -1 399, 1 1136

# Build the fold-change threshold into the test (the approach limma recommends)
tr <- topTreat(treat(contrasts.fit(fit, contr), lfc = 1, trend = TRUE), number = Inf)
sum(tr$adj.P.Val < 0.05)                                # measured here: 166

# Paired samples: if every sample is paired, put patient in the design
# design <- model.matrix(~ patient + group)
# Partly paired (GSE16515 has 16 pairs plus 20 unpaired tumors):
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)

RNA-seq

limma for arrays, counts for RNA-seq: choosing the method

RNA-seq differential expression needs raw integer counts. If GEO only has FPKM or TPM, first check whether NCBI generated raw counts; otherwise re-quantify from SRA fastq files.

Measured on 2026-10-11: getRNASeqData("GSE164073") in GEOquery 2.74.0 fails with invalid 'row.names' length. The download link for the human annotation file Human.GRCh38.p13.annot.tsv.gz has no acc parameter and is still blocked by reCAPTCHA, so GEOquery reads the challenge page HTML. The raw count matrix itself downloads fine; annotating by Entrez ID with org.Hs.eg.db maps 37,692 of 39,376 GeneIDs to gene symbols.

This dataset has 18 samples (3 ocular surface tissues × infected/mock × 3 replicates) with design ~ tissue + infection. filterByExpr keeps 16,933 genes, and voom gives 54 up- and 118 down-regulated genes (adj.P < 0.05, |logFC| > 1) in 2.4 s.

library(GEOquery); library(edgeR); library(limma); library(org.Hs.eg.db)
acc <- "GSE164073"
hasRNASeqQuantifications(acc)   # TRUE means NCBI-generated counts exist
# Download raw counts directly (this URL form is exempt from reCAPTCHA since 2026-09-02)
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))  # row names are Entrez Gene IDs
sym <- mapIds(org.Hs.eg.db, rownames(cnt), "SYMBOL", "ENTREZID")     # replaces the blocked annotation file

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 genes
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)]

Merging datasets

Merging several GSEs: put batch in the design matrix

Before merging, confirm the platforms match and each dataset is on the log2 scale. Batch effects are often larger than the biology: after merging GSE16515 and GSE15471 (both GPL570), PC1 explains 48.9% of the variance and mostly separates the datasets.

The limma help for removeBatchEffect states that the function is intended for plotting and data exploration, not for preparing data for lmFit; for linear modelling the batch should be included in the model so that lmFit can assess standard errors correctly. Correcting first and testing afterwards overstates the residual degrees of freedom and makes p-values too small.

ComBat without mod removes group differences that overlap with batch; in the balanced design here it lost 86 up-regulated genes. Nygaard et al. 2016 showed that with unbalanced batch and group, ComBat with mod makes downstream tests overconfident and increases false positives. The difference was small in the unbalanced case here: with group labels permuted 5 times, ComBat yielded 0 to 1 false-positive genes.

If one dataset contains only tumors and the other only normals, batch and group are completely confounded and no method can separate them. All group differences then come from dataset differences, and such a merged comparison should not be done; a reviewer will spot it from table(batch, group).

library(sva)
# m: merged gene matrix of two same-platform GSEs (each confirmed log2); batch and group are factors
table(batch, group)          # check the design first: does every batch contain both groups?

# Differential expression: put batch in the design
design <- model.matrix(~ group + batch)
fit <- eBayes(lmFit(m, design), trend = TRUE)
tt  <- topTable(fit, coef = "groupTumor", number = Inf)

# Corrected matrix for PCA and heatmaps only (not fed back into limma)
m_plot <- removeBatchEffect(m, batch = batch, design = model.matrix(~ group))
# If you must use ComBat, pass mod to protect the group effect
m_cb <- ComBat(m, batch = batch, mod = model.matrix(~ group))

Plots

Volcano plots and heatmaps: code and practical tips

Use the same p-value on the volcano y-axis as in the filter (adj.P here). In heatmaps, order samples by group and scale genes by row.

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) +   # same P on the axis as in the filter
  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)

# Heatmap: the 25 up- and 25 down-regulated genes with the smallest adj.P
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),   # cap z-scores so outliers do not flatten the scale
         filename = "heatmap.png", width = 7, height = 8)

TCGA comparison

How TCGA differs from GEO, and the state of TCGAbiolinks downloads

GEO is a collection of independent studies submitted by many labs, with different platforms and processing; TCGA is a uniformly processed cancer cohort with full clinical and survival data but very few normal controls.

GDC changes: HTSeq counts were removed in Data Release 32 (2022-03-29), leaving STAR - Counts; the Legacy Archive was retired on 2023-05-04, so old code using legacy = TRUE and hg19 no longer works. A GDC API query on 2026-10-11 reported Data Release 46.0 (2026-08-10).

Known TCGAbiolinks issues (GitHub issues still open in 2025–2026): GDCprepare in 2.37.x and 2.38.0 fails on some projects with Column `disease_response` doesn't exist, while users report 2.34 works; GDCquery_clinic broke in early 2025 after GDC API changes, and the maintainer suggested the GitHub devel branch; GDCdownload occasionally fails with Error in if (ret == 1) break (issue #656, unresolved), and method = "client", which calls the GDC Data Transfer Tool, is an option to retry. If downloads keep failing, UCSC Xena provides curated GDC TCGA expression and clinical matrices.

Because TCGA has few normal samples, a common approach is to add GTEx normal tissue. The two differ in sequencing and processing, so batch and group are completely confounded after merging, as explained in the previous section; conclusions from such comparisons need validation in GEO datasets or experiments.

library(TCGAbiolinks)
q <- GDCquery(project = "TCGA-PAAD",
              data.category = "Transcriptome Profiling",
              data.type = "Gene Expression Quantification",   # without this line you get twice the files
              workflow.type = "STAR - Counts")                # HTSeq - Counts no longer exists
GDCdownload(q, files.per.chunk = 20)   # chunked; rerunning skips files already downloaded
se  <- GDCprepare(q)
cnt <- SummarizedExperiment::assay(se, "unstranded")         # use this assay for DESeq2/edgeR
table(se$sample_type)                  # TCGA-PAAD: Primary Tumor 178, Solid Tissue Normal 4

Peer review

Common reviewer objections and how to answer them

Most rejections of GEO data-mining papers come down to the points below. Address what you can before submission.

Practice in China

Verified tips from Chinese community posts

The points below come from posts on CSDN, Jianshu and Biotrainee, each backed by error output, source code or reproduction on this page.

Hand off to an agent

Hand this workflow off to an agent

Download, annotation, differential expression and plotting follow fixed steps and suit a science agent; choosing datasets, defining groups and interpreting results remain your job.

Example one-line instruction: "Using GSE16515 and GSE15471, compare pancreatic tumor and normal tissue: check each for log2, annotate with hgu133plus2.db keeping the highest-mean probe, run limma on each, then merge the matrices with batch in the design matrix; cutoff adj.P < 0.05, |logFC| > 1; output the DEG overlap between datasets, PCA before and after batch correction, a volcano plot and a heatmap of 25 up- and 25 down-regulated genes."

1. 01

   Download and check

   Download and cache the Series Matrix with GEOquery, switching to FTP if the endpoint is blocked; record sample counts, platform and the data_processing description.

2. 02

   Preprocess

   Apply the GEO2R log2 rule, draw boxplots, annotate probes and report probe and gene counts before and after annotation.

3. 03

   Test and plot

   Extract groups and verify sample order, run limma, and output the all-gene results table, volcano plot and heatmap.

4. 04

   Merge and validate

   Merge datasets with batch in the design matrix; output PCA and the DEG overlap table for the two datasets.

5. 05

   Adversarial review

   Check that groups match the paper, logFC signs are correct, already-logged data were not logged again, and the heatmap includes both up- and down-regulated genes.

• The workspace keeps the cached raw files, R scripts, sessionInfo, results tables and all figures.
• You still need to check that the datasets fit your research question, that groups and sample exclusions match the paper, and that the cutoffs suit your paper.
• Experimental validation such as qPCR and immunohistochemistry is up to you.