Transcriptomic Signatures of Renal Cell Carcinoma: A Visual Analytics Case Study

1. Research Question

Cancer biology is fundamentally a problem of gene expression. When a cell becomes malignant, it does not simply acquire a mutation and stop there, it undergoes a wholesale reorganization of its transcriptional state, activating programs that promote uncontrolled growth while silencing the regulatory mechanisms that would normally restrain it. Understanding exactly which genes are affected, by how much, and in what coordinated patterns is the central task of cancer transcriptomics.

This analysis addresses the following research question:

What transcriptomic changes characterize clear cell renal cell carcinoma, and can a small gene expression signature reliably distinguish tumor from normal kidney tissue?

Three specific objectives follow from this question. First, to identify the individual genes most significantly dysregulated in RCC relative to normal kidney tissue and characterize the magnitude and statistical robustness of those differences. Second, to determine whether differentially expressed genes cluster into coherent biological pathways, particularly those consistent with known RCC molecular biology. Third, to evaluate whether a sparse, LASSO-selected gene panel can serve as a candidate diagnostic signature, establishing a connection between population-level transcriptomics and potential clinical utility.

2. Problem Description and Hypothesis

Renal cell carcinoma (RCC) is the most common malignancy of the kidney, accounting for approximately 90 percent of kidney cancers and an estimated 80,000 new diagnoses annually in the United States. The clear cell subtype, which represents roughly 75 percent of RCC cases, is defined at the molecular level by biallelic inactivation of the VHL (Von Hippel-Lindau) tumor suppressor gene on chromosome 3p. VHL normally targets the transcription factor HIF1\(\alpha\) (Hypoxia Inducible Factor 1-alpha) for ubiquitin-mediated degradation under normoxic conditions. When VHL is lost, HIF1\(\alpha\) constitutively accumulates regardless of oxygen availability, driving a broad transcriptional program that mimics chronic hypoxia. This program activates genes involved in angiogenesis (particularly VEGF and its receptors), anaerobic glucose metabolism, and cell survival, processes that collectively favor tumor growth and invasion.

Despite this well-characterized molecular etiology, RCC is frequently diagnosed at advanced stage because early disease is largely asymptomatic. The “incidentaloma” pattern, discovery during imaging performed for unrelated reasons, accounts for a growing proportion of diagnoses, but many patients still present with metastatic disease, for which five-year survival rates remain below 15 percent. Identifying robust molecular signatures from bulk tumor expression data has the potential to improve early detection, subtype classification, and patient stratification for targeted therapy.

The dataset used in this analysis, GDS507 from the NCBI Gene Expression Omnibus, contains expression profiles from 17 clear cell RCC tumor samples and 10 matched normal kidney cortex samples, measured on the Affymetrix HG-U133A microarray platform (~22,000 probe sets). All patients underwent radical or partial nephrectomy, and tumor and normal samples were collected from the same surgical specimen where possible, reducing inter-patient confounding.

The central hypothesis is threefold: (1) clear cell RCC will produce a large, statistically robust differential expression signature relative to normal kidney; (2) that signature will be enriched for gene sets related to hypoxia, angiogenesis, and metabolic reprogramming consistent with HIF1\(\alpha\) pathway activation; and (3) a small number of genes, fewer than 20, will jointly carry sufficient discriminative information to classify tumor from normal tissue with high accuracy.

3. Related Work

The analysis of bulk tumor gene expression data has a rich history in computational biology, and several landmark studies have shaped both the methods and the visual conventions used here.

The Cancer Genome Atlas (TCGA) kidney renal clear cell carcinoma (KIRC) cohort represents the most comprehensive multi-omic characterization of clear cell RCC to date, integrating expression, methylation, copy number, and somatic mutation data across over 500 patients. The TCGA KIRC analysis confirmed the centrality of VHL pathway dysregulation and identified VEGFA upregulation, CA9 (carbonic anhydrase 9) as a near-universal marker of clear cell RCC, and broad downregulation of oxidative phosphorylation genes as consistent molecular features of the disease. This analysis uses the same biological hypotheses, though on a much smaller dataset, making GDS507 a useful proof-of-concept for the methods developed in BioUtils.

