单个组学玩腻了，来试试20分文章喜欢的多组学联合聚类吧！

碎碎念专栏之联合多种组学数据对样本聚类分析

大家好，我是风。这里是风风的碎碎念专栏，今天我们继续碎碎念的奇妙之旅。在繁如星海的众多SCI期刊中，我最喜欢看的就是《Cancer Cell》，既不会像CNS那样让我望而生畏，又不会像OncoTarget让我扼腕叹息，组会上跟师弟师妹们吹吹水，平常在学员群里跟学员侃（zhaung）一(zhauung)侃(X)，这个期刊文章的内容呢也算比较友好，既不会难以理解，也经常可以带来灵光一现的思路，最重要是我大概还能看懂啊（毕竟再难就消化不良了）！

周二一次偶然的机会（电脑卡机）我翻到了文件夹里躺着的这篇文章，这篇文章有点久了，是2017年发在Cance Cell上面的文章。全篇文章几乎只用了一种思想，没错，就是聚类的思想。一句话总结：利用聚类的思想把同一个肿瘤的多种组学数据分别鉴定Cluster，再探讨这些cluster与特定科学问题之间的内在关系。文章我放在附件给大家，后台回复关键词就可以获取，大家有兴趣可以详细阅读一下，绝对会让你叹为观止。

今天我们的重点不是解读文章，也不是跟大家一起学校单个组学数据的聚类要怎么做，那我们要学习什么呢？我看完文章后就在想：如果把多种组学数据联合起来进行聚类，那么聚类的结果会不会更可信呢？没错，这个就是我们今天要一起学习的内容了，即：联合多种组学数据对样本进行聚类分析！

什么是聚类？

聚类分析是一种数据归约技术，旨在揭露一个数据集中观察值的子集。它可以把大量的观测值归约为若干个类，这里的类被定义为若干个观测值组成的组群，群组内观测值的相似度比群间相似度高。由于这不是一个精确的定义，从而就导致了各种聚类方法的出现。

在数据分析中，曾经有过“十大聚类分析方法”的说法，其中最基础的两种方法就是层次聚类和划分聚类。这两种方法又分别对应了多种不同的聚类算法：单联动、平均联动、K均值（K-means）和中心点划分（PAM）等等。

以上内容都是简单科普，如果想要深入学习聚类分析，建议可以看机器学习类的书籍，我比较推荐《机器学习实战》和《机器学习入门》，如果没啥兴趣的话那就当故事听一听就好了，这些都不是重点。

今天我们将以mRNA数据和miRNA数据为例，结合预后数据对上述两种数据的分析进行筛选，选择有统计学意义（p<0.05）的分子纳入进行聚类分析，所用的算法有两种，分别是高分论文中常见的Consensus Cluster（共识聚类和层次聚类的结合）和SNF-CC（样本相似性网络融合和共识聚类的结合）。

数据处理

同样我们先对读取数据进行处理，这里分别读取mRNA表达数据、miRNA表达数据和临床预后数据，样本顺序需要保持一致：

# 读取输入数据rm(list = ls())feng_mRNA <- read.table("mRNA.txt", header = T, row.names = 1, sep = "\t", check.names = F)dim(feng_mRNA)## [1] 110 208feng_mRNA[1:3, 1:3]## TCGA-EE-A29C TCGA-EE-A29D TCGA-EE-A29E## AACS 1.466886 1.09440680 1.4429054## ABCC9 1.335074 -0.30615893 -0.7257416## AC000068.5 -4.024281 -0.05549948 -0.2322607feng_miRNA <- read.table("miRNA.txt", header = T, row.names = 1, sep = "\t", check.names = F)dim(feng_miRNA)## [1] 50 208feng_miRNA[1:3, 1:3]## TCGA-EE-A29C TCGA-EE-A29D TCGA-EE-A29E## hsa-mir-140 -0.2077052 -1.820723 -0.6128286## hsa-mir-6501 -0.1076748 -1.222510 1.1499973## hsa-mir-320c-1 -1.6874758 1.008698 0.6884818# 读取临床信息clinical <- read.table("survival.txt", header = T, row.names = 1, sep = "\t", check.names = F)dim(clinical)## [1] 208 2clinical[1:3, 1:2]## futime fustat## TCGA-EE-A29C 2402 1## TCGA-EE-A29D 425 1## TCGA-EE-A29E 1940 0

