# Bio Applied Cancer Transcriptomics

> Classify tumor RNA-seq into subtypes (melanoma Tirosh/Harbst on TCGA-SKCM): log1p/z-score, PCA/t-SNE, hierarchical clustering, random forest, Kaplan-Meier survival. Use when subtyping cBioPortal expression data.

- Skill: `pavel-kravchenko/bio-applied-cancer-transcriptomics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-cancer-transcriptomics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-cancer-transcriptomics/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-cancer-transcriptomics

---


# Cancer Transcriptomics: Subtype Classification

## When to Use
- Classifying tumor samples into molecular subtypes from a bulk RNA-seq expression matrix (e.g., melanoma Tirosh 3-class or Harbst 4-class from TCGA-SKCM)
- Exploring cohort structure with PCA/t-SNE/clustermaps before deciding on a subtyping scheme
- Training a marker-gene-based classifier (random forest) or unsupervised clustering (hierarchical) to assign subtype labels
- Testing whether a molecular subtype is associated with patient outcome via Kaplan-Meier / log-rank
- Comparing two independent subtype classifications for concordance (e.g., Tirosh vs Harbst)

## Version Compatibility
- Python ≥3.10, pandas ≥2.0, scikit-learn ≥1.4, seaborn ≥0.13, lifelines ≥0.28
- R ≥4.3, survival ≥3.5, survminer ≥0.4.9 (for KM plotting in R)

## Prerequisites
- `pip install pandas scikit-learn seaborn lifelines matplotlib`
- Expression matrix in genes × samples format (cBioPortal `data_mrna_seq_v2_rsem.txt` style); familiarity with PCA, hierarchical clustering, and censored survival data
- Related: `bio-differential-expression-deseq2-basics` for upstream normalization, `cbioportal-database` for pulling TCGA cohorts

## Key References
- Tirosh et al. (2016) *Science* — single-cell dissection of melanoma intratumoral heterogeneity
- Harbst et al. (2016) *Clin Cancer Res* — four molecular subtypes from bulk RNA-seq (TCGA-SKCM)
- Data: [cBioPortal TCGA-SKCM](https://www.cbioportal.org/study/summary?id=skcm_tcga)

## Melanoma Subtype Marker Genes

| Subtype | Key Markers | Classification |
|---------|-------------|----------------|
| Pigmentation (MITF-high) | MITF, DCT, TYRP1, MLANA | Tirosh + Harbst |
| Keratin | KRT5, KRT14, KRT6A, EGFR | Tirosh |
| Immune | CD3D, CD8A, GZMB, PRF1 | Tirosh + Harbst |
| Proliferative | MKI67, TOP2A, CDK1 | Harbst |
| Normal-like | VIM, CDH2, FN1 | Harbst |

## Preprocessing

**Goal:** turn a raw RSEM/counts matrix into a variance-stabilized, gene-scaled matrix suitable for clustering and classification.
**Approach:** log1p to compress dynamic range, z-score per gene across samples, then filter to the most variable genes (compute variance on log1p, not z-scored, data to avoid circularity).

```python
import numpy as np
import pandas as pd


def preprocess_expression(expr_df: pd.DataFrame, n_top_genes: int = 1500) -> pd.DataFrame:
    """Log-transform, z-score, and variance-filter a genes x samples matrix.

    Args:
        expr_df: raw expression matrix, genes as rows, samples as columns.
        n_top_genes: number of highest-variance genes to retain.

    Returns:
        z-scored matrix restricted to the top-variance genes.
    """
    expr_log = np.log1p(expr_df)
    expr_z = expr_log.subtract(expr_log.mean(axis=1), axis=0).divide(
        expr_log.std(axis=1) + 1e-8, axis=0
    )
    top_genes = expr_log.var(axis=1).nlargest(n_top_genes).index  # variance on log1p, not z-scored
    return expr_z.loc[top_genes]


