Cancer Transcriptomics: Differential Gene Expression Analysis in TCGA-BRCA

End-to-end differential gene expression (DGE) analysis pipeline applied to RNA-Seq data from TCGA-BRCA (The Cancer Genome Atlas — Breast Invasive Carcinoma), comparing tumor and adjacent normal breast tissue to identify genes and biological pathways altered in breast cancer.

Aim and Purpose

Breast cancer is driven by widespread dysregulation of gene expression relative to normal tissue. This project aims to:

• Build a reproducible pipeline that takes raw, publicly available RNA-Seq count data from TCGA and turns it into a statistically sound list of differentially expressed genes (DEGs) between tumor and normal breast tissue.
• Characterize the global structure of the expression data (sample composition, count distribution, PCA) before running formal statistical tests, to catch batch effects, class imbalance, or normalization issues early.
• Identify genes that are significantly up- or down-regulated in tumor tissue using rigorous multiple-testing correction, and visualize them (volcano plot, heatmap).
• Interpret the resulting DEGs biologically through pathway/GO enrichment, connecting statistical output back to known hallmarks of cancer (proliferation, loss of differentiation, cell adhesion, etc.).
• Serve as a portfolio-grade demonstration of practical bioinformatics/data-analysis skills: public data acquisition, large-scale data wrangling, statistical inference, and scientific visualization.

Dataset