共识聚类和层次聚类的结合

联合多种组学数据进行聚类分析有一个专门的R包——CancerSubtypes,听名字就知道这是一个专门用来进行分子分型的R包，这个R包内置了多种算法，并且允许输入多种组学数据，同时也可以自己设置分类数目，并且可以输入聚类图、生存曲线图、柱状图等等图形，在进行聚类分析之前，我们先对输入文件的分子进行单因素cox分析，筛选预后相关的分子，然后再对预后相关的分子进行分子分型：

# 安装R包#if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")# BiocManager::install("CancerSubtypes", version = "3.8")library(stringr)library(CancerSubtypes)# 将多组学数据转化为matrixfeng_mRNA <- as.matrix(feng_mRNA)feng_miRNA <- as.matrix(feng_miRNA)# 提取预后数据并转为数值型数据time <- as.numeric(clinical$futime)status <- as.numeric(clinical$fustat) # 通过单因素cox分析筛选与生存有关的变量feng_mRNA_1 <- FSbyCox(feng_mRNA,time,status,cutoff=0.05)head(feng_mRNA_1)## TCGA-EE-A29C TCGA-EE-A29D TCGA-EE-A29E TCGA-EE-A29G TCGA-EE-A29H## AACS 1.4668856 1.09440680 1.4429054 3.81449474 0.7244045## AC000068.5 -4.0242808 -0.05549948 -0.2322607 0.09395026 0.4682728## ADCYAP1R1 1.5042208 -0.34880485 -1.4297321 1.27292403 0.3197311## AEN 1.8920698 0.99483688 -0.1633129 -2.24824275 1.1652219## AHCY 0.9471709 -0.25826976 1.3093170 -1.15988386 -0.2344228## ALG9 1.0732856 0.09332898 0.9864030 0.60176076 -1.0161857## AC000068.5 -1.16206446 -0.83052259 -0.02954254 0.4587865 0.26982084## ADCYAP1R1 -0.48819925 -0.36269252 -0.69676842 -1.7033973 -0.70887198## AEN -0.09370338 -0.05529836 1.04563126 2.0979534 -0.74635223## AHCY 1.26663063 0.14866871 1.51950756 1.0609323 -1.09843713## ALG9 -0.64215674 -0.55737227 0.43454023 2.3730109 0.39268892## TCGA-Z2-A8RT TCGA-Z2-AA3S TCGA-Z2-AA3V## AACS 0.4515478 -0.6520251 0.7511143## AC000068.5 0.9138454 -1.5999277 -2.6565365## ADCYAP1R1 -0.3489888 0.3658171 0.5891909## AEN -2.1394648 0.6924972 1.7239421## AHCY -0.3933362 0.4666795 0.2020553## ALG9 0.1655339 2.1185737 -1.2014688dim(feng_mRNA_1) # 预后相关mRNA矩阵## [1] 58 208feng_miRNA_1 <- FSbyCox(feng_miRNA,time,status,cutoff=0.05)head(feng_miRNA_1)## TCGA-EE-A29C TCGA-EE-A29D TCGA-EE-A29E TCGA-EE-A29G TCGA-EE-A29H## hsa-mir-140 -0.2077052 -1.8207228 -0.6128286 -0.72715340 -0.23733172## hsa-mir-29c -0.3060311 0.6886091 -1.1958765 -0.43787853 -0.76284173## hsa-mir-3615 -1.1811729 2.5013622 -0.3787177 -0.18641524 0.11262520## hsa-mir-342 0.3266871 0.2512507 -2.3404871 -0.08497941 -0.57767513## hsa-mir-150 0.3948697 -0.7893406 -1.6516306 -0.00912708 0.06285036## hsa-mir-6715a -1.4300915 -0.2108727 -1.4300915 0.01182342 0.14626088## TCGA-EE-A29L TCGA-EE-A29M TCGA-EE-A29N TCGA-EE-A29P TCGA-EE-A29Q## hsa-mir-140 0.09755653 -0.4012395 -1.1014861 -1.4545895 0.35182982## hsa-mir-29c -1.50669180 -0.3747896 1.2331096 -0.6983878 -0.62047055## hsa-mir-3615 0.87996183 0.7572464 0.1790880 0.5933491 0.50246263## hsa-mir-342 -2.1955683 -0.8463859 -0.3670937 -1.3057328 -0.4374958## hsa-mir-150 -1.6870735 -0.8004644 0.1981141 -0.6303854 0.4926376## hsa-mir-6715a 0.4087208 -0.7819786 1.0311021 -0.3264904 -0.4989945## TCGA-GN-A4U7 TCGA-GN-A4U8 TCGA-GN-A4U9 TCGA-GN-A8LK TCGA-GN-A8LL## hsa-mir-140 0.1264486 0.3575642 0.04938648 0.5110366 0.3982933## hsa-mir-29c -0.5115752 0.6207672 -0.74725730 -0.7413471 -1.0342885## hsa-mir-3615 -0.1107457 -0.8100202 -0.05717685 -1.6027796 -0.7239730## hsa-mir-342 -0.6084463 0.8537676 1.04082144 -1.4733982 0.5491078## hsa-mir-150 -0.8621921 0.9759203 -0.54453511 -1.9705866 -1.4052064## hsa-mir-6715a 1.1590540 -0.2220695 -1.43009153 1.4062901 2.8071632## TCGA-GN-A9SD TCGA-LH-A9QB TCGA-OD-A75X TCGA-QB-A6FS TCGA-QB-AA9O## hsa-mir-140 -0.1627846 2.0326293 -0.3213311 -0.4639647 -1.085978296## hsa-mir-29c 0.6297534 -0.1876349 -1.1310561 -1.2804463 -0.618297701## hsa-mir-3615 0.7508907 3.4250712 0.8917478 -1.6341735 -0.810048696## hsa-mir-342 1.2044722 -0.5564898 -0.8469227 1.0796945 0.006063361## hsa-mir-150 1.6244359 -0.3588959 -0.8994351 0.9273331 0.005589674## hsa-mir-6715a 0.8715719 -1.4300915 -0.3152921 0.4024280 -0.465558238## TCGA-RP-A690 TCGA-RP-A693 TCGA-RP-A694 TCGA-W3-A824 TCGA-W3-A825## hsa-mir-140 3.04077333 -0.024647311 -0.42729111 -0.8584655 -0.001319206## hsa-mir-29c 0.03350409 0.005983904 -0.09795682 -0.7038038 -0.424495873## hsa-mir-3615 -0.40482638 -0.698881174 -0.91143669 0.6803995 -0.225036428## hsa-mir-342 0.06728260 1.426602997 1.70095371 0.6615799 0.139623050## hsa-mir-150 -0.49461830 0.106000471 -0.08050157 0.3363975 0.065289692## hsa-mir-6715a 1.08305870 -0.597943618 -0.80955121 -0.5720382 0.455363229## TCGA-W3-A828 TCGA-W3-AA1O TCGA-W3-AA1Q TCGA-W3-AA1R TCGA-W3-AA1V## hsa-mir-140 1.15271429 -0.775830477 0.63841948 -0.92029580 1.0108726## hsa-mir-29c -0.04006374 -0.436159935 0.09398823 0.02933264 -1.0658215## hsa-mir-3615 -2.65102753 1.027526381 -0.64889569 -0.77798941 0.5179123## hsa-mir-342 0.42387614 -0.285329320 0.18156686 0.01951699 -0.7402158## hsa-mir-150 -0.18125457 0.009359856 0.64680582 0.18736026 -0.6888952## hsa-mir-6715a -1.43009153 -0.363195196 -0.18377747 0.01531002 -0.8672379## TCGA-W3-AA1W TCGA-W3-AA21 TCGA-WE-A8K1 TCGA-WE-A8K5 TCGA-WE-A8K6## hsa-mir-140 0.004245305 0.8758883 -0.2480779 -0.69594648 0.5058502## hsa-mir-29c 1.228820947 -0.3953445 0.2987257 -1.55231201 0.9218683## hsa-mir-3615 0.754882764 0.5124453 0.2917846 -0.02029873 -0.1248402## hsa-mir-342 1.633174183 -0.2803311 0.7375962 -0.43718896 0.1411473## hsa-mir-150 1.312661577 0.1342371 0.8905994 0.18944587 -1.3466881## hsa-mir-6715a -0.855274127 0.3804002 -0.7140988 0.29732621 -0.1263051## TCGA-WE-A8ZM TCGA-WE-A8ZT TCGA-WE-A8ZY TCGA-WE-AA9Y TCGA-WE-AAA0## hsa-mir-140 4.1910133 1.8605198 0.3152588 1.1724176 0.2957501## hsa-mir-29c 0.4180620 -0.6780236 -0.9901888 0.2178403 0.2482231## hsa-mir-3615 -0.3932321 0.4734871 -0.9240524 -0.1905569 0.6322070## hsa-mir-342 -0.4537503 -0.7781260 -1.0766778 0.9910582 0.6594354## hsa-mir-150 -1.6488570 0.5330026 -0.3488405 0.6182259 1.2462586## hsa-mir-6715a 0.8721599 -0.4178828 0.7091166 -1.4300915 -0.1457442## TCGA-WE-AAA3 TCGA-WE-AAA4 TCGA-XV-AB01 TCGA-YD-A89C TCGA-YD-A9TA## hsa-mir-140 -0.64952831 -0.44429921 -0.24989871 0.6197881 -0.8516525## hsa-mir-29c -0.47689725 -0.26154911 0.02368495 -2.6673506 -0.3786012## hsa-mir-3615 0.13322659 0.78024160 -2.68799359 0.5214133 -0.8439531## hsa-mir-342 0.12646119 0.80249415 -1.47154679 -0.9195341 0.2323464## hsa-mir-150 0.22422716 0.60700834 -0.27933165 -1.4884894 -0.1460285## hsa-mir-6715a 0.09153248 -0.07661256 -0.55110128 1.3838606 -1.4300915## TCGA-Z2-A8RT TCGA-Z2-AA3S TCGA-Z2-AA3V## hsa-mir-140 -0.9894110 0.07225548 -0.7313071## hsa-mir-29c 0.9240470 -0.82913918 0.7250476## hsa-mir-3615 0.1302950 -0.00081556 -1.1466755## hsa-mir-342 -0.2402631 -0.92905475 0.5239246## hsa-mir-150 0.6511516 -0.11056251 0.8675116## hsa-mir-6715a 0.7442022 1.23585850 0.3801821dim(feng_miRNA_1) # 预后相关miRNA矩阵## [1] 29 208

