| name | scientific-metabolomics |
| description | メタボロミクス解析スキル。Pareto スケーリング、PLS-DA + VIP スコア、置換検定（Q²）、
代謝パスウェイ濃縮解析（Fisher exact test）、代謝物相関ネットワーク、
Volcano プロット/箱ひげ図による差次代謝物同定パイプライン。
Scientific Skills Exp-07 で確立したパターン。
|
| tu_tools | [{"key":"hmdb","name":"HMDB","description":"ヒトメタボロームデータベース"},{"key":"metabolomics_workbench","name":"Metabolomics Workbench","description":"Metabolomics Workbench REST API メタボロームデータ・RefMet"}] |
Scientific Metabolomics Analysis
LC-MS / GC-MS / NMR ベースのメタボロミクスデータを対象に、品質管理→前処理→
単変量解析→多変量解析→パスウェイ解析の標準パイプラインを提供する。
メタボロミクス固有の統計手法（PLS-DA、VIP スコア）に特化している。
When to Use
- メタボロミクスデータの統計解析パイプラインが必要なとき
- PLS-DA による群間判別＋VIP スコアによるバイオマーカー候補を算出するとき
- 代謝パスウェイ濃縮解析が必要なとき
- Pareto スケーリングや代謝物相関ネットワークが必要なとき
Quick Start
1. メタボロミクス前処理
import numpy as np
import pandas as pd
from sklearn.impute import KNNImputer
def metabolomics_preprocessing(df, sample_col="Sample_ID", group_col="Group",
min_detect_pct=0.5):
"""
メタボロミクス標準前処理パイプライン。
1. 低検出率代謝物の除去
2. KNN 欠損値補完
3. log2 変換
4. Pareto スケーリング
"""
metabolite_cols = [c for c in df.columns if c not in [sample_col, group_col]]
detect_rate = df[metabolite_cols].notna().mean()
keep = detect_rate[detect_rate >= min_detect_pct].index.tolist()
removed = len(metabolite_cols) - len(keep)
print(f" Removed {removed} metabolites with <{min_detect_pct*100:.0f}% detection rate")
metabolite_cols = keep
imputer = KNNImputer(n_neighbors=5)
df[metabolite_cols] = imputer.fit_transform(df[metabolite_cols])
df[metabolite_cols] = np.log2(df[metabolite_cols].clip(lower=1e-10) + 1)
for col in metabolite_cols:
mean = df[col].mean()
std = df[col].std()
df[col] = (df[col] - mean) / np.sqrt(std + 1e-10)
return df, metabolite_cols
2. 単変量解析 — 差次代謝物同定
from scipy.stats import mannwhitneyu, ttest_ind
from statsmodels.stats.multitest import multipletests
def univariate_analysis(df, metabolite_cols, group_col, group1, group2,
test="mannwhitneyu", correction="fdr_bh"):
"""
2 群間の差次代謝物を同定する。
Returns:
DataFrame with columns: metabolite, log2FC, pvalue, padj, significant
"""
g1 = df[df[group_col] == group1]
g2 = df[df[group_col] == group2]
results = []
for met in metabolite_cols:
v1 = g1[met].dropna()
v2 = g2[met].dropna()
if test == "mannwhitneyu":
stat, pval = mannwhitneyu(v1, v2, alternative="two-sided")
else:
stat, pval = ttest_ind(v1, v2)
log2fc = v2.mean() - v1.mean()
results.append({"metabolite": met, "log2FC": log2fc, "pvalue": pval})
results_df = pd.DataFrame(results)
reject, padj, _, _ = multipletests(results_df["pvalue"], method=correction)
results_df["padj"] = padj
results_df["significant"] = reject
results_df["neg_log10p"] = -np.log10(results_df["pvalue"] + 1e-300)
results_df = results_df.sort_values("pvalue")
return results_df
3. PLS-DA + VIP スコア
from sklearn.cross_decomposition import PLSRegression
from sklearn.preprocessing import LabelEncoder
def plsda_analysis(X, y, n_components=2):
"""
PLS-DA を実行し、VIP スコアを算出する。
VIP (Variable Importance in Projection):
VIP_j = sqrt(p * Σ(q²_a * w²_ja) / Σ(q²_a))
VIP > 1 の変数がバイオマーカー候補
Parameters:
X: 代謝物データ行列 (n_samples, n_metabolites)
y: グループラベル
Returns:
pls_model, scores, vip_scores
"""
le = LabelEncoder()
y_encoded = le.fit_transform(y).astype(float)
pls = PLSRegression(n_components=n_components, scale=True)
pls.fit(X, y_encoded)
scores = pls.transform(X)
T = pls.x_scores_
W = pls.x_weights_
Q = pls.y_loadings_
p = X.shape[1]
vip = np.zeros(p)
ss_total = np.sum(Q**2 * np.sum(T**2, axis=0))
for j in range(p):
ss_j = np.sum(Q**2 * np.sum(T**2, axis=0) * W[j, :]**2)
vip[j] = np.sqrt(p * ss_j / ss_total)
return pls, scores, vip
def plot_plsda_scores(scores, y, group_names=None, figsize=(8, 6)):
"""PLS-DA スコアプロットを描画する。"""
import matplotlib.pyplot plt
fig, ax = plt.subplots(figsize=figsize)
unique = np.unique(y)
colors = plt.cm.Set1(np.linspace(, , (unique)))
color, group (colors, unique):
mask = y == group
label = group_names[group] group_names (group)
ax.scatter(scores[mask, ], scores[mask, ],
c=[color], label=label, s=, alpha=, edgecolors=)
matplotlib.patches Ellipse
cov = np.cov(scores[mask, ], scores[mask, ])
vals, vecs = np.linalg.eigh(cov)
angle = np.degrees(np.arctan2(vecs[, ], vecs[, ]))
w, h = * np.sqrt(vals * )
ell = Ellipse(xy=(scores[mask, ].mean(), scores[mask, ].mean()),
width=w, height=h, angle=angle,
fill=, color=color, linewidth=, linestyle=)
ax.add_patch(ell)
ax.set_xlabel()
ax.set_ylabel()
ax.set_title(, fontweight=)
ax.legend()
plt.tight_layout()
plt.savefig(, dpi=, bbox_inches=)
plt.close()
4. 置換検定（PLS-DA バリデーション）
def permutation_test_plsda(X, y, n_components=2, n_permutations=100):
"""
PLS-DA モデルの置換検定。
Q² と R²Y の分布を生成し、真のモデルの有意性を評価する。
"""
from sklearn.model_selection import cross_val_predict
le = LabelEncoder()
y_enc = le.fit_transform(y).astype(float)
pls_true = PLSRegression(n_components=n_components, scale=True)
y_pred = cross_val_predict(pls_true, X, y_enc, cv=5)
ss_res = np.sum((y_enc - y_pred.ravel())**2)
ss_tot = np.sum((y_enc - y_enc.mean())**2)
q2_true = 1 - ss_res / ss_tot
q2_perm = []
for _ in range(n_permutations):
y_perm = np.random.permutation(y_enc)
pls_p = PLSRegression(n_components=n_components, scale=True)
y_pred_p = cross_val_predict(pls_p, X, y_perm, cv=5)
ss_res_p = np.sum((y_perm - y_pred_p.ravel())**2)
ss_tot_p = np.sum((y_perm - y_perm.mean())**2)
q2_perm.append(1 - ss_res_p / ss_tot_p)
p_value = np.mean(np.array(q2_perm) >= q2_true)
return {
"Q2_true": q2_true,
"Q2_perm_mean": np.mean(q2_perm),
"Q2_perm_std": np.std(q2_perm),
"p_value": p_value,
"significant": p_value < 0.05,
}
5. 代謝パスウェイ濃縮解析
from scipy.stats import fisher_exact
def pathway_enrichment(significant_metabolites, pathway_annotations,
metabolite_col="Metabolite", pathway_col="Pathway",
total_metabolites=None):
"""
Fisher 正確検定による代謝パスウェイ濃縮解析。
Parameters:
significant_metabolites: list of significant metabolite names
pathway_annotations: DataFrame with metabolite-pathway mapping
total_metabolites: 解析対象の全代謝物数
"""
sig_set = set(significant_metabolites)
all_annotated = set(pathway_annotations[metabolite_col])
if total_metabolites is None:
total_metabolites = len(all_annotated)
pathways = pathway_annotations[pathway_col].unique()
results = []
for pw in pathways:
pw_members = set(
pathway_annotations[pathway_annotations[pathway_col] == pw][metabolite_col]
)
k = len(sig_set & pw_members)
K = len(pw_members)
n = len(sig_set)
N = total_metabolites
table = [[k, K - k], [n - k, N - K - n + k]]
odds_ratio, p_value = fisher_exact(table, alternative="greater")
results.append({
"Pathway": pw,
"Hits": k,
"Pathway_Size": K,
"Significant_Total": n,
"Odds_Ratio": odds_ratio,
"p_value": p_value,
})
results_df = pd.DataFrame(results).sort_values("p_value")
_, padj, _, _ = multipletests(results_df[], method=)
results_df[] = padj
results_df.to_csv(, index=)
results_df
6. 代謝物相関ネットワーク
def metabolite_correlation_network(df, metabolite_cols, method="spearman",
threshold=0.7):
"""
代謝物間の相関からネットワークを構築する。
Parameters:
threshold: |r| ≥ threshold のペアのみエッジとして採用
"""
import networkx as nx
corr = df[metabolite_cols].corr(method=method)
G = nx.Graph()
for i, met_i in enumerate(metabolite_cols):
G.add_node(met_i)
for j, met_j in enumerate(metabolite_cols):
if i < j:
r = corr.iloc[i, j]
if abs(r) >= threshold:
G.add_edge(met_i, met_j, weight=abs(r),
sign="positive" if r > 0 else "negative")
return G, corr
References
Output Files
| ファイル | 形式 |
|---|
results/univariate_results.csv | CSV |
results/vip_scores.csv | CSV |
results/pathway_enrichment.csv | CSV |
figures/plsda_scores.png | PNG |
figures/vip_barplot.png | PNG |
figures/metabolite_network.png | PNG |
利用可能ツール
ToolUniverse SMCP 経由で利用可能な外部ツール。
| カテゴリ | 主要ツール | 用途 |
|---|
| HMDB | HMDB_search | 代謝物検索 |
| HMDB | HMDB_get_metabolite | 代謝物詳細取得 |
| HMDB | HMDB_get_diseases | 代謝物-疾患関連 |
| KEGG | kegg_get_pathway_info | 代謝パスウェイ情報 |
| MetaCyc | MetaCyc_search_pathways | 代謝経路検索 |
| MetabolomicsWB | MetabolomicsWorkbench_search_compound_by_name | 代謝物データベース検索 |
参照実験
- Exp-07: PLS-DA + VIP、Pareto スケーリング、パスウェイ濃縮、相関ネットワーク
