2016年2月18日木曜日

Rで統計:階層クラスター分析(HCA)

どうも.中谷航太です.
備忘録がてら勉強したRの解析方法を記録しておきます.
今回は,Rでの階層クラスター分析付きヒートマップの作り方について.

①まずデータの準備から始めます.


  Wild Type strain1 strain2 strain3 strain4
metabolite1 7234526 6255017 6754156 490024.2 7127350
metabolite2 1970236 1685756 1803086 148639.6 1199489
metabolite3 18111892 16226529 17646758 1254115 18497672
metabolite4 266696 151436 206675.3 57778.57 384859
metabolite5 12992368 11370112 12383234 883364.6 13461877
・・・・・・・・・
・・・・・・・・・
・・・・・・・・・

みたいなデータをエクセルで作って,csvファイルで保存します.
要は分析結果をエクセルにまとめればいいだけです.

今回はデータをわざわざ新しくデータを用意するのも面倒なので,

a1 <- rnorm(100,mean=5000, sd=500)
a2 <- rnorm(100,mean=4000, sd=400)
a3 <- rnorm(100,mean=3000, sd=300)
a4 <- rnorm(100,mean=2000, sd=200)
a5 <- rnorm(100,mean=1000, sd=100)
x <- matrix(c(a1, a2, a3, a4, a5), nrow=5, ncol=100)
colnames(x) <- c(1:100)
rownames(x) <- c("WT", "strain1", "strain2", "strain3", "strain4")

これでデータ行列を作ります.
WT, strain1~4の代謝物を分析して,それぞれ代謝物1~100の分析値をエクセルにまとめたものと考えてもらえばいいです.


②次に,類似度 or 距離の選択方法を指定し,計算します.

階層クラスター分析は,サンプルと代謝物の両方で行います.
不要なら要らない方を削除してください.

d1 <- dist(x, method="euclidean")
d2 <- dist(t(x), method="euclidean")
# methodは "euclidean", "maximum", "manhattan", "canberra", "binary" or "minkowski"から選んで下さい.距離の計算方法が違うので,それぞれ結果が少し異なります.

③次にクラスター合併のアルゴリズムを指定し,計算します.


c1 <- hclust(d1, method="ward.D")
c2 <- hclust(d2, method="ward.D")
# methodは"ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC)から選んで下さい.

これで階層クラスタ解析の計算は終わったので,

④最後に結果をヒートマップで表します.


gplotsというパッケージを使うので,初回の人は以下を実行.

install.packages("gplots")

参照するパッケージを指定してヒートマップを作ります.

library(gplots)
hm <- heatmap.2(as.matrix(x),
                col=colorpanel(256, low = "yellow", mid="black", high= "blue"),
                scale="none",
                margins=c(6,30),
                Rowv=as.dendrogram(c1),
                Colv=as.dendrogram(c2),
                density.info="none",
                trace ="none",
                key ="True")

図の縦横比はmarginsで,色のグラジエントはcolで弄れます.
scaleは予めエクセルにまとめる段階で自分で計算しておくのが良い気がします.
今回はデータを適当に作ったので"column"にしています.

以上を実行するだけでHCA付きヒートマップがゲットできました!
お疲れ様です.


最後にコードだけでまとめておきます.
a1 <- rnorm(100,mean=5000, sd=500)
a2 <- rnorm(100,mean=4000, sd=400)
a3 <- rnorm(100,mean=3000, sd=300)
a4 <- rnorm(100,mean=2000, sd=200)
a5 <- rnorm(100,mean=1000, sd=100)
x <- matrix(c(a1, a2, a3, a4, a5), nrow=5, ncol=100)
colnames(x) <- c(1:100)
rownames(x) <- c("WT", "strain1", "strain2", "strain3", "strain4")
d1 <- dist(x, method="euclidean")
d2 <- dist(t(x), method="euclidean")
c1 <- hclust(d1, method="ward.D")
c2 <- hclust(d2, method="ward.D")
library(gplots)
hm <- heatmap.2(as.matrix(x),
                col=colorpanel(256, low = "yellow", mid="black", high= "blue"),
                scale="column",
                margins=c(6,30),
                Rowv=as.dendrogram(c1),
                Colv=as.dendrogram(c2),
                density.info="none",
                trace ="none",
                key ="True")

読み込みからヒートマップ化する場合はこっちで.
setwd("「データのある場所」")
x <- read.csv("「データの名前」", header=T, row.names = 1)
d1 <- dist(x, method="euclidean")
d2 <- dist(t(x), method="euclidean")
c1 <- hclust(d1, method="ward.D")
c2 <- hclust(d2, method="ward.D")
library(gplots)
hm <- heatmap.2(as.matrix(x),
                col=colorpanel(256, low = "yellow", mid="black", high= "blue"),
                scale="none",
                margins=c(6,30),
                Rowv=as.dendrogram(c1),
                Colv=as.dendrogram(c2),
                density.info="none",
                trace ="none",
                key ="True")

感想:楽しかった.