终于找到发高分生信文章的财富密码了！（附代码！）

分子分型

大家好，我是风风。上一次我们看了机器学习中三种基于线性模型的特征筛选方法，今天我们来聊一聊在高分生信文章中经常见到的生信分析手段——共识聚类。

高分文章中怎么用共识聚类的呢？我们来看两个例子：

第一个例子是2018年发表在Immunity上的"The Immune Landscape of Cancer"，这是非常经典的一篇文章，作者利用TCGA数据库的泛癌数据，将千奇百怪的肿瘤和繁杂的高通量数据联合前来，通过共识聚类将肿瘤分为6种免疫亚型，并且寻找了每一种亚型的分子特征（如下图）：

第二个例子是2019年发表在Cell上的"Integrated Proteogenomic Characterization of HBV Related Hepatocellular Carcinoma"，利用蛋白质组学数据将肝癌分成了三个亚组，通过多组学分析最后找到了两个预后生物标志物PYCR2和ADH1A：

上面两篇文章都是共识聚类应用非常经典也非常值得学习的文章，我把它们放在了资料下载区，有兴趣可以自取阅读。今天我们来看看如何利用ConsensusClusterPlus包进行共识性聚类分析。

共识聚类

共识聚类是一种为确定数据集中可能存在的亚类而提供定量分析结果的方法，可以使用芯片数据或者测序数据进行分析，在疾病中，这种方法可以用来发现了新的分子亚类，从而对疾病分类进行重新定义。接下来我们看看如何对测序数据进行共识性聚类，这里我们使用TCGA的COAD数据进行演示：

rm(list = ls())library(ConsensusClusterPlus)library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.3 v purrr 0.3.4## v tibble 3.1.2 v dplyr 1.0.6## v tidyr 1.1.3 v stringr 1.4.0## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --## x dplyr::filter() masks stats::filter()## x dplyr::lag() masks stats::lag()library(RColorBrewer)library(pheatmap)library(survival)library(survminer)## 载入需要的程辑包：ggpubrlibrary(Rtsne)library(ggplot2)library(ggsci)load("COAD_final.Rda") # 加载数据COAD_final[1:3, 1:3] # 观察数据## TCGA-AA-3518-11A-01R-1672-07 TCGA-A6-5662-11A-01R-1653-07## A1CF 11.923297 11.578310## A2M 15.231503 15.214207## A2ML1 6.083194 6.087948## TCGA-AZ-6603-11A-02R-1839-07## A1CF 10.601730## A2M 16.223713## A2ML1 6.232503

# 从数据上我们可以得知数据行名为基因名，列名为样本名，接下来我们需要提取特定基因集的表达量，这里我们选择m6A基因集

m6A_gene <- read_tsv("gene.txt")## ## -- Column specification --------------------------------------------------------## cols(## m6Agene = col_character()## )# 提取特定基因集的表达m6A_geneExp <- COAD_final[m6A_gene$m6Agene, ]m6A_geneExp <- as.matrix(m6A_geneExp)m6A_geneExp[1:3, 1:3]## TCGA-AA-3518-11A-01R-1672-07 TCGA-A6-5662-11A-01R-1653-07## METTL3 9.381633 9.708965## METTL14 10.592891 10.456011## WTAP 11.313789 11.617279## TCGA-AZ-6603-11A-02R-1839-07## METTL3 9.918261## METTL14 10.449458## WTAP 11.508632# 接下来我们进行共识聚类Cluster <- ConsensusClusterPlus(d = m6A_geneExp, # 分析矩阵 maxK = 7, # 最大聚类数目 reps = 1000, # 重抽样的次数 pItem = 0.8, # 样品的重抽样比例 clusterAlg = "pam", # 使用的聚类算法，可以选择"hc"(hclust), "pam", "km"(k-means) innerLinkage = "ward.D2", finalLinkage = "ward.D2",

 distance = "euclidean", # 计算距离的方法，可以选择pearson、spearman、euclidean、binary、maximum、canberra、minkowski

seed = 123456, # 设置随机种子，方便重复 plot = "png", # 结果图片的导出类型，可以选择"png"或者"pdf" title = "Consensus Cluster")## end fraction## clustered## clustered## clustered## clustered## clustered## clustered

结果解读

接下来我们简单看下结果图片：