接下来我们先看看Consensus Cluster算法的应用，这里使用ExecuteCC：

# cancersubtype支持将多组学数据合并为一个list再进行聚类TCGA_target<- list(Omics1 = feng_mRNA_1,Omics2 = feng_miRNA_1) # 进行Consensus Cluster聚类CC<- ExecuteCC(TCGA_target,clusterNum=2,maxK = 10, clusterAlg = "hc",distance = "pearson",title = "Molecular_Subtype_CC", reps = 500,pItem = 0.8, pFeature = 1, plot = "png", innerLinkage = "average",finalLinkage = "average", verbose = FALSE, corUse = "everything")# 提取聚类结果的样本分类CC_group<- CC$group# 提取聚类结果矩阵CC_distanceMatrix<- CC$distanceMatrix# 提取原始结果CC_originalResult<- CC$originalResult# 进行生存分析并绘图Survival_Analysis<- survAnalysis(mainTitle="Survival Analysis",time,status,CC_group,CC_distanceMatrix,similarity=TRUE)## ## *****************************************************## Survival Analysis Cluster= 2 Call:## survdiff(formula = Surv(time, status) ~ group)## ## N Observed Expected (O-E)^2/E (O-E)^2/V## group=1 114 76 56.3 6.87 13.2## group=2 94 48 67.7 5.72 13.2## ## Chisq= 13.2 on 1 degrees of freedom, p= 3e-04

