肿瘤分子亚型分类¶

# === 基于 NMF 的肿瘤分子亚型发现 ===
library(NMF)                           # 非负矩阵分解包

# 1. 准备基因表达矩阵（非负值，如 TPM）
expr <- as.matrix(read.table(
  "expression_tpm.txt",               # 输入TPM矩阵
  header = TRUE, row.names = 1))

# 筛选高变异基因（减少噪音）
gene_var <- apply(expr, 1, var)        # 计算每个基因的方差
top_genes <- names(sort(gene_var,      # 按方差排序
                        decreasing = TRUE))[1:2000]  # 取前2000
expr_sub <- expr[top_genes, ]          # 子集矩阵

# 2. 运行 NMF（测试不同 k 值）
estim <- nmf(expr_sub,                 # 输入矩阵
             rank = 2:6,               # 测试2-6个亚型
             method = "brunet",        # Brunet算法
             nrun = 30,                # 每个k运行30次
             seed = 42)                # 随机种子

# 3. 选择最优 k 值
plot(estim)                            # 绘制cophenetic等指标
# 选择 cophenetic 系数最高且开始下降前的 k

# 4. 用最优 k 运行最终 NMF
best_k <- 4                           # 假设最优k=4
final_nmf <- nmf(expr_sub,            # 输入矩阵
                 rank = best_k,        # 最优k
                 method = "brunet",    # 算法
                 nrun = 100,           # 100次运行取共识
                 seed = 42)

# 5. 提取亚型分配
subtypes <- predict(final_nmf)         # 每个样本的亚型
table(subtypes)                        # 各亚型样本数

# 6. 可视化
consensusmap(final_nmf)                # 共识矩阵热图

# === 乳腺癌 PAM50 分类 ===
library(genefu)                        # 基因签名包

# PAM50 分类
pam50_result <- molecular.subtyping(
  sbt.model = "pam50",                # PAM50模型
  data = t(expr),                      # 转置（样本为行）
  annot = gene_annotation,            # 基因注释（需要EntrezID）
  do.mapping = TRUE                    # 自动映射基因ID
)

# 查看结果
table(pam50_result$subtype)            # 各亚型人数

# === 共识聚类确定最优分型数 ===
from sklearn.cluster import KMeans     # K-Means聚类
import numpy as np                     # 数值计算
import pandas as pd                    # 数据处理

def consensus_clustering(data, k_range=range(2, 7), n_iter=100):
    """简化版共识聚类"""
    results = {}                       # 存储结果
    n_samples = data.shape[0]          # 样本数

    for k in k_range:                  # 遍历每个k值
        consensus = np.zeros((n_samples, n_samples))  # 共识矩阵
        count = np.zeros((n_samples, n_samples))      # 共现次数

        for _ in range(n_iter):        # 重复n次
            # 随机抽样80%的样本
            idx = np.random.choice(n_samples,
                                   size=int(0.8 * n_samples),
                                   replace=False)
            sub_data = data[idx]       # 子集数据

            # K-Means 聚类
            km = KMeans(n_clusters=k, n_init=10, random_state=None)
            labels = km.fit_predict(sub_data)  # 聚类标签

            # 更新共识矩阵
            for i in range(len(idx)):  # 遍历被选中的样本
                for j in range(i+1, len(idx)):
                    count[idx[i], idx[j]] += 1  # 共同被选中次数
                    if labels[i] == labels[j]:   # 如果分在同一类
                        consensus[idx[i], idx[j]] += 1  # 共识+1

        # 归一化
        consensus = np.divide(consensus, count,
                              where=count > 0)  # 共识/共现
        consensus += consensus.T       # 对称化
        np.fill_diagonal(consensus, 1) # 对角线=1

        # 计算 CDF 下面积（评估聚类稳定性）
        flat = consensus[np.triu_indices(n_samples, k=1)]
        results[k] = {"matrix": consensus, "auc": np.mean(flat)}
        print(f"k={k}, CDF AUC={np.mean(flat):.4f}")

    return results

# 使用
# result = consensus_clustering(expression_matrix.T)