R：ヒートマップを作成する方法

ヒートマップを作成するコードの一つを、忘備録として記す。

いつもながら、working directoryの場所にご注意いただきたい。今回の場合は、CELファイルが置いてあるファイルに設定している。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

library(affy) cel_dat <- ReadAffy() eset <- rma(cel_dat) emat <- exprs(eset) max <- apply(emat, 1, max) min <- apply(emat, 1, min) sub <- max - min emat_slct <- cbind.data.frame(emat, sub) emat_slct <- emat_slct[which(emat_slct$sub > 1),] emat_slct$sub <- NULL library(stringr) array_names <- colnames(emat_slct) array_names <- str_sub(array_names, end = -5) colnames(emat_slct) <- array_names library(gplots) heatmap.2(as.matrix(emat_slct), col=bluered(50), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, cexCol = 0.5) # related to font size

1~5までは、CELファイルの中身を取り出す作業。emat (expression matrix) には、CELファイルから抽出された遺伝子情報がlog2変換されて入っている。

7~9は、各遺伝子プローブの最大値から最小値を引いた値を求め、subに格納している。

11~13では、subをematに結合した後、subが1以下(発現変動が極めて小さい)であるものを除き、subを削除している。

15~18では、サンプル名から".CEL"を除き、見やすくしている。

20~23でヒートマップを描く作業になる。ヒートマップを描くための関数・パッケージは様々にあるが、今回はgplotsのheatmap.2を用いた。

Heatmaplyも試したが、遺伝子プローブの数が多い場合には計算に時間がかかりそうで断念した。

21行目をscale="none"とすると、発現量(log2)がそのまま表示されるようになる。scale="row"とすることで、サンプル間の相対量(z-score)に変えることができる。

他にも良い関数やパッケージがあれば、ぜひご教授ください。