Cancer Multiomics

Cancer Multi-Omics Analysis

Integrated analysis of expression, mutation, copy number, and methylation data from TCGA and GEO for solid tumor characterization.

When to Use This Skill

Activate when the user requests:

• TCGA data download and processing for any cancer type
• Differential expression analysis between tumor conditions
• Pathway and gene set enrichment analysis
• Mutation landscape analysis (oncoplot, signatures, co-occurrence)
• Copy number variation analysis
• DNA methylation analysis (450K/EPIC arrays)
• Integration of two or more omics layers

Inputs

| Data Type | Format | Source | |-----------|--------|--------| | Expression | Raw counts (STAR - Counts) | TCGA via TCGAbiolinks (v2.38+), GEO via GEOquery | | Mutations | MAF | TCGA GDC, cBioPortal | | Copy number | Segment files, GISTIC2.0 output | TCGA GDC, GDAC Firehose | | Methylation | IDAT files or beta-value matrices | TCGA GDC, GEO | | Clinical | Tabular | TCGA GDC, cBioPortal |

---

Expression Analysis

Data Retrieval

library(TCGAbiolinks)  # v2.38.0+, Bioconductor 3.22

query  DE analysis: feed directly to DESeq2 (handles normalization internally)
    -> Visualization (PCA, heatmap): apply vst() or rlog() from DESeq2
    -> Cross-sample comparison: vst() preferred for n > 30 samples (faster than rlog)
  TPM/FPKM (already normalized)
    -> DE analysis: STOP. Go back, get raw counts. TPM/FPKM invalid for DESeq2/edgeR.
    -> Correlation/visualization: log2(TPM + 1), acceptable
    -> Gene set scoring (ssGSEA, GSVA): TPM acceptable as input
  RSEM expected counts (non-integer)
    -> round() before DESeq2: DESeqDataSetFromMatrix(countData = round(counts), ...)

Differential Expression (DESeq2 v1.50+)

library(DESeq2)  # v1.50.2+, Bioconductor 3.22
set.seed(42)

counts_mat = 10 counts in at least 5 samples
# More stringent than rowSums >= 10; reduces false discoveries
keep = 10) >= 5
dds  1)
sig 1): %d up, %d down",
  sum(sig$log2FoldChange > 0), sum(sig$log2FoldChange  Include batch in DESeq2 design formula. Do NOT pre-correct.
  For PCA/heatmap (log-scale data) -> limma::removeBatchEffect on vst/rlog values
  For downstream tools requiring counts -> ComBat-seq (sva package) on raw counts
  NEVER: ComBat (original) on raw counts. NEVER: removeBatchEffect before DE.

# Visualization with batch removal
vsd  1]

ego  1,
  ifelse(df$log2FoldChange > 0, "Up", "Down"), "NS")
top10 1
     10000 DEGs: check for confounders, tumor purity bias
  PCA: tumor and normal should separate on PC1 or PC2
  Known markers: TP53, EGFR, KRAS mutations reflected in expression?
  Direction: oncogenes should trend up in tumor, TSGs down
  Volcano symmetry: roughly symmetric unless strong biological bias expected

Mutation Analysis

Data Retrieval

library(TCGAbiolinks)
library(maftools)  # v2.22+, Bioconductor 3.22

# Method 1: GDCquery_Maf (downloads open-access MAF, hg38-aligned)
# Pipelines: mutect2, muse, varscan2, somaticsniper
maf_df  Most stringent, used in TCGA PanCancer Atlas
    Access via: maftools::tcgaLoad(study = "LUAD")
  Single caller:
    MuTect2 -> Best sensitivity for low-frequency variants
    MuSE -> Good for paired tumor-normal, models tumor heterogeneity
    VarScan2 -> Works with and without matched normal
    SomaticSniper -> Bayesian, good specificity
  For most analyses -> MuTect2 or MC3 consensus

Visualization

plotmafSummary(maf, rmOutlier = TRUE, addStat = "median")

oncoplot(maf, top = 20,
  clinicalFeatures = c("ajcc_pathologic_stage", "gender"),
  sortByAnnotation = TRUE)

# Lollipop: protein domain plot for a single gene
lollipopPlot(maf, gene = "TP53", AACol = "HGVSp_Short",
  showMutationRate = TRUE)

# Rainfall: inter-mutation distance, detects kataegis (localized hypermutation)
rainfallPlot(maf, detectChangePoints = TRUE, pointSize = 0.4)

Tumor Mutation Burden

# TMB = total nonsynonymous mutations / capture region size (MB)
laml.tmb = 30 samples for reliable decomposition
tnm  0.85 is a confident assignment

plotSignatures(sig, contributions = FALSE)  # signature profiles
plotSignatures(sig, contributions = TRUE)   # per-sample contributions

Driver Gene Detection

# oncodrive: identifies genes with clustered mutations (positional clustering)
# Not dN/dS — uses the concept that driver mutations cluster in specific protein regions
oncodrive_res 0 = gain,  0.3   -> Gain
  > 0.7   -> Amplification (high-level gain)
   Loss (heterozygous deletion)
   Deep deletion (likely homozygous)
  These are conventions, not absolute cutoffs. Adjust based on tumor purity.

GISTIC2.0 Analysis

