R语言实现双聚类

聚类大家应该都听说过，但是双聚类想必大家接触的比较少，今天我们就给大家介绍下双聚类，首先看下基础的定义：针对二维数据进行处理的算法。假设给定矩阵M，寻找到矩阵M的多个子矩阵A，对于每一个A满足其指定条件进行聚类，最后得到需要的子矩阵B。目前广泛的模型有四种：矩阵等值模型、矩阵加法模型、矩阵乘法模型和信息共演变模型。接下来在R语言中的实现需要用到包biclust。其安装如下：

install.packages(“biclust”)

然后，我们来看下此包的使用。首先是数据的预处理函数包括discretize、binarize。

1. Discretize主要是将表达矩阵转化为间隔等级的离散矩阵，其中的主要参数是quant主要是用来确定是否用分位数来划分等级。Nof来确定等级的水平。实例：

library(biclust)
data(BicatYeast)

discretize(BicatYeast[1:10,1:10])

2. Binarize 转化为二元的矩阵。通过设定threshold来判断0，1的分割界限，如果为空默认中位数。另外还有函数binarizeByPercentage通过输入percentage来生成矩阵，反应二元矩阵中0，1的密度；densityOnes 返回1对0的百分比。

接下来我们看下此包的主函数biclust(x,method,...)。这里主要是需要选择自己需要的内置算法。其算法包括BCCC, BCXmotifs, BCPlaid, BCSpectral,BCBimax, BCQuest。

BCCC（Cheng, Y.& Church, G.M. Biclustering of Expression Data Proceedings of the EighthInternational Conference on Intelligent Systems for Molecular Biology, 2000, 1,93-103）

聚类定义为均方残差得分小于某个事先指定的阙值δ的子矩阵

BCXmotifs (Murali, T. & Kasif, S.Extracting Conserved Gene Expression Motifs from Gene Expression Data PacificSymposium on Biocomputing, sullivan.bu.edu, 2003, 8, 77-88)

该算法需要离散的矩阵作为输入。被定义为在几列样本中，基因表达几乎不变的行。

BCPlaid (Heather Turner et al,"Improved biclustering of microarray data demonstrated through systematicperformance tests",Computational Statistics and Data Analysis, 2003, vol.48, pages 235-254.)

该算法是对原始的格子数据模型的改进。该算法将数据矩阵模型化为一组层，通过最小化误差来拟合数据。

BCSpectral（Klugeret al., "Spectral Biclustering of Microarray Data: Coclustering Genes andConditions", Genome Research, 2003, vol. 13, pages 703-716）

该算法找到的值高于相应的其他行和列中的值。每行和每列只属于一个双聚类，因此重新排列行和列中的这些高值，使这些分区沿着对角线连续显示。

BCBimax (Prelic, A.; Bleuler, S.;Zimmermann, P.; Wil, A.; Buhlmann, P.; Gruissem, W.; Hennig, L.; Thiele, L.& Zitzler, E. A Systematic Comparison and Evaluation of BiclusteringMethods for Gene Expression Data Bioinformatics, Oxford Univ Press, 2006, 22,1122-1129)

假设每个基因只有两个可能的表达水平：随某一实验条件变化（1）或不变（0）。聚类为元素全为1或0的子矩阵。把原矩阵分为子矩阵U、V和0矩阵：

BCQuest (Murali, T. & Kasif, S.Extracting Conserved Gene Expression Motifs from Gene Expression Data PacificSymposium on Biocomputing, sullivan.bu.edu, 2003, 8, 77-88)

搜索具有相同或类似答案的问题的子组。

每个算法都有自己对应的参数设置，我们就不一一展开了，毕竟有些参数咱也不知道咋设，但是应该默认基本上符合大部分情况的。我们直接看下实例：

est <- matrix(rbinom(400,50, 0.4), 20, 20)
bics1<- biclust(test, method=BCCC(),delta=1.5,  alpha=1, number=10)

最后我们看下数据结果的可视化，我们直接看下实例：

s2=matrix(rnorm(5000),100,50)
drawHeatmap(s2)

  s2[11:20,11:20]=rnorm(100,3,0.3)
 set.seed(1)
 bics <- biclust(s2,BCPlaid(), back.fit = 2, shuffle = 3, fit.model =~m + a + b,
 iter.startup = 5, iter.layer = 30, verbose = TRUE)
 parallelCoordinates(x=s2,bicResult=bics,number=1, plotBoth=TRUE,
plotcol=TRUE, compare=TRUE,info=TRUE,bothlab=c("Genes Bicluster
1","Conditions Bicluster1"), order =TRUE)

上面的图称为平行坐标图，其中图中每个X坐标代表一个特征，上面的Y值代表该特征的值，每个样本表示出来就是一个贯穿的折线图。一般来说，用不同的颜色代表不同的类别，这样可以方便的看出不同特征对分类的影响。

有以下两个原则：

主要看相同颜色的折线是否集中，若在某个属性上相同颜色折线较为集中，不同颜色有一定的间距，则说明该属性对于预测标签类别有较大的帮助。

若某个属性上线条混乱，颜色混杂，则较大可能该属性对于标签类别判定没有价值。

另外一个可视化就是气泡图，我们直接看实例：

data(BicatYeast)
 set.seed(1)
 bics1 <- biclust(BicatYeast,BCPlaid(), back.fit = 2, shuffle = 3,fit.model = ~m + a + b,
 row.release = 0.7, col.release = 0.7,
 verbose = FALSE, max.layers = 10, iter.startup = 5,
 iter.layer = 30)
 bubbleplot(BicatYeast,bics1, showLabels=TRUE)

loma=binarize(BicatYeast,2)
bics2=biclust(loma,BCBimax(), minr=4,minc=4, number=10)
bubbleplot(BicatYeast,bics1,bics2, showLabels=TRUE)

上图为气泡图，其中每个双光泽都表示为一个圆形(气泡)。颜色表示双光泽所属的双光泽集(最多可同时表示三个双光泽集)。亮度表示双光泽均匀性(较暗、较不均匀)。大小表示双光泽的大小，如(基因数量)x(条件数量)。定位是基因和条件的二维投影。