multi-omics-pipeline

多组学流程 — 流程技能

概述

为多组学数据处理任务提供确定性、可直接复制粘贴的代码：ID 映射、批次校正、通过 MOFA2/mofapy2 进行整合、富集分析和共识聚类。

用法

• 在 HGNC、Ensembl 和 UniProt 系统之间映射基因标识符
• 运行批次校正（ComBat/ComBat-seq）和多组学因子分析（MOFA2/mofapy2）
• 对整合数据执行 GSEA 富集分析和共识聚类

核心概念

---

响应格式

• 先给出用户需要的命令或代码，再进行说明
• 按以下结构组织：确认输入 → 可运行的代码 → 关键参数说明 → 注意事项
• 每项任务提供一个完整可运行的示例；不要展示所有替代方案
• 代码注释保持简洁且实用（说明用途，不解释存在原因）
• 目标：50–100 行代码，并附简短说明

1. ID 映射

1a. Python — HGNC 符号 ↔ 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 映射参数表

参数	Python (mygene)	R (biomaRt)
输入 ID 类型	scopes="symbol"	filters = "hgnc_symbol"
输出 ID 类型	fields="ensembl.gene,uniprot.Swiss-Prot"	attributes = c("ensembl_gene_id", "uniprotswissprot")
物种筛选	species="human"	dataset = "hsapiens_gene_ensembl"
多对多处理	展开嵌套列表	自动返回多行

---

2. 批次校正

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
)

批次校正参数表

参数	ComBat（连续数据）	ComBat-seq（计数数据）
输入	对数转换后的矩阵	原始计数矩阵（整数）
参数先验	par.prior = TRUE（默认）	不适用
生物学协变量	mod（模型矩阵）	group（向量）
输出	校正后的连续矩阵	校正后的计数矩阵
适用场景	蛋白质组学、甲基化 M 值、对数表达量	RNA-seq 原始计数

---

3. MOFA2 / mofapy2 整合

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 参数表

参数	默认值	指导建议
num_factors / factors	10	探索性分析可设为 15–25；非活跃因子会自动剪除
convergence_mode	"fast"	最终分析使用 "slow"；探索时使用 "fast"
scale_views	TRUE / True	必须为 True，以均衡各视图间的方差
spikeslab_weights	TRUE / True	启用每个因子的稀疏特征选择
ard_factors	TRUE / True	自动相关性判定；剪除非活跃因子
seed	None	始终设置，以确保可复现性

---

4. GSEA 富集分析

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

富集分析参数表

参数	gseapy（Python）	fgsea（R）
基因集来源	"MSigDB_Hallmark_2020" 或 GMT 文件	msigdbr() + split()
基因集最小大小	min_size=15	minSize = 15
基因集最大大小	max_size=500	maxSize = 500
置换次数	permutation_num=1000	nPermSimple = 10000
显著性阈值	FDR q-val < 0.05	padj < 0.05
排序指标	log2FC、有符号的 -log10(p) 或 MOFA 权重	相同

---

5. 共识聚类

Python — 基于 scikit-learn 的共识聚类

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)

共识聚类参数表

参数	Python（自定义）	R (ConsensusClusterPlus)
k 范围	k_range=range(2, 9)	maxK = 8
迭代次数	n_iterations=100	reps = 1000
子采样比例	0.8（硬编码）	pItem = 0.8
算法	KMeans	clusterAlg = "km"
选择标准	轮廓系数	CDF 曲线下面积图

---

常见错误

• 错误： 使用基于基因符号的连接映射基因 ID，却未处理重复项或已弃用的符号 正确做法： 使用 HGNC ID 或 Ensembl 稳定 ID 作为主键；通过 HGNC 多符号检查器解决有歧义的符号 原因： 基因符号并非唯一（例如，“MARCH1”已更名为“MARCHF1”）——基于符号的连接会在不提示的情况下丢弃或重复行

• 错误： 在模型矩阵中未包含生物学协变量就运行 ComBat 正确做法： 始终传入 mod = model.matrix(~ disease_status + sex, data=meta)，以保护生物学信号 原因： 如果没有受保护的协变量，ComBat 会移除与批次相关的生物学变异

• 错误： 使用 mofapy2 处理原始计数，而非归一化且经过方差稳定变换的数据 正确做法： 将 RNA-seq 计数转换为 log2(CPM+1) 或 VST 后再传入 MOFA；蛋白质组学数据使用 log2 强度值 原因： MOFA 假设数据是服从高斯分布的连续数据——原始计数不符合这一假设，会产生毫无意义的因子

• 错误： 使用仅按倍数变化排序的基因列表运行 gseapy/fgsea（未考虑显著性权重） 正确做法： 按 -log10(p) × sign(log2FC) 排序，或使用 DESeq2/limma 的 Wald 统计量 原因： 仅按倍数变化排序会将噪声较大的高 FC 基因与显著基因同等对待，导致假性富集增加

• 错误： 在 fgsea 中使用默认的 min_size=1，导致极小的基因集也可能显著 正确做法： 设置 minSize=15, maxSize=500，将基因集筛选到具有生物学可解释性的大小 原因： 极小的基因集会产生不稳定的富集分数；极大的基因集则过于宽泛，无法提供有用信息