The limma/voom framework used here for differential expression was introduced by Smyth (2004) and has become the field standard for microarray analysis specifically because of its empirical Bayes moderation of variance estimates. Its outputs, particularly the volcano plot, have become the canonical visualization for differential expression results across thousands of published papers. The volcano plot is compelling as a visual because it encodes two distinct dimensions of evidence simultaneously: the x-axis captures biological magnitude (fold change) and the y-axis captures statistical confidence (adjusted p-value), so that genes of genuine interest naturally cluster in the visually distinctive upper corners. An interactive version using plotly extends this by allowing hover-based gene identification, which in a static figure would require labeling and quickly becomes cluttered.

Gene Set Enrichment Analysis (GSEA), introduced by Subramanian et al. (2005), addresses a fundamental limitation of single-gene testing: the most statistically significant individual gene is not necessarily the most biologically informative. GSEA operates on the full ranked gene list and asks whether genes from a known biological pathway are systematically enriched at either end of the ranking. The NES bar chart used here to display GSEA results follows the visual convention established by the Broad Institute’s GSEA software and is effective as a part-to-whole visualization because it shows, simultaneously, the direction and relative magnitude of pathway-level responses across the entire Hallmark collection.

The LASSO biomarker approach draws on Tibshirani’s original 1996 paper and its subsequent application to genomic data by Zou and Hastie. The cross-validated LASSO has been used extensively in cancer biomarker discovery because it naturally handles the high-dimensional, correlated predictor structure of expression data and produces interpretable sparse solutions. Importantly, the genes selected by LASSO are not simply the most differentially expressed, they are the genes with non-redundant predictive value given each other, which is the clinically relevant criterion for a diagnostic panel.

4. Data and Methods

library(BioUtils)
library(plotly)
library(msigdbr)
library(ggplot2)

geo <- extract.expression(
  load.geo.soft(accession = "GDS507", log.transform = TRUE)
)

All expression values are log2-transformed at import. Log transformation is standard for Affymetrix arrays because raw intensities span several orders of magnitude and are right-skewed; log2 transformation compresses this range, stabilizes variance, and means that arithmetic differences correspond to fold changes, a one-unit difference represents a doubling of expression.

The dataset contains 17 samples profiled across 22645 probe sets. Sample composition is 9 RCC tumor samples and 8 normal kidney samples.

The analytical workflow proceeds through six stages: multivariate dimensionality reduction for quality control, genome-wide differential expression testing, distribution analysis of top candidate genes, correlation analysis of co-expression structure, pathway enrichment analysis, and multi-gene predictive modeling. Each stage addresses a distinct analytical question and employs a distinct visualization type.

5. Multivariate Analysis: Principal Component Analysis

Principal Component Analysis (PCA) is applied first, before any filtering or testing, as a global quality control step. PCA finds the linear combinations of the 22,000 probe variables that explain the most variance across samples and projects each sample into this reduced space. If the primary axis of variation in the data corresponds to the biological comparison of interest, disease state, it suggests the transcriptional signal is large relative to technical noise and that downstream statistical tests will have sufficient power.

pca.plot(geo$expression, geo$phenotype, color.by = "disease.state")

The PCA plot reveals strong separation between RCC and normal samples along PC1, which explains the majority of variance in the dataset. This clean separation has two important implications. First, it confirms that the dominant source of variation in the expression matrix is biology rather than technical factors such as batch effects, RNA quality differences, or processing date, if technical confounders were present, we would expect samples to cluster by processing order or batch rather than disease state. Second, the magnitude of the PC1 separation suggests that RCC produces a pervasive and consistent transcriptional reprogramming across tumor samples, rather than heterogeneous, sample-specific changes. This bodes well for the subsequent differential expression analysis and for the hypothesis that a small gene panel can reliably classify tumor from normal tissue.

The tight clustering of normal samples relative to the somewhat wider spread of RCC samples along PC2 is also biologically meaningful: normal kidney tissue is a relatively homogeneous cell type (primarily proximal tubule epithelium), whereas clear cell RCC tumors can vary in grade, stage, and the degree of HIF1\(\alpha\) pathway activation, producing slightly more variable expression profiles within the tumor group.

6. Ranking: Genome-Wide Differential Expression

