407_多组学降维方法比较¶

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

np.random.seed(42)

# ========== 模拟多组学数据 ==========
n_samples = 200   # 样本数
n_rna = 1000      # RNA特征数
n_protein = 300   # 蛋白特征数

# 生成3个亚群（模拟不同细胞类型）
labels = np.repeat([0, 1, 2], [70, 70, 60])  # 3组样本

def make_omics_data(n_samples, n_features, labels, effect_size=3):
    """生成有组别差异的模拟多组学数据"""
    data = np.random.randn(n_samples, n_features)
    # 对前50个特征加入组间差异信号
    for group in range(3):
        mask = labels == group
        data[mask, :50] += effect_size * group  # 不同组在前50个特征有差异
    return data

rna_data = make_omics_data(n_samples, n_rna, labels)
protein_data = make_omics_data(n_samples, n_protein, labels, effect_size=2)

# ========== 1. PCA降维 ==========
print("===== PCA =====")
scaler = StandardScaler()
rna_scaled = scaler.fit_transform(rna_data)  # Z-score标准化

pca = PCA(n_components=50)  # 先降到50维（中间步骤）
pca_result = pca.fit_transform(rna_scaled)

# 方差解释量（Scree plot）
var_explained = pca.explained_variance_ratio_
print(f"PC1解释方差: {var_explained[0]:.3f}")
print(f"PC2解释方差: {var_explained[1]:.3f}")
print(f"前10个PC累计解释方差: {var_explained[:10].sum():.3f}")

# ========== 2. tSNE降维 ==========
print("\n===== tSNE =====")
from sklearn.manifold import TSNE

# 先用PCA降到50维，再做tSNE（节省时间，推荐做法）
tsne = TSNE(
    n_components=2,    # 降到2维（用于可视化）
    perplexity=30,     # 邻域大小，通常5-50
    learning_rate=200, # 学习率
    max_iter=1000,    # 迭代次数
    random_state=42
)
tsne_result = tsne.fit_transform(pca_result[:, :50])  # 输入前50个PC
print(f"tSNE输出形状: {tsne_result.shape}")

# ========== 3. UMAP降维（推荐！） ==========
print("\n===== UMAP =====")
# pip install umap-learn
import umap

umap_model = umap.UMAP(
    n_components=2,    # 输出2维
    n_neighbors=15,    # 邻居数：小=更局部，大=更全局
    min_dist=0.1,      # 点间最小距离：小=更紧凑的聚类
    metric="euclidean", # 距离度量
    random_state=42
)
umap_result = umap_model.fit_transform(pca_result[:, :50])
print(f"UMAP输出形状: {umap_result.shape}")

# ========== 4. 多组学联合降维（数据拼接策略） ==========
print("\n===== 多组学联合UMAP =====")
# 简单策略：将多种组学数据拼接后一起降维
rna_scaled_pca = PCA(n_components=30).fit_transform(
    StandardScaler().fit_transform(rna_data)
)
protein_scaled_pca = PCA(n_components=15).fit_transform(
    StandardScaler().fit_transform(protein_data)
)

# 水平拼接（保持样本数不变，特征数相加）
combined = np.hstack([rna_scaled_pca, protein_scaled_pca])
print(f"拼接后形状: {combined.shape}")  # (200, 45)

combined_umap = umap.UMAP(n_components=2, random_state=42).fit_transform(combined)

# ========== 5. 可视化对比四种方法 ==========
colors = {0: "#E74C3C", 1: "#2ECC71", 2: "#3498DB"}
color_list = [colors[l] for l in labels]

fig, axes = plt.subplots(2, 2, figsize=(14, 12))
axes = axes.flatten()

datasets = [
    (pca_result[:, :2], "PCA (PC1 vs PC2)"),
    (tsne_result, "tSNE (perplexity=30)"),
    (umap_result, "UMAP (RNA only)"),
    (combined_umap, "UMAP (RNA+Protein)")
]

for ax, (data, title) in zip(axes, datasets):
    for group in range(3):
        mask = labels == group
        ax.scatter(
            data[mask, 0], data[mask, 1],
            c=colors[group], label=f"Group {group}",
            alpha=0.7, s=30, edgecolors="none"
        )
    ax.set_title(title, fontsize=12, fontweight="bold")
    ax.legend(loc="upper right", markerscale=1.5)
    ax.set_xlabel("Dim 1")
    ax.set_ylabel("Dim 2")

plt.suptitle("Dimensionality Reduction Methods Comparison", fontsize=14)
plt.tight_layout()
plt.savefig("dim_reduction_comparison.png", dpi=150)
print("\n对比图已保存: dim_reduction_comparison.png")