Survival_Analysis# 生存分析P值，发现两组分型之间存在显著预后差异## [1] 0.0002872837

简单明了直接，一次性解决分子分型并绘制生存曲线绘制分型的差异性，不用再自己单独写代码合并数据和分型结果，并且这种分型是综合了多个组学数据的结果，比单个组学结果更加准确明了。接下来我们看第二种聚类方法：样本相似性网络融合和共识聚类的结合

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

样本相似性网络融合SNF是以样本相似性算法为基础开发的一种聚类算法，而CancerSubtypes提供了一种新的算法，是将SNF与前面的ConsensusCluster算法联合运用的算法，具体的算法原理我们不做展开，这里主要使用ExecuteSNF.CC，我们来看一下：

SNFCC<- ExecuteSNF.CC(TCGA_target, clusterNum=2, K = 20, alpha = 0.5, t = 20,maxK = 10, pItem = 0.8, reps = 500, title = "ConsensusClusterResult",plot = "png", finalLinkage = "average") SNFCC_group<- SNFCC$group #得到的类可用作其他亚型表征分析，比如差异表达，GSEA等等，下同SNFCC_distanceMatrix<- SNFCC$distanceMatrixSNFCC_originalResult<- SNFCC$originalResultSurvival_Analysis_p_value<- survAnalysis(mainTitle="Survival Analysis",time,status,SNFCC_group,SNFCC_distanceMatrix,similarity=TRUE)## ## *****************************************************## Survival Analysis Cluster= 2 Call:## survdiff(formula = Surv(time, status) ~ group)## ## N Observed Expected (O-E)^2/E (O-E)^2/V## group=1 113 80 55.4 10.97 20.6## group=2 95 44 68.6 8.85 20.6## ## Chisq= 20.6 on 1 degrees of freedom, p= 6e-06