de.results <- run.limma.de(geo, condition.col = "disease.state")

sig <- de.results[de.results$adj.P.Val < 0.05 & abs(de.results$logFC) > 1, ]
cat("Significant DE probes (FDR < 0.05, |logFC| > 1): ", nrow(sig), "\n")
#> Significant DE probes (FDR < 0.05, |logFC| > 1):  1375
cat("Upregulated in RCC:                              ", sum(sig$logFC > 0), "\n")
#> Upregulated in RCC:                               759
cat("Downregulated in RCC:                            ", sum(sig$logFC < 0), "\n")
#> Downregulated in RCC:                             616

Differential expression testing using the limma empirical Bayes framework identifies a substantial number of significantly dysregulated probes at a 5 percent FDR threshold with a minimum 2-fold change. The asymmetry between upregulated and downregulated genes is itself biologically informative. RCC is notable among cancers for producing extensive downregulation of normal kidney function genes alongside the activation of oncogenic programs, the proximal tubule cells that give rise to clear cell RCC normally express high levels of metabolic, transport, and detoxification genes that are systematically silenced during malignant transformation. This pattern of “lineage infidelity,” the loss of the transcriptional identity of the cell of origin, is reflected in the balance of up- and downregulated probes.

7. Deviation Analysis: Interactive Volcano Plot

The volcano plot is the canonical visualization for differential expression results because it simultaneously encodes two orthogonal dimensions of evidence. The horizontal axis represents the magnitude of change (log2 fold change), and the vertical axis represents the statistical confidence in that change (-log10 adjusted p-value). Genes of genuine biological interest occupy the upper corners, large magnitude changes with high statistical confidence, while uninformative probes cluster near the origin. The dashed threshold lines at |logFC| = 1 and FDR = 0.05 define the boundaries used throughout this analysis.

volcano.plot(de.results, fc.threshold = 1, fdr.threshold = 0.05)

The static volcano plot communicates the global structure of the differential expression result, the sheer number of significantly dysregulated probes and the approximate symmetry (or asymmetry) of upregulation versus downregulation. The interactive version below extends this by enabling identification of specific genes by hovering, which is particularly valuable for exploring the biologically interesting probes in the upper corners that would be too numerous to label in a static figure.

de.plot <- de.results
de.plot$neg.log10.fdr <- -log10(de.plot$adj.P.Val)
de.plot$status <- "Not Significant"
de.plot$status[de.plot$adj.P.Val < 0.05 & de.plot$logFC > 1] <- "Upregulated"
de.plot$status[de.plot$adj.P.Val < 0.05 & de.plot$logFC < -1] <- "Downregulated"

de.plot$gene <- ""
top.idx <- order(de.plot$adj.P.Val)[1:50]
de.plot$gene[top.idx] <- get.gene.name(
  geo$gene, rownames(de.plot)[top.idx], use.symbols = TRUE
)

plotly::plot_ly(
  data = de.plot,
  x = ~logFC,
  y = ~neg.log10.fdr,
  color = ~status,
  colors = c("Upregulated" = "firebrick",
             "Downregulated" = "steelblue",
             "Not Significant" = "grey70"),
  text = ~paste0("Gene: ", gene,
                 "<br>logFC: ", round(logFC, 2),
                 "<br>FDR: ", signif(adj.P.Val, 3)),
  hoverinfo = "text",
  type = "scatter",
  mode = "markers",
  marker = list(size = 4, opacity = 0.6)
) %>%
  plotly::layout(
    title = "Volcano Plot: RCC vs Normal Kidney (Interactive, hover for gene labels)",
    xaxis = list(title = "Log2 Fold Change", zeroline = TRUE),
    yaxis = list(title = "-Log10 Adjusted P-Value", zeroline = FALSE),
    shapes = list(
      list(type = "line", x0 = 1,  x1 = 1,
           y0 = 0, y1 = max(de.plot$neg.log10.fdr, na.rm = TRUE),
           line = list(dash = "dot", color = "black")),
      list(type = "line", x0 = -1, x1 = -1,
           y0 = 0, y1 = max(de.plot$neg.log10.fdr, na.rm = TRUE),
           line = list(dash = "dot", color = "black"))
    )
  )