# GISTIC identifies recurrent focal CNV events across a cohort
# Pre-computed GISTIC results available from Firehose/GDAC for most TCGA cancers:
# https://gdac.broadinstitute.org/

# Read GISTIC output with maftools
gistic  0.01
keep_samples  0.01 in > 5% of samples)
keep 
  (ncol(grSet) * 0.95)
grSet  ~420-440K retained
# EPIC: ~866K -> ~780-810K retained

Alternative: ChAMP Pipeline

library(ChAMP)

# All-in-one loading with built-in filtering
# Handles detection p, SNPs, cross-reactive, sex chromosomes in one call
myLoad  Fine-grained control, Funnorm (best for cancer), EPIC v2 support
  ChAMP -> Rapid standardized pipeline, built-in filtering + batch correction
  Both use limma internally for DMP calling

Beta-Values vs M-Values

Beta-values: ratio = M / (M + U + 100), range [0, 1]
  Use for: biological interpretation, visualization, reporting delta-beta
  Problem: heteroscedastic (variance compressed at 0 and 1)

M-values: log2(M / U), range (-Inf, Inf)
  Use for: ALL statistical testing (limma, t-tests, linear models)
  Reason: approximately normal, satisfies linear model assumptions

Conversion: M = log2(beta / (1 - beta))
Always test on M-values, report delta-beta for biological effect size.

Differentially Methylated Positions (DMPs)

library(limma)

# Extract M-values for statistical testing
mVals  0.2, ]
message(sprintf("Significant DMPs: %d hyper, %d hypo",
  sum(sig_dmps$deltaBeta > 0), sum(sig_dmps$deltaBeta  BRAF V600E, MLH1 hypermethylation, MSI-H
    CIMP-L -> KRAS mutations
  Glioblastoma (GBM): G-CIMP defined by IDH1 mutation + global hypermethylation
    G-CIMP -> proneural subtype, younger patients, better survival
  Pan-cancer: 19/26 TCGA cancer types have confirmed CIMP subtypes

Integration: Methylation-Expression Correlation

# --- Match samples across platforms ---
met_patients  0.2) |
| `dmr_ranges.rds` | RDS | GRanges of differentially methylated regions |
| `grSet.rds` | RDS | Filtered, normalized GenomicRatioSet |
| `cimp_clusters.csv` | CSV | Sample-to-cluster assignments from consensus clustering |
| `meth_expr_cor.csv` | CSV | Per-gene Spearman rho and adjusted p-values |
| `dmp_heatmap.pdf` | PDF | Top DMPs across samples |

### Validation Checks

After running methylation analysis, verify: DMP count: Expect 10,000-100,000 DMPs for tumor vs normal at adj.P.Val 0.7 at many probes: check for sample swaps or contamination Known markers: MGMT promoter methylation in GBM (predicts temozolomide response) MLH1 methylation in MSI-H colorectal (silences mismatch repair) BRCA1 methylation in basal-like breast cancer Probe filtering: verify 40,000-65,000 probes removed (450K) or 55,000-85,000 (EPIC) Integration: promoter methylation-expression correlation should be predominantly negative Median rho around -0.05 to -0.15 across all genes Gene body methylation should show weak positive correlation

## Common Pitfalls

### Expression
1. **DESeq2 input**: Using FPKM/TPM as input. DESeq2 requires raw integer counts.
2. **Pre-filtering GSEA input**: GSEA requires the full ranked gene list. Never pre-filter to significant genes only (that is ORA, not GSEA).
3. **LFC shrinkage for hypothesis testing**: apeglm shrinkage is for ranking and visualization. Use unshrunken results for strict significance calls unless using `lfcThreshold > 0`.
4. **Gene ID mismatch**: TCGA uses Ensembl IDs with version suffixes (ENSG00000141510.18). Strip versions before mapping: `sub("\\..*", "", ensembl_ids)`.

### Mutation
5. **Mixing callers**: Different variant callers produce different mutation sets. Pick one pipeline (or use MC3 consensus) and stick with it across the cohort.
6. **Signature sample size**: Extracting signatures from 0.3 gain,  60%) or include as covariate.
18. **TCGA barcode matching**: When integrating data types, match on patient barcode (first 12 characters), not full barcode. A patient may have multiple aliquots.

## Related Skills

- [`immune-deconvolution`](../immune-deconvolution/SKILL.md): Estimate immune cell composition from the expression data produced here
- [`survival-analysis`](../survival-analysis/SKILL.md): Build prognostic models using mutation, expression, and methylation features as Cox covariates

## Public Datasets for Testing

| Dataset | Samples | Use Case |
|---------|---------|----------|
| TCGA-LUAD | 585 | Lung adenocarcinoma, well-characterized drivers (EGFR, KRAS, ALK) |
| TCGA-BRCA | 1098 | Breast cancer, strong molecular subtypes (Luminal A/B, HER2, Basal) |
| TCGA-GBM | 617 | Glioblastoma, classic for survival and methylation subtyping |
| TCGA-COAD | 521 | Colorectal, MSI-H vs MSS comparison |

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [zamushwani2](https://github.com/zamushwani2)
- **Source:** [zamushwani2/biomedical-ai-skills](https://github.com/zamushwani2/biomedical-ai-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.