multi-omics-pipeline

マルチオミクスパイプライン — パイプラインスキル

概要

マルチオミクスデータ処理タスク向けに、決定論的でコピー＆ペーストしてすぐに使えるコードを提供します。ID マッピング、バッチ補正、MOFA2/mofapy2 による統合、エンリッチメント解析、コンセンサスクラスタリングを扱います。

使用方法

• HGNC、Ensembl、UniProt の各システム間で遺伝子識別子をマッピングします
• バッチ補正（ComBat/ComBat-seq）とマルチオミクス因子分析（MOFA2/mofapy2）を実行します
• 統合データに対して GSEA エンリッチメント解析とコンセンサスクラスタリングを実施します

主な概念

---

応答形式

• 必要なコマンドまたはコードを最初に示し、その後で説明します
• 次の構成にします：入力の確認 → 動作するコード → 主なパラメーターの説明 → 注意点
• 各タスクにつき、動作する完全な例を1つ示します。すべての代替案は示しません
• コードコメントは最小限にし、実用的な内容にします（存在理由ではなく、何をするか）
• 目安：コード50～100行と、簡潔な補足説明

1. ID マッピング

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 マッピングのパラメータ表

パラメーター	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	なし	再現性を確保するため、必ず設定してください

---

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 — clusterProfiler を使った fgsea

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 を実行する 正しい方法： MOFA に渡す前に RNA-seq カウント値に log2(CPM+1) または VST を適用し、プロテオミクスには log2 強度を使用します 理由： MOFA はガウス分布に従う連続データを前提としています。生カウント値はこの前提に反するため、意味のない因子が生成されます。

• 誤り： 有意性による重み付けなしで、フォールドチェンジのみで順位付けした遺伝子リストを使って gseapy/fgsea を実行する 正しい方法： -log10(p) × sign(log2FC) で順位付けするか、DESeq2/limma の Wald 統計量を使用します 理由： フォールドチェンジのみの順位付けでは、ノイズの多い大きな FC の遺伝子が有意な遺伝子と同じように扱われ、偽のエンリッチメントが増加します。

• 誤り： fgsea でデフォルトの min_size=1 を使用し、小さな遺伝子セットが有意になるのを許容する 正しい方法： minSize=15, maxSize=500 を設定し、生物学的に解釈可能なサイズの遺伝子セットに絞り込みます 理由： 非常に小さな遺伝子セットでは不安定なエンリッチメントスコアが生じ、非常に大きなセットは一般的すぎて有益な情報が得られません。