Overview
Provide deterministic, copy-paste-ready code for multi-omics data processing tasks: ID mapping, batch correction, integration via MOFA2/mofapy2, enrichment analysis, and consensus clustering.
Usage
- Map gene identifiers across HGNC, Ensembl, and UniProt systems
- Run batch correction (ComBat/ComBat-seq) and multi-omic factor analysis (MOFA2/mofapy2)
- Perform GSEA enrichment and consensus clustering on integrated data
Core Concepts
Response Format
- Lead with the command or code the user needs — explain after
- Structure as: confirm inputs → working code → key parameters explained → gotchas
- One complete working example per task; do not show every alternative
- Keep code comments minimal and functional (what, not why-it-exists)
- Target: 50-100 lines of code with brief surrounding explanation
1. ID Mapping
1a. Python — HGNC Symbol ↔ Ensembl ↔ UniProt
import pandas as pd
import mygene
mg = mygene.MyGeneInfo()
# Input: list of HGNC symbols
genes = ["TP53", "BRCA1", "EGFR", "KRAS"]
# Query multiple ID types at once
result = mg.querymany(
genes,
scopes="symbol",
fields="ensembl.gene,uniprot.Swiss-Prot",
species="human",
returnall=True,
)
df = pd.DataFrame(result["out"])
# Handle many-to-many: explode nested Ensembl IDs
df_ensembl = df.explode("ensembl")
df_ensembl["ensembl_id"] = df_ensembl["ensembl"].apply(
lambda x: x["gene"] if isinstance(x, dict) else None
)
print(df_ensembl[["query", "ensembl_id"]].dropna())
1b. R — biomaRt
library(biomaRt)
ensembl <- useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")
# Map HGNC symbols to Ensembl gene IDs and UniProt accessions
mapping <- getBM(
attributes = c("hgnc_symbol", "ensembl_gene_id", "uniprotswissprot"),
filters = "hgnc_symbol",
values = c("TP53", "BRCA1", "EGFR", "KRAS"),
mart = ensembl
)
# Remove rows with empty UniProt
mapping <- mapping[mapping$uniprotswissprot != "", ]
print(mapping)
ID Mapping Parameter Table
| Parameter | Python (mygene) | R (biomaRt) |
|---|---|---|
| Input ID type | scopes="symbol" |
filters = "hgnc_symbol" |
| Output ID types | fields="ensembl.gene,uniprot.Swiss-Prot" |
attributes = c("ensembl_gene_id", "uniprotswissprot") |
| Species filter | species="human" |
dataset = "hsapiens_gene_ensembl" |
| Many-to-many handling | Explode nested lists | Multiple rows returned automatically |
2. Batch Correction
2a. Python — ComBat (pycombat)
from combat.pycombat import pycombat
import pandas as pd
import numpy as np
# expression_df: genes (rows) x samples (columns), log-transformed
# batch: list of batch labels, same order as columns
# covariates: DataFrame with biological covariates to protect
expression_df = pd.read_csv("expression_matrix.csv", index_col=0)
batch = [1, 1, 1, 2, 2, 2, 3, 3, 3] # batch labels per sample
# Protect biological covariate (e.g., disease status)
covar_df = pd.DataFrame({"disease": [0, 0, 1, 0, 1, 1, 0, 1, 1]})
corrected = pycombat(
expression_df,
batch,
mod=covar_df, # biological covariates to protect
par_prior=True, # parametric adjustment (default)
)
2b. R — ComBat / ComBat-seq
library(sva)
# For continuous data (log-transformed expression, protein intensity)
# expr_mat: genes x samples matrix
# batch: vector of batch labels
# mod: model matrix with biological covariates to protect
mod <- model.matrix(~ disease_status, data = sample_info)
corrected <- ComBat(
dat = expr_mat,
batch = sample_info$batch,
mod = mod,
par.prior = TRUE # parametric
)
# For count data (RNA-seq raw counts) — use ComBat-seq
corrected_counts <- ComBat_seq(
counts = count_mat,
batch = sample_info$batch,
group = sample_info$disease_status # biological variable to protect
)
Batch Correction Parameter Table
| Parameter | ComBat (continuous) | ComBat-seq (counts) |
|---|---|---|
| Input | Log-transformed matrix | Raw count matrix (integers) |
| Parametric prior | par.prior = TRUE (default) |
Not applicable |
| Biological covariate | mod (model matrix) |
group (vector) |
| Output | Corrected continuous matrix | Corrected count matrix |
| When to use | Proteomics, methylation M-values, log-expression | RNA-seq raw counts |
3. MOFA2 / mofapy2 Integration
3a. R — MOFA2
library(MOFA2)
# Prepare data: list of matrices, one per omic layer
# Each matrix: features x samples
data_list <- list(
transcriptomics = expr_mat, # genes x samples
proteomics = prot_mat, # proteins x samples
methylation = meth_mat # CpG sites x samples
)
# Create MOFA object
mofa <- create_mofa(data_list)
# Set data options
data_opts <- get_default_data_options(mofa)
# Set model options
model_opts <- get_default_model_options(mofa)
model_opts$num_factors <- 20 # start with 20; inactive factors pruned automatically
# Set training options
train_opts <- get_default_training_options(mofa)
train_opts$convergence_mode <- "slow" # more iterations for better convergence
train_opts$seed <- 42
# Prepare and run
mofa <- prepare_mofa(
mofa,
data_options = data_opts,
model_options = model_opts,
training_options = train_opts
)
mofa <- run_mofa(mofa, outfile = "mofa_model.hdf5")
# Inspect variance explained per factor per view
plot_variance_explained(mofa)
# Extract top features for factor 1
weights <- get_weights(mofa, views = "transcriptomics", factors = 1, as.data.frame = TRUE)
top_genes <- weights[order(abs(weights$value), decreasing = TRUE), ][1:50, ]
3b. Python — mofapy2
from mofapy2.run.entry_point import entry_point
import numpy as np
# Prepare data: list of numpy arrays [views][groups]
# Each array: samples x features
data = [
[expr_array], # view 0: transcriptomics, group 0
[prot_array], # view 1: proteomics, group 0
[meth_array], # view 2: methylation, group 0
]
ent = entry_point()
ent.set_data_options(scale_groups=False, scale_views=True)
ent.set_data_matrix(data, views_names=["transcriptomics", "proteomics", "methylation"])
ent.set_model_options(factors=20, spikeslab_weights=True, ard_factors=True)
ent.set_train_options(iter=1000, convergence_mode="slow", seed=42)
ent.build()
ent.run()
ent.save("mofa_model.hdf5")
MOFA Parameter Table
| Parameter | Default | Guidance |
|---|---|---|
num_factors / factors |
10 | Set 15–25 for exploratory analysis; inactive factors auto-pruned |
convergence_mode |
"fast" |
Use "slow" for final analysis; "fast" for exploration |
scale_views |
TRUE / True |
MUST be True to equalize variance across views |
spikeslab_weights |
TRUE / True |
Enables sparse feature selection per factor |
ard_factors |
TRUE / True |
Automatic relevance determination; prunes inactive factors |
seed |
None | ALWAYS set for reproducibility |
4. GSEA Enrichment
4a. Python — gseapy
import gseapy as gp
import pandas as pd
# ranked_genes: DataFrame with columns ["gene", "rank_metric"]
# rank_metric can be: log2FC, -log10(p) * sign(FC), or MOFA weight
ranked_genes = pd.read_csv("ranked_genes.csv")
# Pre-ranked GSEA against MSigDB Hallmark gene sets
result = gp.prerank(
rnk=ranked_genes,
gene_sets="MSigDB_Hallmark_2020",
min_size=15,
max_size=500,
permutation_num=1000,
seed=42,
outdir="gsea_results",
)
# Filter significant results
sig = result.res2d[result.res2d["FDR q-val"] < 0.05]
print(sig[["Term", "NES", "FDR q-val"]].sort_values("NES", ascending=False))
4b. R — fgsea with clusterProfiler
library(fgsea)
library(msigdbr)
# Get Hallmark gene sets for human
hallmark <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(hallmark$gene_symbol, hallmark$gs_name)
# ranked_genes: named numeric vector (names = gene symbols, values = rank metric)
ranked_genes <- setNames(de_results$log2FoldChange, de_results$gene_symbol)
ranked_genes <- sort(ranked_genes, decreasing = TRUE)
# Run fgsea
fgsea_res <- fgsea(
pathways = pathways,
stats = ranked_genes,
minSize = 15,
maxSize = 500,
nPermSimple = 10000
)
# Filter and sort
sig_pathways <- fgsea_res[fgsea_res$padj < 0.05, ]
sig_pathways <- sig_pathways[order(sig_pathways$NES, decreasing = TRUE), ]
print(sig_pathways[, c("pathway", "NES", "padj")])
Enrichment Parameter Table
| Parameter | gseapy (Python) | fgsea (R) |
|---|---|---|
| Gene set source | "MSigDB_Hallmark_2020" or GMT file |
msigdbr() + split() |
| Min gene set size | min_size=15 |
minSize = 15 |
| Max gene set size | max_size=500 |
maxSize = 500 |
| Permutations | permutation_num=1000 |
nPermSimple = 10000 |
| Significance threshold | FDR q-val < 0.05 |
padj < 0.05 |
| Rank metric | log2FC, signed -log10(p), or MOFA weight | Same |
5. Consensus Clustering
Python — scikit-learn based consensus
import numpy as np
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
def consensus_clustering(data, k_range=range(2, 9), n_iterations=100, seed=42):
"""Run consensus clustering over a range of k values.
data: samples x features array. Returns dict mapping k to results.
"""
rng = np.random.default_rng(seed)
n_samples = data.shape[0]
results = {}
for k in k_range:
co_occurrence = np.zeros((n_samples, n_samples))
co_sampled = np.zeros((n_samples, n_samples))
for _ in range(n_iterations):
# Subsample 80% of samples
idx = rng.choice(n_samples, size=int(0.8 * n_samples), replace=False)
km = KMeans(n_clusters=k, n_init=10, random_state=rng.integers(1e6))
labels = km.fit_predict(data[idx])
for i in range(len(idx)):
for j in range(i + 1, len(idx)):
co_sampled[idx[i], idx[j]] += 1
co_sampled[idx[j], idx[i]] += 1
if labels[i] == labels[j]:
co_occurrence[idx[i], idx[j]] += 1
co_occurrence[idx[j], idx[i]] += 1
# Consensus matrix
with np.errstate(divide="ignore", invalid="ignore"):
consensus = np.where(co_sampled > 0, co_occurrence / co_sampled, 0)
# Final clustering on consensus matrix
km_final = KMeans(n_clusters=k, n_init=10, random_state=seed)
final_labels = km_final.fit_predict(1 - consensus) # distance = 1 - consensus
sil = silhouette_score(1 - consensus, final_labels, metric="precomputed")
results[k] = {"consensus_matrix": consensus, "silhouette": sil, "labels": final_labels}
print(f"k={k}: silhouette={sil:.3f}")
return results
# Usage with MOFA factors
# factors: samples x factors array from MOFA
results = consensus_clustering(factors, k_range=range(2, 9))
# Select optimal k by maximum silhouette score
optimal_k = max(results, key=lambda k: results[k]["silhouette"])
print(f"Optimal k: {optimal_k}")
R — ConsensusClusterPlus
library(ConsensusClusterPlus)
# input_mat: features x samples
results <- ConsensusClusterPlus(
d = input_mat,
maxK = 8,
reps = 1000,
pItem = 0.8, # proportion of samples to subsample
pFeature = 1, # use all features
clusterAlg = "km", # k-means
distance = "euclidean",
seed = 42,
plot = "pdf",
title = "consensus_clustering"
)
# Inspect consensus CDF and delta area plots in output PDF
# Select k where CDF curve is flattest (highest consensus)
Consensus Clustering Parameter Table
| Parameter | Python (custom) | R (ConsensusClusterPlus) |
|---|---|---|
| k range | k_range=range(2, 9) |
maxK = 8 |
| Iterations | n_iterations=100 |
reps = 1000 |
| Subsample fraction | 0.8 (hardcoded) | pItem = 0.8 |
| Algorithm | KMeans | clusterAlg = "km" |
| Selection criterion | Silhouette score | CDF delta area plot |
Common Mistakes
-
Wrong: Mapping gene IDs using symbol-based joins without handling duplicates or deprecated symbols Right: Use HGNC ID or Ensembl stable ID as the primary key; resolve ambiguous symbols via the HGNC multi-symbol checker Why: Gene symbols are not unique (e.g., "MARCH1" was renamed to "MARCHF1") — symbol joins silently drop or duplicate rows
-
Wrong: Running ComBat without including biological covariates in the model matrix Right: Always pass
mod = model.matrix(~ disease_status + sex, data=meta)to protect biological signal Why: Without protected covariates, ComBat removes biological variation that correlates with batch -
Wrong: Using
mofapy2with raw counts instead of normalized, variance-stabilized data Right: Apply log2(CPM+1) or VST to RNA-seq counts before passing to MOFA; use log2 intensity for proteomics Why: MOFA assumes Gaussian-distributed continuous data — raw counts violate this and produce meaningless factors -
Wrong: Running gseapy/fgsea with a gene list ranked by fold change alone (no significance weighting) Right: Rank by
-log10(p) × sign(log2FC)or use the Wald statistic from DESeq2/limma Why: Fold-change-only ranking treats noisy large-FC genes the same as significant ones, inflating false enrichments -
Wrong: Using default
min_size=1in fgsea, allowing tiny gene sets to appear significant Right: SetminSize=15, maxSize=500to filter gene sets to biologically interpretable sizes Why: Very small gene sets produce unstable enrichment scores; very large sets are too generic to be informative