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

Melanoma Subtype Marker Genes

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).

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.

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.

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).

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.

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"])

 # 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