expr_df = pd.read_csv("data_mrna_seq_v2_rsem.txt", sep="\t", index_col=0)
expr_top = preprocess_expression(expr_df)  # shape: (1500, n_samples)
```

## Exploratory Analysis (PCA / t-SNE / Clustermap)

**Goal:** check for batch effects and subtype-driven structure before modeling.
**Approach:** PCA for a linear overview, t-SNE (seeded from PCA for large cohorts) for local structure, and a clustermap of top variable genes colored by known subtype.

```python
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import seaborn as sns

X = expr_top.values.T  # (n_samples, n_genes)

pca = PCA(n_components=2, random_state=42)
X_pca = pca.fit_transform(X)
print(f"PC1+PC2 variance explained: {pca.explained_variance_ratio_.sum():.1%}")

tsne = TSNE(n_components=2, perplexity=30, random_state=42)
X_tsne = tsne.fit_transform(X_pca[:, :50] if X.shape[0] > 500 else X)

subtype_palette = {"MITF-low": "#4C72B0", "Keratin": "#DD8452", "Immune": "#55A868"}
col_colors = pd.Series(subtype_labels, index=expr_top.columns).map(subtype_palette)
sns.clustermap(expr_top.iloc[:50], col_colors=col_colors, cmap="RdBu_r",
               center=0, vmin=-3, vmax=3, xticklabels=False, figsize=(12, 8))
```

## Subtype Classification

**Goal:** assign a subtype label to each sample either via a supervised classifier (when ground-truth labels exist, e.g. Tirosh 3-class) or via unsupervised marker-based clustering (e.g. Harbst 4-class).
**Approach:** random forest on the top-variance genes for the supervised case; Ward-linkage hierarchical clustering on a small curated marker panel for the unsupervised case, with clusters relabeled by centroid marker expression.

```python
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
from sklearn.cluster import AgglomerativeClustering

 # Tirosh 3-class: supervised random forest
X = expr_top.values.T
X_train, X_test, y_train, y_test = train_test_split(
    X, subtype_labels, test_size=0.30, random_state=42, stratify=subtype_labels
)
rf = RandomForestClassifier(n_estimators=200, random_state=42, n_jobs=-1)
rf.fit(X_train, y_train)
print(classification_report(y_test, rf.predict(X_test)))
importance = pd.Series(rf.feature_importances_, index=expr_top.index).nlargest(20)

 # Harbst 4-class: unsupervised clustering on curated markers
harbst_markers = ["MITF", "DCT", "TYRP1", "MLANA", "MKI67", "TOP2A", "CDK1",
                   "VIM", "CDH2", "FN1", "CD3D", "CD8A", "GZMB", "PRF1"]
X_harbst = expr_z.loc[harbst_markers].values.T
cluster_ids = AgglomerativeClustering(n_clusters=4, linkage="ward").fit_predict(X_harbst)

centroids = pd.DataFrame(X_harbst, columns=harbst_markers)
centroids["cluster"] = cluster_ids
ctrs = centroids.groupby("cluster").mean()
subtype_map = {
    ctrs[["MITF", "DCT", "TYRP1", "MLANA"]].mean(axis=1).idxmax(): "Pigmentation",
    ctrs[["MKI67", "TOP2A", "CDK1"]].mean(axis=1).idxmax(): "Proliferative",
    ctrs[["CD3D", "CD8A", "GZMB", "PRF1"]].mean(axis=1).idxmax(): "High-immune",
}
remaining = [c for c in range(4) if c not in subtype_map]
if remaining:
    subtype_map[remaining[0]] = "Normal-like"
harbst_labels = np.array([subtype_map.get(c, f"Cluster_{c}") for c in cluster_ids])
```

## Subtype Concordance

**Goal:** quantify agreement between two independently derived subtype schemes (e.g. Tirosh vs Harbst) on the same cohort.
**Approach:** cross-tabulate labels and summarize agreement with normalized mutual information (0 = independent, 1 = perfect agreement).

```python
from sklearn.metrics import normalized_mutual_info_score

