如何通过一致性聚类实现对表达谱数据的亚型分类

image

​通过ConsensusClusterPlus包对基因表达谱执行一致性聚类（Consensus Clustering）

在大样本的组学分析中，经常需要讨论样本的分子分型。文章中最常见的方法，就是使用一种叫一致性聚类（Consensus Clustering）的方法，对转录组、蛋白组谱等数据进行聚类，最后可以将样本划分为不同的聚类群，不同聚类群的样本之间在转录组、蛋白组等分子模式上存在明显差异，但各聚类群内的样本的分子模式则较为相似，如此则实现了对大样本队列的分子分型的目的。

例如，在“Proteogenomic landscape of squamous cell lung cancer”这篇文献中，作者基于108个肺鳞状细胞癌样本的定量蛋白组数据，通过一致性聚类将108个肿瘤组织划分为5种分子亚型，即（1）炎症亚型A、（2）炎症亚型B、（3）氧化还原亚型A、（4）氧化还原亚型B以及（5）混合亚型。并在获得了分子亚型后，后续展开的详细亚型特征的描述和讨论。

image

好了，现在假设我们也有类似的一批组学样本队列，我们也期望通过一致性聚类实现对转录组、蛋白组或者蛋白修饰谱的分型，该如何实现呢？本节教程简介通过R语言ConsensusClusterPlus包实现一致性聚类的方法。

1 准备表达谱数据集

首先准备一个基因表达矩阵，将它读入到R中。可以是转录组（如RNA-seq，芯片数据等），也可以是定量蛋白组，或者蛋白磷酸化、糖基化等修饰。具体以哪种类型的数据为主，根据实际关注的问题来。如果您更期望使用转录组进行分型，就使用RNA表达谱；如果您更期望使用蛋白组进行分型，就使用蛋白定量谱；如果您更关注蛋白修饰的分型，就使用表观修饰的组学。

本次我们以Biobase包的芯片数据集为例，展示如何对基因表达谱执行一致性聚类分析。首先来看一下示例数据，该数据集一共包含26个样本，500个基因的表达谱。

#以Biobase包的表达谱芯片数据集为例
library(Biobase)

data(geneData)
geneData[1:6,1:6]

image

2 通过ConsensusClusterPlus包执行聚类

能够执行一致性聚类的R包很多，但基本原理都是差不多的。这里我们以ConsensusClusterPlus包的方法为例作为展示。

#通过Bioconductor安装ConsensusClusterPlus包
#BiocManager::install("ConsensusClusterPlus")

#加载ConsensusClusterPlus包
library(ConsensusClusterPlus)

#对基因表达数据执行中位数中心化，以用于后续聚类
dc <- sweep(geneData,1, apply(geneData,1,median))

#使用函数ConsensusClusterPlus()执行聚类
cluster <- ConsensusClusterPlus(
  d = dc, 
  maxK = 4, 
  pItem = 0.8, 
  pFeature = 1, 
  clusterAlg = "hc", 
  distance = "pearson", 
  seed = 1234, 
  innerLinkage = "complete", 
  finalLinkage = "complete", 
  corUse = "pairwise.complete.obs",
  plot = 'png', 
  title = "example"
)

过程大致包括两步，（1）对表达谱数据的标准化，（2）对标准化后的数据执行聚类。函数执行完毕后，R语言的当前工作目录下会生成一个文件夹（本示例以“example”未命名），里面存放了聚类结果图。

3 结果说明

在本示例中，我们通过参数“maxK = 4”指定了期望划分的最大聚类群数量，即从聚2类开始，逐一增加类别数量进行尝试，直到到达指定的最大聚类群数量4为止结束。

结果中，这两张图可以帮助我们评估选择最合适的聚类数量。我们需要同时考虑这两个标准，根据左图选择CDF下降坡度更小的曲线，根据右图选择CDF更高的值。通常而言这两个标准趋势相反，也就是我们需要尽量保证CDF下降不那么剧烈且CDF值也不能太小。在这里，似乎当k=3时是最合适的。当然，有时不一定要遵守这个方法，一些研究也按照自己的研究目的选择了其它最优k值标准。

image

这些图展示了在不同的k值（划分的聚类群数量）下，样本的聚类情况。在示例中，似乎k=3时相对更好一些。

image

该图展示了在不同的k值（划分的聚类群数量）下，各样本所属的聚类群。

#查看样本所属的聚类群
cluster

#例如，查看当划分3个聚类群时，各样本所属的聚类情况
#其它聚类群的情况，更改数字即可
cluster[[3]]$consensusClass 

image