読者です 読者をやめる 読者になる 読者になる

Data Science by R and Python

統計学を、広く、深く、わかりやすく。

ヒートマップで行列を可視化する! -ggplot2を用いたmicroarray発現量の可視化をしよう!-

今日のテーマは「ヒートマップ」

ヒートマップというのは、行列を可視化する時などに便利な手法です。有名なのは「遺伝子発現」というのを可視化する場合です。この他、店舗の一人の人が買った商品を可視化するのに使ってみてもおもしろいかもしれません。同じような購買パターンの人が視覚的にわかるようになります。一応、ヒートマップについてWeblioから引用しておきます。

ヒートマップ(英: Heat map)は個々の値のデータ行列を色として表現した可視化グラフの一種である。フラクタル図や樹形図は、変数によるヒエラルキー値を表現するため同様に色分ける事が多い。また、この用語は主題図の一種である階級区分図を示す。

用いるデータ

今回は、遺伝子発現データを用います。データは(ヒートマップ heatmap | R,マイクロアレイなどの遺伝子発現量の図示などによく使われる)で用いられているものをお借りしました。実際のデータを見ておくと、

> data <- read.table("http://stat.biopapyrus.net/data/arraydata.txt", header = TRUE)
> data <- as.matrix(data)
> head(data)
            AF       AX       AC       AT       AN       BF       BX       BC       BT       BN
gene1 -0.01183  1.37744  0.67014  0.82689  0.87908  0.88640  0.49236  0.88016  0.10693  0.38530
gene2 -0.00600  0.19820 -0.17458 -0.87388  0.19304 -0.13020 -1.56783  0.33063 -1.19051 -0.72618
gene3 -0.41678  0.52872  0.01384  0.08311  0.30768 -0.27890  0.48269  0.21181  0.37127 -0.32058
gene4 -0.69614  0.32774 -1.21520 -0.31868 -0.70187 -0.53080 -0.12315 -0.48518 -0.71129 -1.02066
gene5  0.91733 -1.23051 -0.32244  0.95149 -1.40095 -0.82365 -1.56783 -0.89023 -0.48405 -1.39034
gene6  1.46852  1.81342  3.12429  3.60012  1.93483  2.08088  1.63321  3.06953  2.32866  2.04482

さて、このデータをぱっと見ても、数値が並んでいるだけでわかりません。このデータの数値を可視化したいという課題に今回はアタックします。それがヒートマップです。実際、データ引用元でも普通にヒートマップをggplot2を使わないでもかけるよとあるんですが、まぁggplot2でも書いてみたいので、書いてみましょう。

heatmap by ggplot2

ggplot2でヒートマップを書くのは、そんなに単純ではありません。普通のheatmap関数を用いるときは、ただデータを引数に取ればいいんですが、ggplot2では関数が使えるようにデータを加工する必要があります。その加工はreshape2のmelt関数を用いて、データを以下のように展開するという操作に対応します。

> #library#
> library(ggplot2)
> library(reshape2)
> library(ggthemes)
> 
> #data arrange#
> data2=melt(data,id.var="AF")
> head(data2)
   Var1 Var2    value
1 gene1   AF -0.01183
2 gene2   AF -0.00600
3 gene3   AF -0.41678
4 gene4   AF -0.69614
5 gene5   AF  0.91733
6 gene6   AF  1.46852

このように展開できたら、ggplot2を用いてヒートマップを書くことができます。ここで、heatmapを書くための関数としてgeom_tilesを用いています。X軸とY軸はggplotの中のaesで指定して、中の値を何で埋める(fill)かを選ぶわけですが、そのとき上でいうvalueの値で埋めますという意味です。最後のscale_fill_gradientはヒートマップの色の指定です。普通に書いてしまうと、ヒートマップは「濃い青〜薄い青」という非常にみえにくい感じになるので、小さい値には”white"、大きな値は"red"で塗りつぶして!というのが最後の意味です。

> #heatmap#
> (p<-ggplot(data2,aes(as.factor(Var1),as.factor(Var2)))+geom_tile(aes(fill=value))+scale_fill_gradient(low="white",high="red"))

図を書いてみると、以下のようになります。
f:id:tomoshige_n:20140815234406p:plain

さて、しかしながら、X軸のラベルがみえにくいですよね。これを解消しましょう。それから、右側にある凡例も消してしまいましょう。余白も小さくして。そのためには以下のようにコマンドを書きます。

> #-heatmap-#
> p + labs(x = "",y = "") + scale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) +  theme(legend.position = "none",axis.ticks = element_blank(), axis.text.x = element_text(size = 7, angle = 280, hjust = 0, colour = "grey50"))


ここで、余白の指定はscale_x_discrete(expand = c(0, 0)) + scale_y_discrete(expand = c(0, 0)) で行っています。これは余白を0にせよという意味です。次のthemeでlegend.position="none"と書いているのが凡例を消す操作。axis.ticks=element_blank()は軸を消す操作。axis.text.x = element_text(size = 7, angle = 280, hjust = 0, colour = "grey50"))でラベルに表示する文字の大きさ、その回転角などを指定しています。このようにして書くと、図は...

f:id:tomoshige_n:20140815234845p:plain

とそこそこ見やすくなりました。以上でggplot2でheatmapを書けたことになります。ただ、もうちょっとわかりやすくしたいんですけどね。実際にRに入ってるheatmap関数を使った図を見てみると、あー、こっちの方がいいなと思います。

heatmap関数によるヒートマップ作成

以下のようにコマンドを打ちます。すると、ヒートマップが作成されます。簡単です。しかも、クラスタリングした結果まで書くことができます。これggplot2で書こうと思うときっと大変なんだと思います...

> heatmap(
+   data,                     # 行列型データ
+   main = "Heat colors",     # グラフタイトル
+   Rowv = TRUE,              # 列にクラスタ図を描く
+   Colv = TRUE,                # 行にクラスタ図を描かない
+   distfun = dist,           # distance calculation
+   hclustfun = hclust,       # clustering method
+   col = heat.colors(256),    # 色彩の指定
+ )

f:id:tomoshige_n:20140815235044p:plain

まとめ

今回は、ヒートマップを用いてデータを可視化してみました。最初にもいいましたが、これでマーケティングデータを可視化することができれば、ユーザーが共通して買っているものや、この人は買っているけど、この人は買っていないなどわかりそうですね。ちなみに、エクセルでもこれはできるようです...やる気にはなりませんけどwww

追記...こんなこともできます。

ヒートマップの中に値を書き込みたいなという場合には、次のようにすることで図を書くことができます。ただし、行・列数が多い場合にはきつそうです。。。

> dat <- matrix(rnorm(100, 3, 1), ncol=10)
> names(dat) <- paste("X", 1:10)
> dat2 <- melt(dat, id.var = "X1")
> ggplot(dat2, aes(as.factor(Var1), Var2, group=Var2)) +
+     geom_tile(aes(fill = value)) + 
+     geom_text(aes(fill = dat2$value, label = round(dat2$value, 1))) +
+     scale_fill_gradient(low = "white", high = "red") 

f:id:tomoshige_n:20140815235501j:plain