tracking plot ：图片下方的黑色条纹表示样品（这里因为样本太多，看起来像是一个黑色图形），展示的是样品在k取不同的值时，归属的分类亚型情况，不同颜色的色块代表不同的亚型。取不同k值前后经常改变颜色分类的样品表示这个样品分类不稳定。若分类不稳定的样本太多，则说明该k值下的分类不稳定。

一致性累积分布函数图：简称为CDF图，这张图展示了k取不同数值时的累积分布函数，用于判断k值的最佳取值，当CDF达到一个近似最大值时，此时的聚类分析结果最可靠，对应的k值就是最佳k值，得到的分类就是最佳分类。

Delta Area Plot：这张图展示了 k 和 k-1 相比CDF曲线下面积的相对变化。当k = 2时，因为没有k=1，所以第一个点表示的是k=2时CDF曲线下总面积，而当k = 7时，曲线下面积趋近于平稳，所以7为最合适的k值。

提取数据重新画图

通过上面的结果解读我们可以选取出初步分类的亚组的个数，但是这种方法有时候带有主观性，我们也可以用PAC法确定最佳聚类数目。此外，如果对于Consensus Cluster Plus包得到的图片的颜色不太喜欢，我们还可以重新进行绘图：

# 用PAC的方法确定最佳聚类数Kvec <- 2:7x1 = 0.1x2 = 0.9 PAC <- rep(NA,length(Kvec)) names(PAC) <- paste("K =", Kvec, sep="") for(i in Kvec){ M = Cluster[[i]]$consensusMatrix Fn = ecdf(M[lower.tri(M)]) # 计算出共识矩阵 PAC[i-1] = Fn(x2) - Fn(x1)} optK <- Kvec[which.min(PAC)] # 理想的K值optK # 7## [1] 7# 通过PAC法，我们确定了k = 7是比较合适的聚类数目，接下来我们进行重新画图# 重新绘制热图annCol <- data.frame(Cluster = paste0("Cluster", Cluster[[7]][["consensusClass"]]), row.names = colnames(m6A_geneExp))head(annCol)## Cluster## TCGA-AA-3518-11A-01R-1672-07 Cluster1## TCGA-A6-5662-11A-01R-1653-07 Cluster1## TCGA-AZ-6603-11A-02R-1839-07 Cluster1## TCGA-AA-3496-11A-01R-1839-07 Cluster1## TCGA-AZ-6605-11A-01R-1839-07 Cluster1## TCGA-AA-3516-11A-01R-A32Z-07 Cluster1mycol <- brewer.pal(7, "Set3")annColors <- list(Cluster = c("Cluster1" = mycol[1],"Cluster2" = mycol[2],"Cluster3" = mycol[3],"Cluster4" = mycol[4],"Cluster5" = mycol[5],"Cluster6" = mycol[6],"Cluster7" = mycol[7]))heatdata <- Cluster[[7]][["consensusMatrix"]]dimnames(heatdata) <- list(colnames(m6A_geneExp), colnames(m6A_geneExp))heatdata[1:3, 1:3]## TCGA-AA-3518-11A-01R-1672-07## TCGA-AA-3518-11A-01R-1672-07 1.0000000## TCGA-A6-5662-11A-01R-1653-07 0.9984399## TCGA-AZ-6603-11A-02R-1839-07 0.8520249## TCGA-A6-5662-11A-01R-1653-07## TCGA-AA-3518-11A-01R-1672-07 0.9984399## TCGA-A6-5662-11A-01R-1653-07 1.0000000## TCGA-AZ-6603-11A-02R-1839-07 0.8495298## TCGA-AZ-6603-11A-02R-1839-07## TCGA-AA-3518-11A-01R-1672-07 0.8520249## TCGA-A6-5662-11A-01R-1653-07 0.8495298## TCGA-AZ-6603-11A-02R-1839-07 1.0000000# 绘制热图pheatmap(mat = heatdata, color = colorRampPalette((c("white", "steelblue")))(100), border_color = NA, annotation_col = annCol, annotation_colors = annColors, show_colnames = F, show_rownames = F)

dev.copy2pdf(file = "Cluster.pdf") # 保存图片## png ## 2

去提取分型结果