crosstab = pd.crosstab(pd.Series(tirosh_labels, name="Tirosh"),
                        pd.Series(harbst_labels, name="Harbst"))
nmi = normalized_mutual_info_score(tirosh_labels, harbst_labels)
print(crosstab)
print(f"NMI = {nmi:.4f}")
```

## Survival Analysis

**Goal:** test whether a molecular subtype stratifies patient outcome.
**Approach:** Kaplan-Meier curves per subtype with `events` as a boolean censoring indicator, compared via log-rank test. Use `lifelines` in Python if available; fall back to a manual estimator otherwise.

```python
import matplotlib.pyplot as plt


def km_estimate(times: np.ndarray, events: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Manual Kaplan-Meier estimator, used when lifelines is unavailable.

    Args:
        times: observed time-to-event or censoring, per sample.
        events: boolean array, True = event observed, False = censored.

    Returns:
        (time_points, survival_probability) step-function arrays.
    """
    order = np.argsort(times)
    t_sorted, e_sorted = times[order], events[order]
    unique_event_times = np.unique(t_sorted[e_sorted])
    surv, t_out = [1.0], [0.0]
    for t in unique_event_times:
        deaths = np.sum((t_sorted == t) & e_sorted)
        at_risk = np.sum(t_sorted >= t)
        surv.append(surv[-1] * (1 - deaths / at_risk))
        t_out.append(t)
    return np.array(t_out), np.array(surv)


try:
    from lifelines import KaplanMeierFitter
    from lifelines.statistics import logrank_test

    fig, ax = plt.subplots(figsize=(9, 5))
    for subtype in ["MITF-low", "Keratin", "Immune"]:
        mask = subtype_labels == subtype
        kmf = KaplanMeierFitter()
        kmf.fit(survival_times[mask], events[mask], label=subtype)
        kmf.plot_survival_function(ax=ax, ci_show=True)
    res = logrank_test(t1, t2, e1, e2)
    print(f"p = {res.p_value:.4f}")
except ImportError:
    t_km, s_km = km_estimate(survival_times[subtype_labels == "Immune"],
                              events[subtype_labels == "Immune"])
```

```r
 # Equivalent Kaplan-Meier + log-rank in R
library(survival)
library(survminer)

fit <- survfit(Surv(survival_time, event) ~ subtype, data = clinical_df)
survdiff(Surv(survival_time, event) ~ subtype, data = clinical_df)  # log-rank test
ggsurvplot(fit, data = clinical_df, pval = TRUE, conf.int = TRUE,
           risk.table = TRUE, palette = c("#4C72B0", "#DD8452", "#55A868"))
```

## Pitfalls

- **log1p before z-score, not after**: log1p compresses dynamic range first; z-scoring after ensures comparable gene scales; reversing the order gives distorted variance estimates
- **Variance filter on z-scored data is circular**: compute gene variance on log1p (not z-scored) data for filtering, then z-score the filtered matrix
- **t-SNE is not deterministic**: always set `random_state`; t-SNE axes have no biological meaning — distances between clusters are not interpretable, only local neighborhood structure is
- **Ward linkage requires Euclidean distance**: `AgglomerativeClustering(linkage='ward')` assumes Euclidean; for correlation-based clustering use `linkage='average'` with a precomputed distance matrix
- **Kaplan-Meier requires a censoring indicator**: `events` must be boolean (True = event observed); patients lost to follow-up are censored (False), not dead — mislabeling censored as events inflates mortality estimates
- **Batch effects in TCGA**: TCGA-SKCM has plate/batch effects; always check PC1/PC2 coloring by batch variable before biological interpretation; use ComBat or limma's `removeBatchEffect` if needed

## See Also
- `bio-differential-expression-deseq2-basics` — upstream DE testing before subtyping
- `bio-machine-learning-survival-analysis` — Cox regression and more general survival modeling
- `bio-data-visualization-heatmaps-clustering` — clustermap/heatmap styling details
- `cbioportal-database` — pulling TCGA-SKCM expression and clinical data