664 微生物组临床转化¶

# 微生物组诊断标志物发现流程
import pandas as pd  # 数据处理
import numpy as np  # 数值计算
from sklearn.ensemble import RandomForestClassifier  # 随机森林
from sklearn.model_selection import (
    RepeatedStratifiedKFold, cross_val_score  # 交叉验证
)
from sklearn.metrics import roc_auc_score, classification_report  # 评估
from sklearn.preprocessing import StandardScaler  # 标准化
import shap  # SHAP解释性分析

# 1. 数据准备
otu = pd.read_csv("otu_table.csv", index_col=0).T  # OTU表（样本×物种）
meta = pd.read_csv("metadata.csv")  # 元数据
y = (meta["diagnosis"] == "CRC").astype(int)  # 标签：结直肠癌vs健康

# 2. CLR转换
from scipy.stats import gmean
otu_pseudo = otu + 0.5  # 加伪计数
otu_clr = np.log(otu_pseudo.div(gmean(otu_pseudo, axis=1), axis=0))  # CLR

# 3. 多中心交叉验证（模拟真实临床验证）
cv = RepeatedStratifiedKFold(n_splits=5, n_repeats=10, random_state=42)  # 重复CV

# 4. 训练分类器
clf = RandomForestClassifier(
    n_estimators=500,  # 500棵树
    max_depth=10,  # 最大深度
    min_samples_leaf=5,  # 叶节点最少样本
    class_weight="balanced",  # 处理类别不平衡
    random_state=42
)

# 5. 评估
scores = cross_val_score(clf, otu_clr, y, cv=cv, scoring="roc_auc")
print(f"AUC: {np.mean(scores):.3f} ± {np.std(scores):.3f}")

# 6. SHAP分析找关键标志物
clf.fit(otu_clr, y)  # 全数据训练
explainer = shap.TreeExplainer(clf)  # SHAP解释器
shap_values = explainer.shap_values(otu_clr)  # 计算SHAP值

# Top20标志物
feature_importance = np.abs(shap_values[1]).mean(axis=0)  # 平均绝对SHAP值
top20 = pd.Series(feature_importance, index=otu_clr.columns).nlargest(20)
print("Top20诊断标志物：")
print(top20)

# 7. SHAP可视化
shap.summary_plot(shap_values[1], otu_clr, max_display=20)

# 外部队列验证——临床转化的金标准
# 训练集：队列A，验证集：队列B（独立医院/国家）

# 训练
train_otu = pd.read_csv("cohort_A_otu.csv", index_col=0).T
train_y = pd.read_csv("cohort_A_labels.csv")["disease"]

clf.fit(train_otu_clr, train_y)  # 在队列A训练

# 外部验证
test_otu = pd.read_csv("cohort_B_otu.csv", index_col=0).T
test_y = pd.read_csv("cohort_B_labels.csv")["disease"]

# 确保特征一致
common_features = train_otu.columns.intersection(test_otu.columns)
test_pred = clf.predict_proba(test_otu_clr[common_features])[:, 1]

# 评估
ext_auc = roc_auc_score(test_y, test_pred)
print(f"外部验证AUC: {ext_auc:.3f}")
print(classification_report(test_y, (test_pred > 0.5).astype(int)))

# 临床实用性：要求外部验证AUC > 0.80才有临床价值
if ext_auc > 0.80:
    print("达到临床应用标准")
else:
    print("需要更多数据或改进模型")

# 标准化粪便采集流程（临床研究SOP）
# 1. 采集要求
#    - 使用DNA/RNA Shield稳定液
#    - 室温保存≤72小时 或 -80°C长期保存
#    - 避免freeze-thaw循环

# 2. DNA提取标准
#    - 推荐：Qiagen QIAamp PowerFecal Pro DNA Kit
#    - 必须包含机械裂解步骤（bead-beating）
#    - 提取空白对照（每批次1个）

# 3. 测序标准
#    - 16S: V3-V4区，双端300bp，>10,000 reads/样本
#    - 宏基因组: >5 Gbp/样本
#    - 每批次包含模拟群落（mock community）

# 4. 分析流程
#    - 使用版本锁定的分析管线
#    - DADA2/Deblur（16S ASV推断）
#    - MetaPhlAn4/Kraken2（宏基因组分类）

