いろいろ倉庫

KNIME、EXCEL、R、Pythonなどの備忘録

【R】Heatmap.2のクラスタリング情報を扱ってみたい

・お題:遺伝子発現を比較する際に、ヒートマップというカラフルな図をheatmap.2を使って作ってみた。この図では階層型クラスタリングとセットで使われることが多い。図示するだけではなく、階層型クラスタリングの情報も扱ってみたい。

 

・正確には、クラスタリングしてから、クラスタリング情報をもとにheatmapを描いてもらう感じ。とりあえず、ライブラリを読み込んでデータセットを作る。

> library(tidyverse) #Rを便利に使うパッケージ
> library(gplots) #heatmap.2を含むパッケージ

> df <- data.frame(
+   sample1 = c(1, 2, 3, 4, 5),
+   sample2 = c(1.5, 2.5, 3.5, 4.5, 5.5),
+   sample3 = c(2, 3, 4, 5, 6),
+   sample4 = c(5, 4, 3, 2, 1),
+   sample5 = c(5.5, 4.5, 3.5, 2.5, 1.5),
+   sample6 = c(6, 5, 4, 3, 2),
+   row.names = c("geneA", "geneB", "geneC", "geneD", "geneE")
+   )

> df

 

・サンプル構成情報と処置条件情報もつけておく。

> samples <- names(df)
> groups <- 
+   c("DMSO","DMSO","DMSO","TEST","TEST","TEST") %>% 
+   as.factor(.)

 

・図示するだけなら、heatmap.2だけでできる。

> col <- colorpanel(1000,"blue","white","red") #色を指定しておく。

> col

以下略

> heatmap.2(x = as.matrix(df),
+           scale = "row",
+           labRow = row.names(df),
+           ColSideColors = c("green","orange")[groups],
+           labCol = samples, 
+           col = mycol,
+           trace = "none",
+           density.info = "none", 
+           margin = c(8,6),
+           dendrogram = "row")

クラスタリングの情報を使いたい場合、階層型クラスタリングを別にやる。

> hcl <- hclust(dist(df))

> plot(hcl) #図示してみる


・このクラスタリング情報を使ってヒートマップを描くには、Heatmap.2関数に引数Rowvにas.dendrogram(hcl)を与えれば良い。今回は行なのでRowvだけれど、列ならColvを与えてやれば良い。

> heatmap.2(x = as.matrix(df),
+           scale = "row",
+           labRow = row.names(df),
+           ColSideColors = c("green","orange")[groups],
+           labCol = samples, 
+           Rowv = as.dendrogram(hcl),
+           col = mycol,
+           trace = "none",
+           density.info = "none", 
+           margin = c(8,6),
+           dendrogram = "row")

 

・階層型クラスタリングなので、切るHeightによってできるクラスタの数が変わる。例えば、Heightを6や8にすれば、geneA, B / gene E, C, Dの2つのクラスタに分けられる。

> cutree(hcl, h=6)

・Heightを3や4に設定すれば、geneA, B / geneE / geneC,Dの3つのクラスタに分けられる。

> cutree(hcl, h=3)

 

クラスタ数を指定したいなら、引数kを与える。

> cutree(hcl, k=2)

 

・これはベクトルなので、適当に抜き出せる。

> cutree(hcl, h=3) %>% is.vector(.)
[1] TRUE

> row.names(df)[cutree(hcl, h=3) == 1]
[1] "geneA" "geneB"

 

・hclust関数の引数methodにいろいろ与えると、他のアルゴリズムでの階層型クラスタリングもできる。

 

 

おわり