共识聚类癌症亚型

生信喵

背景

输入TCGA数据，复现原图。

出自PMID: 29948145文章

应用场景

Consensus Clustering for cancer subtype identification，使用SNFCC +与HC和NMF算法进行分子分型，并进行比较。

CancerSubtypes包含5种计算方法对基因组数据进行癌症分子分型鉴定：

1.CC 2.CNMF 3.iCluster 4.SNF 5.WSNF

环境设置

用国内镜像安装包

options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("CancerSubtypes", version = "3.8")

加载包

library(stringr)
library(CancerSubtypes)

输入数据预处理

申明：

• 聚类本身有随机性，对数据的呈现形式没有严格要求。当数据形式变化时（或聚类方式变化时）聚类结果即出现差异，因此要得到有生物学意义的聚类结果需要对聚类数据乃至方法反复调试，包括数据标准化类型（如FPKM、 TPM、RSEM）、特征筛选过程（取高变异、生存相关等等）、距离度量方式、聚类方法等。
• 本代码不涉及具体某种类型的数据获取和统一的预处理方式。
• 此外，非负矩阵分解NMF非常耗时，因此使用的数据仅包括筛选过的少量mRNA和miRNA。

数据下载

TCGA数据的多种下载渠道：

• 所有的数据（组学和临床病理学）都可以通过TCGAbiolink下载获得，包括临床信息，参考FigureYa18（临床和突变），FigureYa22、23或34（表达数据read count、FPKM，或转换为TPM）
• 或者在xena下载你感兴趣的癌症的，除生存信息以外的所有组学数据，下载地址：https://xenabrowser.net/datapages/?hub=https://tcga.xenahubs.net:443
• 或者在xena下载处理好批次效应的TCGA和GTEx的pan-cancer数据，下载Cell处理优化过的生存资料并跟RNA-seq数据对应，参考https://bioinfoer.com/forum.php?mod=viewthread&tid=616&extra=page%3D1。作者提供的下载和预处理代码见How to get easy input from xena pan-cancer file.Rmd。
• 临床病理学和生存信息还可以在firehose上下载，下载地址：http://gdac.broadinstitute.org/
• 数据也可以在cbioportal上下载，包括组学数据和临床病理学信息，下载地址：http://www.cbioportal.org/datasets

数据格式要求

所有组学数据保持样本一致且顺序一致即可。

• 输入：列为样本，行为观测
• 临床生存信息必须和组学数据样本顺序完全一致，基于生存分析的变量筛选以及KM曲线绘制，CancerSubtypes仅利用生存时间和生存状态两个变量，因此请提前检查样本是否匹配！

读取mRNA表达数据和miRNA表达数据