• Source: UCSC Xena GDC hub (https://gdc.xenahubs.net), mirroring NCI Genomic Data Commons (GDC) — GDC TCGA Breast Cancer (BRCA) cohort.
• Expression data: TCGA-BRCA.star_counts.tsv — gene expression RNA-Seq, STAR-aligned counts. 60,661 gene identifiers × 1,226 samples. Values are provided by Xena already in log2(count + 1) units. Gene/ID mapping uses the GENCODE v36 annotation GTF probemap.
• Clinical/metadata: TCGA-BRCA.clinical.tsv — phenotype data, 1,255 samples × 85 identifiers, used to derive each sample's tissue type (Tumor vs. Normal) from the tissue_type.samples field.
• Cohort size (post cleaning): 1,118 tumor samples and 137 normal samples available in the full clinical file. For computational efficiency, an analysis subset of 200 samples (100 tumor / 100 normal) was randomly selected (random_state=42); after matching against the columns actually present in the counts matrix, the final working cohort used was 100 tumor vs. 83 normal samples (183 total).

Data Acquisition

This analysis uses TCGA-BRCA RNA-Seq data pulled from UCSC Xena's GDC hub. Because of file size, raw data is not included in this repository — download it yourself before running the notebook:

1. Download the two required files from UCSC Xena:
   • TCGA-BRCA.star_counts.tsv (gzipped, ~737 MB) — https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.star_counts.tsv.gz
   • TCGA-BRCA.clinical.tsv (gzipped, ~2 MB) — https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.clinical.tsv.gz
   • Both datasets can also be browsed directly on the GDC TCGA Breast Cancer (BRCA) cohort page on UCSC Xena, which links to full metadata for each dataset.
2. Decompress the .gz files and place both .tsv files in a data/ folder at the repository root:

   data/
   ├── TCGA-BRCA.star_counts.tsv
   └── TCGA-BRCA.clinical.tsv
   

3. Run dge_analysis.ipynb top to bottom. All figures and tables are written automatically to results/figures/, results/tables/, and results/enrichr/.

Pipeline / Methodology

1. Data loading & sample selection

   • Load clinical metadata and derive condition (Tumor/Normal) per sample.
   • Randomly subsample up to 100 samples per group for a balanced, tractable cohort.
   • Load only the matching columns from the (large) counts matrix, using float32 to control memory usage.

2. Filtering & normalization

   • Remove genes with zero counts across all samples.
   • Keep only genes with count > 10 in at least 10% of samples (low-expression filtering).
   • Apply a log2(count + 1) transformation ahead of parametric testing.

3. Exploratory Data Analysis (EDA)

   • Sample composition plot — verifies group balance between Tumor and Normal before testing.
   • Expression distribution plot — confirms the expected right-skewed, negative-binomial-like shape of RNA-Seq counts and justifies the log-transform.
   • PCA (2 components) on standardized, log-transformed expression — checks whether Tumor and Normal samples separate along the top principal components and screens for batch effects.

4. Differential gene expression testing

   • Welch's two-sample t-test (unequal variances) run independently per gene, tumor vs. normal.
   • Effect size measured as log2FC = mean(log2 tumor) − mean(log2 normal).
   • Multiple-testing correction via Benjamini-Hochberg FDR (statsmodels.multipletests).
   • Genes classified as Upregulated / Downregulated / Not Significant using thresholds |log2FC| > 1.0 and FDR (adjusted p) < 0.05.
   • Results exported to results/tables/tcga_brca_degs.csv.

5. Visualization

   • Volcano plot (log2FC vs. −log10 adjusted p-value) highlighting up/downregulated genes.
   • Heatmap of the top 20 DEGs (by absolute log2FC), z-score normalized across samples and ordered Tumor | Normal, to visually confirm the genes stratify patients into the correct groups.

6. Pathway / functional enrichment

   • Up- and down-regulated gene lists submitted separately to gseapy's Enrichr interface against the GO_Biological_Process_2021 gene set to interpret DEGs in terms of biological processes.

Key Results

Metric	Result	Notes
Working cohort	100 Tumor / 83 Normal (183 samples)	Subsampled from 1,118 Tumor / 137 Normal available
Genes tested (post-filtering)	11,245	Filtered from 60,660 total Ensembl gene IDs
Significant DEGs	45 total	log2FC > 1.0, FDR < 0.05
— Upregulated	5 genes
— Downregulated	40 genes
PCA separation	Clear separation	Tumor vs. Normal cluster distinctly along PC1
Heatmap stratification	Clean	Top 20 DEGs correctly separate tumor from normal samples
GO pathway enrichment	Inconclusive	No significant GO Biological Process terms returned at the 0.10 cutoff, likely due to the small number of DEGs (esp. only 5 upregulated genes)

Biological interpretation

The identified DEGs align with established hallmarks of cancer:

• Upregulated genes are consistent with pathways driving uncontrolled proliferation, mitosis, and DNA replication.
• Downregulated genes are consistent with loss of normal breast tissue differentiation, reduced cell adhesion, and loss of tumor-suppressive function.

Repository Structure

.
├── data/                             # not included in repo — see "Data Acquisition"
│   ├── TCGA-BRCA.clinical.tsv        # Clinical/sample metadata from UCSC Xena / GDC
│   └── TCGA-BRCA.star_counts.tsv     # STAR gene-level counts matrix (genes x samples)
├── dge_analysis.ipynb                # Main analysis notebook (this pipeline)
├── results/
│   ├── figures/
│   │   ├── sample_composition.png
│   │   ├── expression_distribution.png
│   │   ├── pca.png
│   │   ├── volcano_plot.png
│   │   ├── heatmap_top20_degs.png
│   │   ├── pathways_upregulated.png      # generated only if enrichment succeeds
│   │   └── pathways_downregulated.png    # generated only if enrichment succeeds
│   ├── tables/
│   │   └── tcga_brca_degs.csv        # full DEG results (all tested genes)
│   └── enrichr/
│       ├── up/                       # raw Enrichr output, upregulated genes
│       └── down/                     # raw Enrichr output, downregulated genes
└── README.md

Technical Stack

• Language: Python 3
• Data handling: pandas, numpy
• Statistics: scipy.stats (Welch's t-test), statsmodels (Benjamini-Hochberg FDR correction)
• Dimensionality reduction: scikit-learn (PCA, StandardScaler)
• Visualization: matplotlib, seaborn
• Pathway enrichment: gseapy (Enrichr API, GO Biological Process 2021 gene set)
• Environment: Jupyter Notebook

How to Reproduce

1. Follow the Data Acquisition steps above to download and place the two required files in data/.
2. Install dependencies:

   pip install pandas numpy matplotlib seaborn scipy statsmodels scikit-learn gseapy

3. Launch and run dge_analysis.ipynb top to bottom. All figures and tables are written automatically to results/figures/, results/tables/, and results/enrichr/.

Skills Demonstrated

• Acquisition and wrangling of large public genomics datasets (TCGA / GDC / UCSC Xena)
• RNA-Seq count filtering and normalization
• High-throughput statistical hypothesis testing with proper multiple-testing correction (FDR)
• Exploratory data analysis: PCA, distribution diagnostics, hierarchical/heatmap visualization
• Biological interpretation of omics results via pathway/GO enrichment
• Building a clean, reproducible, portfolio-ready analysis pipeline

Limitations & Future Work

• The Xena star_counts file is distributed already in log2(count + 1) units; the notebook applies an additional log2(count + 1) transform on top of that during filtering. This double-transformation doesn't change the DEG calls qualitatively but means reported log2FC values aren't on a strictly "single log2 of raw counts" scale — worth revisiting if exact fold-change magnitudes matter downstream.
• The analysis subsamples 200 (100/100) of the available samples for tractability; scaling to the full cohort (1,118 Tumor / 137 Normal) would increase statistical power, especially for the underrepresented Normal group.
• GO enrichment did not return significant terms at the chosen cutoff — likely due to the small upregulated gene set (n=5); relaxing DEG thresholds, using a ranked GSEA approach instead of an overrepresentation test, or increasing sample size could yield more robust pathway-level conclusions.
• Per-gene Welch's t-tests are simpler than count-based models (e.g., DESeq2, edgeR) that explicitly model RNA-Seq's mean-variance relationship; adopting one of these could improve DEG calling accuracy.