多界微生物组机器学习助力多发性硬化症诊断¶

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, classification_report
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
import seaborn as sns

# ---------- 1. 加载数据 ----------
# 假设我们有三种界的数据框：bacteria_df, fungi_df, virus_df 每行为样本，每列为特征
# 以及标签 labels (1 = MS, 0 = Healthy)
bacteria_df = pd.read_csv('bacteria_abundance.csv', index_col=0)
fungi_df = pd.read_csv('fungi_abundance.csv', index_col=0)
virus_df = pd.read_csv('virus_abundance.csv', index_col=0)
labels = pd.read_csv('labels.csv', index_col=0)['diagnosis']  # 0/1

# ---------- 2. 特征融合 ----------
# 为每种特征加上界的前缀，以便后续追踪重要性
bacteria_df.columns = ['Bact_' + col for col in bacteria_df.columns]
fungi_df.columns = ['Fungi_' + col for col in fungi_df.columns]
virus_df.columns = ['Virus_' + col for col in virus_df.columns]

# 合并所有微生物特征，按样本索引对齐
multikingdom_df = pd.concat([bacteria_df, fungi_df, virus_df], axis=1, join='inner')

# ---------- 3. 预处理 ----------
# 处理缺失值（缺失的物种丰度填充 0）
multikingdom_df = multikingdom_df.fillna(0)
# 可选：筛选方差极低或零值过多的特征
from sklearn.feature_selection import VarianceThreshold
selector = VarianceThreshold(threshold=0.01)   # 去除方差低于 0.01 的特征
multikingdom_selected = selector.fit_transform(multikingdom_df)
multikingdom_selected = pd.DataFrame(multikingdom_selected,
                                     index=multikingdom_df.index,
                                     columns=multikingdom_df.columns[selector.get_support()])
# 标准化（针对线性模型，对树模型不必须但保留）
scaler = StandardScaler()
X = pd.DataFrame(scaler.fit_transform(multikingdom_selected),
                 index=multikingdom_selected.index, columns=multikingdom_selected.columns)
y = labels[X.index]

# ---------- 4. 划分训练集与测试集 ----------
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2,
                                                    stratify=y, random_state=42)

# ---------- 5. 训练随机森林模型 ----------
rf = RandomForestClassifier(n_estimators=500, max_depth=10, min_samples_leaf=5,
                            random_state=42, class_weight='balanced')
rf.fit(X_train, y_train)

# ---------- 6. 评估模型 ----------
y_pred_prob = rf.predict_proba(X_test)[:, 1]
auc = roc_auc_score(y_test, y_pred_prob)
print(f"多界模型测试集 AUC: {auc:.4f}")
print(classification_report(y_test, rf.predict(X_test)))

# ---------- 7. 跨验证评估（5折） ----------
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
auc_scores = cross_val_score(rf, X, y, cv=cv, scoring='roc_auc')
print(f"5折交叉验证 AUC (mean ± std): {np.mean(auc_scores):.4f} ± {np.std(auc_scores):.4f}")

# ---------- 8. 特征重要性可视化 ----------
feature_importances = pd.Series(rf.feature_importances_, index=X.columns)
top_features = feature_importances.nlargest(30)
plt.figure(figsize=(10, 8))
sns.barplot(x=top_features.values, y=top_features.index)
plt.title('Top 30 多界微生物特征重要性 (Random Forest)')
plt.xlabel('Mean Decrease in Impurity')
plt.tight_layout()
plt.savefig('feature_importance_multikingdom.png', dpi=300)
plt.show()

# ---------- 9. 比较单界模型 ----------
def evaluate_single_kingdom(abundance_df, prefix, X_test, y_test):
    """训练并评估单个界的随机森林模型"""
    X_single = pd.DataFrame(scaler.fit_transform(abundance_df),
                            index=abundance_df.index, columns=abundance_df.columns)
    # 确保与测试集样本对齐
    common_idx = X_single.index.intersection(X_test.index)
    X_train_single = X_single.loc[X_train.index, :]
    X_test_single = X_single.loc[common_idx, :]
    y_test_single = y_test.loc[common_idx]

    rf_single = RandomForestClassifier(n_estimators=500, max_depth=10,
                                       min_samples_leaf=5, class_weight='balanced')
    rf_single.fit(X_train_single, y_train.loc[X_train_single.index])
    y_prob = rf_single.predict_proba(X_test_single)[:, 1]
    auc = roc_auc_score(y_test_single, y_prob)
    print(f"{prefix} 界模型测试集 AUC: {auc:.4f}")
    return auc

print("\n单界模型性能:")
bact_auc = evaluate_single_kingdom(bacteria_df, "细菌", X_test, y_test)
fungi_auc = evaluate_single_kingdom(fungi_df, "真菌", X_test, y_test)
virus_auc = evaluate_single_kingdom(virus_df, "病毒", X_test, y_test)
print(f"\n多界融合相较最优单界的提升: {auc - max(bact_auc, fungi_auc, virus_auc):.4f}")