# 读取mRNA
mRNA <- read.table("GeneExp.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,quote = "",row.names = 1)
mRNA[1:3, 1:3]

# 读取miRNA
miRNA <- read.table("miRNAExp.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,quote = "",row.names = 1)
miRNA[1:3, 1:3]

# 读取临床信息
surv <- read.table("survinfo.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,quote = "",row.names = 1)
surv[1:3, ]

Analysis the raw data by check the data distribution.

#所有样本及样本顺序必须完全一致，因为cancersubtype在提取生存信息的时候不关心样本名，所以要提前检查！！！

#检查组学数据样本名一致性
if(!(all(colnames(mRNA)==colnames(miRNA)) & all(colnames(mRNA)==rownames(surv)))) {
  cat("Samples mismatched! Processing now...\n")
  commonSam <- intersect(intersect(colnames(mRNA),colnames(miRNA)),rownames(surv))
  mRNA <- mRNA[,commonSam]
  miRNA <- mRNA[,commonSam]
  surv <- surv[commonSam,]
}

if((all(colnames(mRNA)==colnames(miRNA)) & all(colnames(mRNA)==rownames(surv)))) {
  cat("Samples successfully matched!\n") }

time <- as.numeric(surv$futime)
status <- as.numeric(surv$fustat) # 1: dead 0: alive

mRNA <- as.matrix(mRNA)
miRNA <- as.matrix(miRNA)

###data distribution
data.checkDistribution(mRNA)

Cancer subtype鉴定

均使用Consensus Clustering

build cancersubtype data set TCGA_target

分子分型的具体做法还是要多看文献，中心思想就是聚类

# 这里通过CancerSubtypes本身的函数通过单因素cox分析筛选与生存有关的变量
data1 <- FSbyCox(mRNA,time,status,cutoff=0.05)
# 也可以筛选top1000高变异的变量，变量的调整会对聚类结果产生影响
# data1 <- FSbyVar(mRNA, cut.type="topk",value=1000)

data2 <- FSbyCox(miRNA,time,status,cutoff=0.05)

# 也可以通过其他方式筛选变量，作为分型的输入数据
# 输入文件为列表，所有感兴趣的组学数据均可作为列表成分。
# 注意：若使用icluster聚类，最多支持5种组学数据
TCGA_target <- list(GeneExp=data1,miRNAExp=data2) 

友情科普：

共识聚类的原理是通过重抽样的思想，对原始数据集做双维度的重抽样，在扰动子集的基础上反复聚类

多次聚类结果累加求平均得到聚类共识，再在共识矩阵的基础上再聚类得最终聚类结果。

方法1. ExecuteSNF.CC

该方法是共识聚类和样本相似性网络融合的结合

result_SNFCC <- ExecuteSNF.CC(TCGA_target, clusterNum=3, K=20, alpha=0.5, t=20,
                    maxK = 5, pItem = 0.8,reps=500, 
                    title = "Molecular_Subtype_SNFCC", plot = "png", 
                    finalLinkage ="average") #基于共识矩阵的谱聚类采用平均联动方法，可换用其他如Ward.D

SNFCC_group <- result_SNFCC$group #得到的类可用作其他亚型表征分析，比如差异表达，GSEA等等，下同
SNFCC_distanceMatrix <- result_SNFCC$distanceMatrix
p_value <- survAnalysis(mainTitle="Molecular_Subtype_SNFCC",time,status,SNFCC_group,
                     SNFCC_distanceMatrix,similarity=TRUE)

# pdf("SNFCC.pdf") # 若要输出pdf，请允许下面的四行，下同。
# p_value <- survAnalysis(mainTitle="Molecular_Subtype_SNFCC",time,status,SNFCC_group,
#                      SNFCC_distanceMatrix,similarity=TRUE)
# invisible(dev.off())

方法2.CC

该方法是共识聚类和层次聚类的结合，默认使用pearson相关系数度量距离，average作为联动方法，均可换用其他

result_CC <- ExecuteCC(TCGA_target, clusterNum=3,
                    maxK = 5, pItem = 0.8,reps=500, 
                    distance = "pearson",innerLinkage = "average",
                    title = "Molecular_Subtype_CC", plot = "png", 
                    finalLinkage ="average")
CC_group <- result_CC$group
CC_distanceMatrix <- result_CC$distanceMatrix
p_value <- survAnalysis(mainTitle="Molecular_Subtype_CC",time,status,CC_group,
                     CC_distanceMatrix,similarity=TRUE)

# pdf("CC.pdf")
# p_value <- survAnalysis(mainTitle="Molecular_Subtype_CC",time,status,CC_group,
#                      CC_distanceMatrix,similarity=TRUE)
# invisible(dev.off())

生信喵

方法 3.ExecuteCNMF

该方法是层次聚类和非负矩阵分解的结合

**注意：**NMF不允许出现负值，但是在CancerSubtype包的处理下，负值会利用Kim & Tidor’s trick 法来处理

这个trick会将所有负值拉至0，所有正值扩大两倍，简单来说就是自己加上自己。

NMF的方法有许多的数学要求，对于稀疏矩阵（0较多的矩阵），往往很难进行分解而报error，所以突变数据，甚至突变签名数据请尽量避免出现在多组学聚类中，可换用离散CNA数据，或甲基化数据

此外，大n大p情况运行共识非负矩阵分解，请保持足够的耐心...建议先进行降维操作

# bug in the original code, recode and source the function
source("ExcuteCNMF.R")
result_NMF <- ExecuteCNMF(TCGA_target, clusterNum=3,nrun=5) #为了快速运行，共识仅运行5次，一般>=30次

NMF_group <- result_NMF$group
NMF_distanceMatrix <- result_NMF$distanceMatrix
p_value <- survAnalysis(mainTitle="Molecular_Subtype_NMF",time,status,NMF_group, NMF_distanceMatrix,similarity=TRUE)

# pdf("CNMF.pdf")
# p_value <- survAnalysis(mainTitle="Molecular_Subtype_NMF",time,status,NMF_group, NMF_distanceMatrix,similarity=TRUE)
# invisible(dev.off())

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250502”获取