ヒートマップ image

行列型のデータをグラフに描く関数として image 関数がある。ある区域における草丈の長さやある区域における海抜を記録したデータなどが行列型になる。(マイクロアレイのデータをヒートマップに描く場合は、一般に heatmap 関数を利用する。)

以下は 100 行 4 列のサンプルデータを利用して image 関数でプロットする例である。

# サンプルデータの生成
x <- cbind(
       rnorm(100, 10, 3),
       rnorm(100, 20, 4),
       rnorm(100, 30, 5),
       rnorm(100, 40, 6)
     )
colnames(x) <- c("A", "B", "C", "D")   # カラム名を付ける

# サンプルデータの確認
head(x)
##              A        B        C        D
## [1,]  7.888219 14.10776 22.53229 32.95763
## [2,] 13.206582 26.30063 28.77801 44.26600
## [3,] 13.088457 21.53165 35.23720 36.51519
## [4,] 10.892165 17.07894 30.14704 28.69056
## [5,] 11.595863 14.04913 32.43846 48.71336
## [6,] 12.903338 23.65794 24.89183 44.12978

ヒートマップは image 関数をそのまま利用して描いてもよいが、ここではスケールバーを付けるため、プロットする色と位置を最初に調整する。

# 画面を 4:1 で分割
layout(matrix(data = c(1, 2), nrow = 1, ncol = 2),
       widths = c(4, 1), heights = c(1, 1))


# カラースケールを作成する
max.value <- ceiling(max(x))    # 最大値(小数の場合は切り上げて整数にする)
min.value <- ceiling(min(x))    # 最小値(小数の場合は切り捨てて整数にする)
colorRamp <- rgb(
  seq(1, 1, length=256),   # 赤成分
  seq(1, 0, length=256),   # 緑成分
  seq(1, 1, length=256)    # 青成分
)
colorLevels <- seq(min.value, max.value, length = length(colorRamp))


# ヒートマップを描きます
par(mar = c(5.5, 4.5, 2.5, 2))     # 余白調整
image(1:ncol(x), 1:nrow(x), t(x[rev(1:nrow(x)), ]),
  col = colorRamp,
  axes = FALSE,
  main = "Expression Levels",
  xlab = "genes",
  ylab = "samples",
  zlim = c(min.value, max.value)
)
axis(1, at = 1:ncol(x), labels = colnames(x))
box()                              # ヒートマップの周りの枠線を描く


# カラースケールを描く
par(mar = c(5.5, 2.5, 2.5, 2))
image(1, colorLevels,
      matrix(colorLevels, ncol = length(colorLevels), nrow = 1),
      col = colorRamp, xlab="", ylab="", xaxt="n")
box()
layout(1)
image関数で描いたヒートマップ関数

比率をヒートマップ描く

比率をヒートマップを描く例。0 から 1 未満の値をシアン色で、1 よりも大きい値をマゼンタ色で描き、1 に近づくと白くなるようなグラデーションでヒートマップを描く。

# サンプルデータの生成
y <- cbind(
         abs(rnorm(100, 3, 2)),
         abs(rnorm(100, 5, 4)),
         1 / abs(rnorm(100, 1, 4)),
         abs(rnorm(100, 3, 2))
     )
colnames(y) <- c("A", "B", "C", "D")


# 画面を分割
layout(matrix(data = c(1, 2), nrow = 1, ncol = 2), widths = c(4, 1), heights=c(1, 1))

# 最大値と最小値を求めてカラーパレットを作成
x <- y
maxlevel <- ceiling(max(x))       # 最大値
minlevel <- ceiling(1 / min(x))   # 最小値の逆数(比率を扱うためにとりあえず変換)
x[x < 1] <- - 1 / x[x < 1] + 2
if (min(x) >= 1) {
  colorRamp <- c("#FFFFFFFF", cm.colors((maxlevel - 1) * 32)[((maxlevel - 1) * 16 + 1):((maxlevel - 1) * 32)])
} else if (max(x) <= 1) {
  colorRamp <- c(cm.colors((minlevel - 1) * 32)[1:((minlevel - 1) * 16 - 1)], "#FFFFFFFF")
} else {
  colorRamp <- c(cm.colors((minlevel - 1) * 32)[1:((minlevel - 1) * 16 - 1)], "#FFFFFFFF",
                 cm.colors((maxlevel - 1) * 32)[((maxlevel - 1) * 16 + 1):((maxlevel - 1) * 32)])
}
# カラーパレット
colorLevels <- seq(2 - minlevel, maxlevel, length = length(colorRamp))

# ヒートマップ
par(mar = c(5.6, 4.5, 2.5, 2))
image(1:ncol(x), 1:nrow(x), t(x[rev(1:nrow(x)), ]),
      col = colorRamp, ylab = "samples", xlab = "genes", main = "", axes = FALSE,
      zlim = range(2 - minlevel, maxlevel))
axis(1, at = 1:ncol(x), labels = colnames(x), cex.axis = 0.8)
box()

# カラーパレット
par(mar = c(5.6, 2.5, 2.5, 2))
image(1, 0:length(colorRamp),
      matrix(colorLevels, ncol = length(colorRamp), nrow = 1),
      col = colorRamp, xlab = "", ylab = "",
      xaxt = "n", yaxt="n")
lb <- seq(from = - minlevel + 2, to = maxlevel, by = 1)
lc <- lb
lc[lc < 1] <- 1 / (2 - lc[lc < 1])
axis(2,
     at = seq(from = 0, to = length(colorRamp),
     by = length(colorRamp) / (length(lb) - 1)),
     labels =  c(rev(paste("1/", 1:minlevel, sep = "")[-1]), sprintf("%d", 1:maxlevel)),
     cex.axis = 0.8)
box()
layout(1)
比率のヒートマップ