RNA-seq分析之共识聚类分析

这个主要是用来对RNA-seq的数据进行划分（例如肿瘤细分），根据划分结果分组，进行分析，常见于冷热肿瘤

# 肿瘤样本共识聚类分析（ConsensusClusterPlus）
setwd("G:/cold_hot tumor/parting")  # 设置工作目录# 安装并加载所需包
if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
if (!require("ConsensusClusterPlus", quietly = TRUE)) {BiocManager::install("ConsensusClusterPlus")
}
library(ConsensusClusterPlus)# 读取表达数据（行=基因，列=样本）
expr_data <- read.delim("G:/cold_hot tumor/Preliminarily processed data/Count.txt",header = TRUE,check.names = FALSE,  # 保留原始样本名格式row.names = 1         # 第一列作为行名（基因名）
)# 数据预处理
## 检查数据基本信息
cat("原始数据维度（基因数×样本数）：", dim(expr_data), "\n")
if (any(is.na(expr_data))) {stop("数据中存在缺失值，请先处理缺失值！")
}## 筛选高变基因（取前5000个MAD最大的基因，提高聚类效率）
mads <- apply(expr_data, 1, mad)  # 计算每行（基因）的中位数绝对偏差
expr_filtered <- expr_data[rev(order(mads))[1:5000], ]  # 按MAD降序取前5000个基因## 行中心化（每个基因减去其中位数，减少基线差异影响）
expr_scaled <- sweep(expr_filtered, 1, apply(expr_filtered, 1, median, na.rm = TRUE),  # 按行减中位数FUN = "-"
)# 共识聚类分析
## 转换为矩阵（ConsensusClusterPlus要求输入矩阵格式）
expr_matrix <- as.matrix(expr_scaled)## 运行聚类（关键参数：maxK=20尝试最多20个聚类，reps=1000次重复）
cc_result <- ConsensusClusterPlus(d = expr_matrix,          # 输入矩阵maxK = 20,                # 最大聚类数reps = 1000,              # 重复次数（越大结果越稳定）pItem = 0.8,              # 每次抽样保留80%样本pFeature = 1,             # 保留所有特征基因clusterAlg = "hc",        # 聚类算法：层次聚类distance = "pearson",     # 距离度量：皮尔逊相关系数title = "consensus_clusters",  # 结果保存文件夹名plot = "png",             # 输出聚类相关图表（png格式）writeTable = FALSE        # 不重复写入中间文件
)# 计算并可视化聚类一致性（帮助确定最佳聚类数k）
icl_result <- calcICL(cc_result, title = "icl_results", plot = "png")# 提取聚类结果（根据实际情况选择最佳k值，这里以k=3为例）
best_k <- 3  # 请根据一致性矩阵和ICL图调整最佳k值
cluster_info <- cc_result[[best_k]]$consensusClass  # 样本-聚类对应关系
cat("\n最佳聚类数k=", best_k, "时的样本分布：\n")
print(table(cluster_info))  # 展示每个聚类的样本数量# 整理并保存结果
## 构建样本-分组数据框
group_df <- data.frame(sample = names(cluster_info),cluster = paste0("Cluster", cluster_info),  # 聚类命名为Cluster1,2,...row.names = names(cluster_info),check.names = FALSE
)## 合并聚类信息与原始表达数据（可选，按需使用）
expr_with_cluster <- rbind(Cluster = group_df$cluster,  # 首行添加聚类信息expr_data  # 原始表达数据
)## 保存结果文件
# 1. 样本聚类信息
write.table(group_df,file = "sample_clusters.txt",sep = "\t",row.names = TRUE,col.names = TRUE,quote = FALSE
)# 2. 带聚类信息的表达矩阵（可选）
write.table(expr_with_cluster,file = "expr_with_clusters.txt",sep = "\t",row.names = TRUE,col.names = NA,  # 避免列名偏移quote = FALSE
)cat("\n分析完成！结果文件已保存至工作目录：\n")
cat("1. 样本聚类信息：sample_clusters.txt\n")
cat("2. 带聚类信息的表达矩阵：expr_with_clusters.txt\n")