有了分型之后，我们需要把分型结果提取出来，这样可以探索不同亚型的分析结果，比如生存分析等等，接下来我们看看如何提取分型结果并进行后续分析：

# 取出分型结果Cluster_res <- annCol %>% as.data.frame() %>% rownames_to_column("patient_ID")table(Cluster_res$Cluster)## ## Cluster1 Cluster2 Cluster3 Cluster4 Cluster5 Cluster6 Cluster7 ## 47 79 75 85 87 48 50# 合并预后数据进行生存分析COAD_clinical <- read_tsv("COAD_clinical.txt")## ## -- Column specification --------------------------------------------------------## cols(## patient_ID = col_character(),## time = col_double(),## status = col_double()## )Cluster_res$patient_ID <- str_sub(Cluster_res$patient_ID, 1, 12)Cluster_res <- Cluster_res %>% inner_join(COAD_clinical, "patient_ID")head(Cluster_res)## patient_ID Cluster time status## 1 TCGA-AA-3518 Cluster1 2.0849315 0## 2 TCGA-A6-5662 Cluster1 3.9671233 0## 3 TCGA-AZ-6603 Cluster1 2.4630137 1## 4 TCGA-AA-3496 Cluster1 2.0849315 0## 5 TCGA-AZ-6605 Cluster1 0.4356164 1## 6 TCGA-AA-3516 Cluster1 1.0849315 1# 进行生存分析fit <- survfit(Surv(time, status) ~ Cluster, data = Cluster_res)ggsurvplot(fit, pval = TRUE, linetype = "solid", palette = mycol, surv.median.line = "hv", title = "Overall survival", ylab = "Cumulative survival (percentage)", xlab = " Time (Years)", legend.title = "Survival Plot", legend = c(0.85,0.70), risk.table = T, tables.height = 0.2, ggtheme = theme_bw(), tables.theme = theme_cleantable())

# 我们还可以对分型结果进行PCA分析和tSNE分析，PCA分析我们前面的推文已经讲过很多了，这里我们采用tSNE分析

tSNE_res<- Rtsne(t(m6A_geneExp), dims=2, perplexity=10, verbose=F, max_iter=500,check_duplicates=F)tsne<- data.frame(tSNE1 = tSNE_res[["Y"]][,1], tSNE2 = tSNE_res[["Y"]][,2], cluster = annCol$Cluster)ggplot(tsne,aes(x = tSNE1,y = tSNE2,color = cluster)) +geom_point(size = 4.5,alpha = 0.5)+scale_color_npg()+theme_classic()+stat_ellipse(level = 0.85, show.legend = F) + theme(legend.position = "top")

ggsave(filename = "tSNE.pdf")## Saving 10 x 8 in image

# 从tSNE的结果可以看出，亚型1、3、4区分度相对较好，后续分析就可以聚焦于这三种亚型进行分析，探索不同亚型的分子特征

好了，本期的内容到此结束！在本期中，我们一起探讨了共识聚类，包括代码实现、提取数据重绘图、生存分析和tSNE分析。其实怎么分型并不重要，重要的是分型结果的解读以及相关分子特征的探索，这也和我们发文章一样，分析方法来来去去就那些，就算换了一个可视化的图片，也还是那些内容，同样的分析方法，同样的内容，有的文章发了10+，有的却投3分也被拒。究其差异，在于思路和科研问题的不同。有一个好的idea才最重要，至于怎么实现，那是后面一步一步分析的事。

后台回复“feng 57”获取代码和数据，我们下期见！

学到了！别人发了CNS的数据，我还能发吗？教你一招，做出新意！（附代码）

单细胞分析专栏传送门

大年初一，单细胞分析从nature communication开始

8+单细胞文章从解读到套路总结全部送给你

30分多的单细胞分析文章是怎么做的？

不用R，点点点就能get的单细胞公共数据下载，你“学废”了吗

单细胞都有两个“对象”，你一个都没有

8成生信审稿人必问！单细胞第一张可视化图，你漏了吗？

我从Science中，“偷”到了这个实用的聚类分析技能！真舍不得分享！（附代码）

不骗你！学了这个R包，搞定九成高大上的单细胞分析！包治百病！一包封神！小白赶紧收起来！

Cell的生信分析能有多高级？单细胞研究必备分析！不会做么？简单，我来教你！（附代码）