Survival_Analysis_p_value# 生存分析P值，发现两组分型之间存在显著预后差异## [1] 5.748777e-06

同样我们就得到了SNF+CC两种算法联合聚类的分型结果，这个结果看起来比单纯使用ConsensusCluster更好一点，后续可以根据分型结果再展开分析。

需要注意的是，尽管我们学习了多种组学数据的共同聚类结果，但是我们仍然要对这种做法慎之又慎，首先要考虑的是纳入联合分析的组学数据，如何对数据进行矫正，不同的数据类型有不同的矫正方法；其次，多种组学数据纳入分型的意义是什么，能不能解释清楚相应结果？再者，哪种算法最适合纳入的组学数据，既不会过度拟合，又能够防止拟合结果不好？最后最后最后，聚类结果解读和验证，这是最重要的一点，也是被大部分文章所忽略的一点，如何保证你得到的“类”是“真实的类”并且其稳健又可重复？这就需要结合具体情况再进行分析验证了，不过这是另外一个故事啦。

最近在写毕业论文和对一篇文章进行返修，提出意见的几位审稿人专业知识让我叹为观止，整整提了40多页word文档的问题，从统计到分析细节，简直让我怀疑人生，不过也让我受益匪浅。果然不能总是差不多，不然关键时刻就总是差一点，是时候学习下强迫症，一个标点符号都不要放过了！

后台回复“feng35”获得本次代码和数据，我们下次再见吧，晚安！*^_^*