代谢组学多元统计（PLS-DA / OPLS-DA）¶

# === PLS-DA和OPLS-DA分析（使用scikit-learn和pyopls） ===
import numpy as np
import pandas as pd
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_predict, LeaveOneOut
from sklearn.preprocessing import LabelEncoder
import matplotlib.pyplot as plt

# 1. 准备数据
data = pd.read_csv("metabolomics_normalized.csv", index_col=0)  # 读取归一化后的数据
X = data.iloc[:, 1:].values  # 代谢物矩阵（去掉分组列）
y_labels = data.iloc[:, 0].values  # 分组标签（如"Disease", "Control"）
le = LabelEncoder()
y = le.fit_transform(y_labels)  # 将标签转为0/1编码

# 2. PLS-DA建模
pls = PLSRegression(n_components=2, scale=True)  # 2个主成分，自动缩放
pls.fit(X, y)  # 拟合模型

# 3. 获取得分（Scores）用于绘制得分图
scores = pls.x_scores_  # 每个样本在主成分上的投影坐标

# 4. 绘制PLS-DA得分图
fig, ax = plt.subplots(figsize=(8, 6))
for label in np.unique(y):
    mask = y == label  # 筛选当前组的样本
    ax.scatter(scores[mask, 0], scores[mask, 1],
               label=le.inverse_transform([label])[0], s=60, alpha=0.7)
ax.set_xlabel("Component 1")
ax.set_ylabel("Component 2")
ax.legend()
ax.set_title("PLS-DA Score Plot")
plt.savefig("plsda_scores.png", dpi=150)

# 5. 交叉验证评估Q2
loo = LeaveOneOut()  # 留一法交叉验证
y_pred_cv = cross_val_predict(pls, X, y, cv=loo)  # 交叉验证预测
ss_res = np.sum((y - y_pred_cv) ** 2)  # 残差平方和
ss_tot = np.sum((y - y.mean()) ** 2)  # 总平方和
q2 = 1 - ss_res / ss_tot  # Q2计算公式
print(f"Q2 = {q2:.3f}")  # Q2>0.5说明模型有预测能力

# 6. 计算VIP值
def calculate_vip(pls_model):
    """计算PLS-DA的VIP值"""
    t = pls_model.x_scores_        # 得分矩阵
    w = pls_model.x_weights_       # 权重矩阵
    q = pls_model.y_loadings_      # Y载荷
    p, h = w.shape                 # p=变量数, h=主成分数
    vips = np.zeros(p)
    s = np.diag(t.T @ t @ q.T @ q).reshape(h, -1)  # 各成分贡献
    total_s = np.sum(s)
    for i in range(p):
        weight = np.array([(w[i, j] / np.linalg.norm(w[:, j])) ** 2 for j in range(h)])
        vips[i] = np.sqrt(p * np.sum(s.T * weight) / total_s)
    return vips

vip_scores = calculate_vip(pls)
# 筛选VIP>1的代谢物
important_features = data.columns[1:][vip_scores > 1]
print(f"VIP>1的代谢物数量: {len(important_features)}")