This vignette demonstrates a complete workflow for extracting, filtering, downloading, and analyzing urothelial cancer expression data from GEO. We cover dataset discovery, quality-based filtering, data retrieval, and exploratory analysis with publication-quality visualizations.
Urothelial carcinoma, the most common histological subtype of bladder cancer, is one of the most common urological malignancies. Meta-analysis of gene expression studies can reveal consistent molecular alterations across patient cohorts. This workflow shows how to:
Search for bladder and urothelial cancer datasets using multiple related terms to maximize coverage. Note that “urothelial carcinoma” is the standard histopathological term and matches records the other two terms miss.
uc_searches <- list(
bladder = geo_search(
"bladder cancer[ALL] AND Homo sapiens[ORGN] AND GSE[ETYP]"
),
urothelial = geo_search(
"urothelial cancer[ALL] AND Homo sapiens[ORGN] AND GSE[ETYP]"
),
carcinoma = geo_search(
"urothelial carcinoma[ALL] AND Homo sapiens[ORGN] AND GSE[ETYP]"
)
)
#> → Parsing GEO records
#> → Parsing GEO records
#> → Parsing GEO records
uc_gse <- unique(dplyr::bind_rows(uc_searches))
cat(sprintf("Found %d unique urothelial cancer datasets\n", nrow(uc_gse)))
#> Found 795 unique urothelial cancer datasets
head(uc_gse[, 1:3])
#> Title
#> 1 Transcriptomic profiling of BCG-resistant bladder cancer cells and paired bladder tumor tissues before and after BCG therapy
#> 2 E2F3 amplification primes bladder cancer cells for premature mitosis
#> 3 RIP-seq and N6-methyladenosine RIP-seq (MeRIP-seq) in UM-UC-3 cells
#> 4 ATHENA: Augmenting Therapeutic Effectiveness through Novel Analytics
#> 5 Bulk RNA sequencing of human bladder cancer clinical specimens for transcriptomic analysis.
#> 6 Morphological and Single-Cell Analyses of En Bloc Resected Non-Muscle invasive Bladder Cancer Reveal Structural and functional disruption of Tertiary lymphoid Structures
#> Summary
#> 1 This study investigated transcriptomic alterations associated with Bacillus Calmette-Guérin (BCG) resistance in bladder cancer. Bulk RNA sequencing was performed using paired bladder tumor tissues collected before and after BCG therapy from two patients, as well as parental T24 bladder cancer cells and a BCG-resistant T24-derived cell line (T24-BCG). Differential gene expression analyses were performed to identify transcriptional changes associated with BCG exposure and resistance and to compare molecular alterations between clinical tumor samples and the BCG-resistant cell model.
#> 2 The E2F–RB pathway controls the G1/S checkpoint, and tumors often bypass it through RB1 or CDKN2A loss. Unlike these lesions, E2F3 amplification drives persistent excessive E2F-dependent transcription through S and G2 phases. This oncogene is frequently amplified in for example bladder cancer, but its impact on the cancer cell cycle remains unclear. Using isogenic bladder cancer models and patient data, we show that E2F3 amplification hyperactivates a mitotic gene expression program, including cyclin B1. more...
#> 3 Cisplatin resistance and tumor relapse are the leading causes of the mortality in muscle-invasive bladder cancer (MIBC). Preceding evidence supported that basal/squamous subtype, N6-methyladenosine (m6A) modifications and hypoxia are associated with a reduced chemotherapy response. To explore the m6a modification in UM-UC-3 cells, we performed RNA immunoprecipitation with an antibody against m6a (No. more...
#> 4 With a yearly incidence of about 2500 new cases in Belgium, bladder cancer (BC) is a commonly-diagnosed cancer type. Furthermore, it is accompanied with high mortality rates, with ca. 70% of patients not surviving past five years (Richters et al., 2020). In order to reduce the risk of recurrence and progression, patients with non-muscle invasive bladder cancer (NMIBC) are monitored frequently, undergoing cystoscopy examination and cytology. more...
#> 5 This study characterized transcriptomic features of human bladder cancer using bulk RNA sequencing of tumor specimens. Processed gene-expression data were analyzed to investigate biological features associated with the tumor microenvironment and tumor characteristics.
#> 6 Tertiary lymphoid structures (TLS) have emerged as critical immune niches within the tumor microenvironment across various cancers. However, their structural organization and functional roles in non–muscle-invasive bladder cancer (NMIBC) remain poorly characterized. In this study, we provide the first detailed characterization of TLS in NMIBC using en bloc resected specimens, enabling high-resolution spatial mapping and quantitative histological assessment across a large cohort of primary tumors (n=97). more...
#> Organism
#> 1 Homo sapiens
#> 2 Homo sapiens
#> 3 Homo sapiens
#> 4 Homo sapiens
#> 5 Homo sapiens
#> 6 Homo sapiensThe count parsed below covers every sample in a series regardless of group — cases, controls, and technical replicates — so treat it as an upper bound. Group sizes and balance (e.g., tumor vs. normal) must be verified from the sample annotations before a dataset is actually included in a meta-analysis.
uc_gse <- uc_gse |>
dplyr::mutate(
number_of_samples = str_match(
Contains, "(\\d+) Samples?"
)[, 2L, drop = TRUE],
number_of_samples = as.integer(number_of_samples)
)
cat("\nSample size summary:\n")
#>
#> Sample size summary:
summary(uc_gse$number_of_samples)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 1.0 6.0 12.0 83.9 45.5 16190.0Apply stringent filters to select candidate datasets suitable for meta-analysis.
uc_gse_filtered <- uc_gse |>
dplyr::filter(
# At least 20 samples for robust analysis
number_of_samples >= 20,
# Focus on expression profiling
str_detect(Type, regex("expression profiling", ignore_case = TRUE)),
# Prefer microarray studies (exclude RNA-seq for now)
str_detect(Type, regex("array", ignore_case = TRUE)),
# Exclude methylation studies (defensive: no current GEO Type string
# passes the filters above, but this guards against future variants)
!str_detect(Type, regex("methylation", ignore_case = TRUE)),
# Exclude compound SuperSeries so each study has a single assay type
!str_detect(Type, ";")
) |>
dplyr::arrange(desc(number_of_samples))
cat(sprintf(
"\nFiltered to %d candidate datasets (%.1f%% of original)\n",
nrow(uc_gse_filtered),
100 * nrow(uc_gse_filtered) / nrow(uc_gse)
))
#>
#> Filtered to 95 candidate datasets (11.9% of original)
head(
uc_gse_filtered[, c("Series Accession", "Title", "number_of_samples")],
10
)
#> Series Accession
#> 1 GSE203024
#> 2 GSE2109
#> 3 GSE36133
#> 4 GSE5479
#> 5 GSE68950
#> 6 GSE68928
#> 7 GSE198269
#> 8 GSE48277
#> 9 GSE245953
#> 10 GSE32894
#> Title
#> 1 Gene Expression Data from Human Peripheral Blood
#> 2 Expression Project for Oncology (expO)
#> 3 Expression data from the Cancer Cell Line Encyclopedia (CCLE)
#> 4 Gene expression signatures predict outcome in non-muscle invasive bladder carcinoma - a multi-center validation study
#> 5 caArray_golub-00327: Sanger cell line Affymetrix gene expression project
#> 6 caArray_golub-00236: Multiclass cancer diagnosis using tumor gene expression signatures
#> 7 Predictive biomarkers for survival benefit with ramucirumab in urothelial cancer in the RANGE trial
#> 8 Identification of distinct basal and luminal subsets of human bladder cancers with different sensitivities to frontline chemotherapy
#> 9 Gene expression data from muscle-invasive bladder cancer samples II
#> 10 A Molecular Taxonomy of Urothelial Carcinoma
#> number_of_samples
#> 1 2845
#> 2 2158
#> 3 917
#> 4 808
#> 5 798
#> 6 539
#> 7 394
#> 8 349
#> 9 310
#> 10 308For demonstration, we work with a subset of datasets. In practice,
you would download all filtered datasets, ideally into a persistent
directory instead of tempdir().
demo_gse <- head(uc_gse_filtered[["Series Accession"]], 3)
cat("Demo datasets:", paste(demo_gse, collapse = ", "), "\n")
#> Demo datasets: GSE203024, GSE2109, GSE36133expression_sets <- list()
for (gse_id in demo_gse) {
message(sprintf("Downloading %s...", gse_id))
expression_sets[[gse_id]] <- tryCatch(
geo_matrix(gse_id, odir = tempdir()),
error = function(e) {
message(sprintf("Failed to download %s: %s", gse_id, conditionMessage(e)))
NULL
}
)
}
# Drop failed downloads; names are kept as the GSE accessions
expression_sets <- Filter(Negate(is.null), expression_sets)
cat(sprintf("Successfully downloaded %d datasets\n", length(expression_sets)))Apply log transformation to ensure data comparability across
datasets. log_trans() applies log2(x + 1) only
when a GEO2R-style heuristic indicates the data are not already on a log
scale — check the printed ranges to confirm each dataset was handled as
intended.
Examine the distribution of expression values across samples using boxplots.
eset <- expression_sets_log[[1]]
expr_matrix <- Biobase::exprs(eset)
expr_df <- as.data.frame(expr_matrix) |>
dplyr::mutate(gene = rownames(expr_matrix)) |>
tidyr::pivot_longer(
cols = -gene,
names_to = "sample",
values_to = "expression"
)
ggplot(expr_df, aes(x = sample, y = expression)) +
geom_boxplot(fill = "lightblue", outlier.size = 0.5, outlier.alpha = 0.3) +
labs(
title = sprintf(
"Expression Distribution - %s",
names(expression_sets_log)[1]
),
x = "Sample",
y = "Expression (log2 scale)"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.major.x = element_blank()
)Perform PCA to visualize sample relationships and identify potential batch effects.
eset <- expression_sets_log[[1]]
expr_matrix <- Biobase::exprs(eset)
# Remove genes with missing values
expr_matrix <- expr_matrix[complete.cases(expr_matrix), ]
pca_result <- stats::prcomp(t(expr_matrix), scale. = TRUE)
pca_df <- as.data.frame(pca_result$x[, 1:2])
pca_df$sample <- rownames(pca_df)
var_explained <- summary(pca_result)$importance[2, 1:2] * 100
ggplot(pca_df, aes(x = PC1, y = PC2)) +
geom_point(size = 3, alpha = 0.7, color = "darkred") +
geom_text(aes(label = sample), size = 2.5, vjust = -1) +
labs(
title = sprintf("PCA - %s", names(expression_sets_log)[1]),
x = sprintf("PC1 (%.1f%% variance)", var_explained[1]),
y = sprintf("PC2 (%.1f%% variance)", var_explained[2])
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid = element_line(color = "gray90")
)Calculate and visualize pairwise sample correlations to assess data quality and identify outliers.
eset <- expression_sets_log[[1]]
expr_matrix <- Biobase::exprs(eset)
expr_matrix <- expr_matrix[complete.cases(expr_matrix), ]
cor_matrix <- cor(expr_matrix, method = "pearson")
cor_df <- as.data.frame(as.table(cor_matrix))
colnames(cor_df) <- c("Sample1", "Sample2", "Correlation")
ggplot(cor_df, aes(x = Sample1, y = Sample2, fill = Correlation)) +
geom_tile() +
scale_fill_gradient2(
low = "blue", mid = "white", high = "red",
# Midpoint at the typical sample-correlation level on log expression;
# adjust to the observed range of your data
midpoint = 0.8,
limits = c(min(cor_df$Correlation), 1),
name = "Pearson\nCorrelation"
) +
labs(
title = sprintf(
"Sample Correlation Heatmap - %s",
names(expression_sets_log)[1]
),
x = "Sample",
y = "Sample"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
axis.text.y = element_text(size = 8),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid = element_blank()
)Generate comprehensive summary statistics across datasets.
summary_stats <- lapply(names(expression_sets_log), function(name) {
eset <- expression_sets_log[[name]]
expr_matrix <- Biobase::exprs(eset)
data.frame(
dataset = name,
n_samples = ncol(expr_matrix),
n_genes = nrow(expr_matrix),
mean_expr = mean(expr_matrix, na.rm = TRUE),
median_expr = median(expr_matrix, na.rm = TRUE),
sd_expr = sd(as.vector(expr_matrix), na.rm = TRUE),
missing_pct = mean(is.na(expr_matrix)) * 100
)
})
summary_table <- dplyr::bind_rows(summary_stats)
print(summary_table)
write.csv(summary_table, "uc_expression_summary.csv", row.names = FALSE)Probe identifiers are not gene identifiers.
geo_matrix() returns series matrices whose rows are
platform probe IDs (check Biobase::annotation(eset) and the
featureData slot). Intersecting rownames is only valid when
all datasets come from the same — or a nested — platform. Across mixed
platforms (e.g., Affymetrix and Illumina), map probes to genes with the
platform annotation first, then intersect on gene identifiers.
# 1. Inspect the platform annotation of each dataset, then find common
# features. Remember: intersecting rownames only works across the same
# (or nested) platforms; for mixed platforms, map probes to genes via
# the platform annotation (featureData / Bioconductor annotation
# packages) first.
sapply(expression_sets_log, Biobase::annotation)
common_genes <- Reduce(intersect, lapply(expression_sets_log, function(eset) {
rownames(Biobase::exprs(eset))
}))
cat(sprintf("Common features across all datasets: %d\n", length(common_genes)))
# 2. Subset to common features
expression_sets_common <- lapply(expression_sets_log, function(eset) {
eset[common_genes, ]
})
# 3. Combine expression matrices
combined_expr <- do.call(cbind, lapply(expression_sets_common, Biobase::exprs))
cat(sprintf(
"Combined matrix dimensions: %d genes x %d samples\n",
nrow(combined_expr),
ncol(combined_expr)
))
# 4. Create study labels for batch correction
study_labels <- unlist(lapply(seq_along(expression_sets_common), function(i) {
rep(
names(expression_sets_common)[i],
ncol(Biobase::exprs(expression_sets_common[[i]]))
)
}))
# At this point, you would apply batch correction methods
# such as ComBat (sva package) or limma's removeBatchEffectCompare expression patterns across multiple datasets.
mean_expr_per_dataset <- sapply(expression_sets_log, function(eset) {
rowMeans(Biobase::exprs(eset), na.rm = TRUE)
})
mean_expr_df <- as.data.frame(mean_expr_per_dataset) |>
dplyr::mutate(gene = rownames(mean_expr_per_dataset)) |>
tidyr::pivot_longer(
cols = -gene,
names_to = "dataset",
values_to = "mean_expression"
)
ggplot(mean_expr_df, aes(x = mean_expression, color = dataset)) +
geom_density(linewidth = 1, alpha = 0.7) +
labs(
title = "Mean Expression Distribution Across Datasets",
x = "Mean expression (log2 scale)",
y = "Density",
color = "Dataset"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "bottom"
)With the downloaded and processed data, you can proceed to:
This workflow demonstrated:
geo_matrix()log_trans()These steps provide a solid foundation for conducting rigorous gene expression meta-analysis in urothelial cancer research.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.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: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tidyr_1.3.2 ggplot2_4.0.3 stringr_1.6.0 dplyr_1.2.1 geokit_0.0.2
#> [6] rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 jsonlite_2.0.0 compiler_4.6.1 tidyselect_1.2.1
#> [5] rentrez_1.2.4 jquerylib_0.1.4 scales_1.4.0 yaml_2.3.12
#> [9] fastmap_1.2.0 R6_2.6.1 generics_0.1.4 curl_8.1.0
#> [13] knitr_1.52 XML_3.99-0.25 tibble_3.3.1 maketools_1.3.2
#> [17] bslib_0.12.0 pillar_1.11.1 RColorBrewer_1.1-3 rlang_1.3.0
#> [21] cachem_1.1.0 stringi_1.8.9 xfun_0.62 S7_0.2.2
#> [25] sass_0.4.10 sys_3.4.3 otel_0.2.0 cli_3.6.6
#> [29] withr_3.0.3 magrittr_2.0.5 digest_0.6.39 grid_4.6.1
#> [33] lifecycle_1.0.5 vctrs_0.7.3 evaluate_1.0.5 glue_1.8.1
#> [37] farver_2.1.2 codetools_0.2-20 buildtools_1.0.0 purrr_1.2.2
#> [41] httr_1.4.9 tools_4.6.1 pkgconfig_2.0.3 htmltools_0.5.9