# 临床微生物组数据质量控制
def clinical_qc_check(otu_table, metadata, mock_expected):
    """临床微生物组数据质量控制"""
    qc_results = {}  # 存储QC结果

    # 1. 测序深度检查
    read_counts = otu_table.sum(axis=0)  # 每个样本的总reads
    qc_results["min_reads"] = read_counts.min()
    qc_results["samples_below_10k"] = (read_counts < 10000).sum()
    print(f"最低reads: {read_counts.min()}")
    print(f"低于10K reads的样本: {qc_results['samples_below_10k']}")

    # 2. 模拟群落验证
    mock_samples = metadata[metadata["sample_type"] == "mock"].index
    if len(mock_samples) > 0:
        mock_data = otu_table[mock_samples]
        # 与预期组成比较
        from scipy.stats import pearsonr
        for sample in mock_samples:
            observed = mock_data.loc[sample]
            r, p = pearsonr(observed, mock_expected)
            qc_results[f"mock_{sample}_r"] = r
            print(f"模拟群落 {sample}: Pearson r = {r:.3f}")

    # 3. 提取空白检查
    blank_samples = metadata[metadata["sample_type"] == "blank"].index
    if len(blank_samples) > 0:
        blank_reads = otu_table[blank_samples].sum(axis=0)
        qc_results["blank_max_reads"] = blank_reads.max()
        print(f"空白对照最大reads: {blank_reads.max()}")
        if blank_reads.max() > 1000:
            print("警告：空白对照污染较高！")

    # 4. 批次效应检查
    from scipy.spatial.distance import braycurtis
    # 计算批次间vs批次内距离比
    # ... (参见662_微生物组可重复性分析.md)

    return qc_results

临床流程：
1. 患者提供粪便样本
2. 宏基因组测序
3. ML模型分析菌群特征
4. 输出CRC风险评分（0-100）
5. 高风险者建议肠镜检查

关键标志物（多项研究共识）：
- Fusobacterium nucleatum ↑（核梭杆菌增加）
- Parvimonas micra ↑
- Peptostreptococcus stomatis ↑
- Bifidobacterium ↓（双歧杆菌减少）
- Roseburia ↓

AUC通常在0.80-0.85，联合FIT可达0.90+

# FMT供体微生物组质量评估
def evaluate_fmt_donor(donor_profile):
    """评估FMT供体的微生物组质量"""
    score = 0  # 初始评分

    # 1. 多样性检查
    shannon = -sum(p * np.log(p) for p in donor_profile if p > 0)
    if shannon > 3.5:  # Shannon多样性>3.5
        score += 25
        print(f"多样性: {shannon:.2f} (优秀)")
    elif shannon > 3.0:
        score += 15
        print(f"多样性: {shannon:.2f} (良好)")

    # 2. 有益菌检查
    beneficial = ["Faecalibacterium", "Roseburia", "Bifidobacterium",
                   "Akkermansia", "Eubacterium"]
    for genus in beneficial:
        if donor_profile.get(genus, 0) > 0.01:  # 相对丰度>1%
            score += 5
            print(f"有益菌 {genus}: 存在")

    # 3. 病原菌排除
    pathogens = ["Clostridioides_difficile", "Salmonella",
                  "Campylobacter", "ESBL_producing"]
    for pathogen in pathogens:
        if donor_profile.get(pathogen, 0) > 0:
            score -= 50  # 有病原菌直接扣分
            print(f"警告：检测到病原菌 {pathogen}！")

    # 4. 耐药基因检查
    # 使用AMRFinderPlus/RGI检测耐药基因
    # 高风险耐药基因存在 → 供体不合格

    print(f"供体评分: {score}/100")
    return score

# 临床微生物组分析流程
样本采集(标准化SOP) → DNA提取(含对照)
  → 测序(含mock) → 质量控制
  → 生物信息分析 → 标志物发现
  → 内部交叉验证 → 外部队列验证
  → 临床试验 → 监管审批

# 关键标准和指南
STORMS: 微生物组研究报告检查清单
2024 Lancet共识: 临床微生物组检测国际共识
IHMS: 人体微生物组研究标准

# 当前FDA获批的微生物组产品
REBYOTA (2022): 首个获批FMT产品(CDI)
VOWST (2023): 首个口服FMT胶囊(CDI)

# 临床转化的AUC标准
AUC < 0.70: 无临床价值
AUC 0.70-0.80: 有一定价值，需与现有方法联合
AUC 0.80-0.90: 有临床应用价值
AUC > 0.90: 优秀的诊断性能