【Plot】[pheatmap] 外れ値があっても色割り当てを工夫してヒートマップを見やすくする

https://zenn.dev/rchiji/articles/451d2458877665

簡単にヒートマップを描くことができるpheatmapだが、外れ値があると色が偏ってしまい、急に表現力が落ちてしまう。
ここでは色割り当てを工夫するtipsを紹介する。

ライブラリ読み込み

install.packages("pheatmap")
library(pheatmap)

デモデータ

irisの数字列のデータのみを使用。irisデータの5列目は文字列で、文字列が残っているとpheatmapに使用でない。

demo <- iris[,1:4]
pheatmap(demo, show_rownames = F)

【Scaling】

irisデータをそのままヒートマップに使用すると、各項目でスケール感が違うものを表示しているので、項目ごとに色が偏ってしまっており、サンプル間の差ではなく、項目間の差が目立ってしまう。
項目内でどのサンプルが平均より離れているのかを見やすくするには、scaleオプションが有用。

pheatmap(demo, show_rownames = F, scale = "column")

項目内での平均値を0, 分散を1に標準化して表示しているのでデータの分布が確認しやすくなる。
行方向にscaleしたい場合はscale = "row"。

最大値を1にする

各列の最大値で各列の値を割れば列の最大値が1になるように補正できる。

hoge <- apply(X = demo, MARGIN = 2, FUN = function(x){
     x/max(x)
     })
pheatmap(hoge, show_rownames = F)     

0~1の範囲にする

最小値を引いてから、上記の処理をすると、0~1の範囲に補正できる。

hoge <- apply(X = demo, MARGIN = 2, FUN = function(x){
　　　x <- x - min(x)
   　 x/max(x)
      # (x - min(x))/(max(x) - min(x)) でも同じ
     })
pheatmap(hoge, show_rownames = F)     

値をclipする

特定の値の範囲に制限する。範囲外の値は強制的に範囲の最小値/最大値に上書きされる。

# scaling。中心が0。SDが1
hoge <- scale(demo)
# 値が-2 ~ 2になるようにclipする
hoge <- pmax(pmin(hoge, 2),-2)

pheatmap(hoge, show_rownames = F) 

※ こういうこともできるというだけで推奨しているわけではない。

【plotの色】

pheatmapのデフォルト色にはRColorBrewerのRdYlBuパレットが使用されている。このパレットから100色分の色コードを作って値に色を割り当てている。

このパレットは本来100色も無いのだが、colorRampPalette()機能により100色に補完している。

RdYlBuパレットの7色を100色に補完

colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))(100)

100個の色コード

このようにcolorRampPalette()機能を使えば簡単にグラデーションを作ることができる。

pheatmap(demo, show_rownames = F, scale = "column",
	 color = colorRampPalette(c("blue","white","red"))(n=100))

色を減らすと最小値から最大値までに割り当てられる色の数が減るため、表現力が低下することがある。

pheatmap(demo, show_rownames = F, scale = "column",
	 color = colorRampPalette(c("blue","white","red"))(n=10))

10色のみ与えた場合

ただ、color引数はデフォルトの100色より増やしても特に表現力の向上は無いように思われる。

【breaks引数で色を割り当てる範囲を決める】

外れ値を含むような極端な例で見てみる。

demo2 <- demo^5
pheatmap(demo2, show_rownames = F)

0~30000の範囲を100色で表現しており、値が高い箇所に色が引っ張られてる。

ここまで極端な場合、scale引数を使ってもうまく表現することができなくなる。（heatmapの下半分はサンプル間の違いが見にくい。）

pheatmap(demo2, show_rownames = F, scale = "column")

対策として、ある一定の値を頭打ちになるような色割り当てを行う。

試しに-2から2の間に色を割り当ててみると、見易さが大きく改善できるのがわかる。

pheatmap(demo2, show_rownames = F, scale = "column",
	 breaks = seq(-2,2,length.out=100))

※ 上述の値のclipと同様に見えるが、clipしてからクラスタリングと元の値のクラスタリングは結果が異なる。ここでは元の値でクラスタリングして色割り当てだけを工夫している。

seq関数はデフォルトでは第一引数から第二引数まで1ずつ値を増減させたベクトルを作成できる。その時にlength.outオプションを使用すると、第一引数から第二引数までを指定した要素数に分割したベクトルを作成できる。

pheatmapのbreaks引数に数字100個のベクトルを指定すると、100色の各色がどの値に割り当てられるのか指定することができる。

ちなみに、color引数の色数とbreaks引数の値の数は揃える必要がある。

50色指定なら、breaksも50

pheatmap(demo2, show_rownames = F, scale = "column",
         breaks = seq(-2,2,length.out=50),
         color = colorRampPalette(c("blue","white","red"))(50))

色は100色なのに、breaksが50。

pheatmap(demo2, show_rownames = F, scale = "column",
         breaks = seq(-2,2,length.out=50),
         color = colorRampPalette(c("blue","white","red"))(100))

最初の50色しか使用されない

色は50色なのに、breaksが100

pheatmap(demo2, show_rownames = F, scale = "column",
         breaks = seq(-2,2,length.out=100),
         color = colorRampPalette(c("blue","white","red"))(50))

色が足りず灰色で表示される

【色ベクトルを途中から同色にする】

上記のbreaksの例だと、legendが真の最小値/最大値を示さなくなり、データの分布の情報を失ってしまう。

例えば、100色あるとして、1-50色までグラデーションをつけて、51-100色目までを50色目と同じにすることでlegendの値の範囲を変えずに頭打ちを表現することができる。

# 50色のグラデーション色ベクトルを作成
color <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))(50)
# 50色目を50回繰り返して、グラデーション色ベクトルの後ろにつなげる
color <- c(color, rep(color[50], times=50))

51色目から同じ色コードの色ベクトル

デモデータをそのまま表示した場合

demo3 <- demo^3
pheatmap(demo3, show_rownames = F)

中央値以降を同じ色で表示

pheatmap(demo3, show_rownames = F, color = color)

このようにしてやれば第一四分位数以下、第三四分位数以上の値に色を引っ張られないようにすることもできる。

# 50色のグラデーション色ベクトルを作成
color <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))(50)
# 1色目*25回 + 50色グラデーション + 50色目*25回
color <- c(rep(color[1], times=25), color, rep(color[50], times=25))

pheatmap(demo2, show_rownames = F, scale = "column",
          color = color)

外れ値がわずかにしかない場合は、95色ぐらいまでグラデーションで96-100色を同色にするでもいいかもしれない。

【分位数を使った色割り当て】

quantile関数はデフォルトでデータの最小値、第一四分位数、中央値、第三四分位数、最大値を返す。

quantile関数

quantile(as.matrix(demo2)) # data.frameだとエラーが出たのでmatrixに変換している。

0% 25% 50% 75% 100%
0.00001 14.19857 335.54432 3450.25251 30770.56399

probs引数で分位を指定することができる。

probs引数

quantile(as.matrix(demo2), probs = c(0,0.33,0.66,1))

0% 33% 66% 100%
0.00001 64.36343 2059.62976 30770.56399

seqとlength.outを使えば欲しい数だけの分位を作成できる。

n = 10
quantile(as.matrix(demo2), probs = seq(0,1,length.out=n))

0% 11.11111% 22.22222% 33.33333% 44.44444% 55.55556% 66.66667%
0.00001 2.48832 7.59375 64.36343 243.00000 525.21875 2293.45007
77.77778% 88.88889% 100%
4546.12854 8445.96301 30770.56399

これを応用すれば、データの順位に応じた色割り当てが可能。

# 20分位作成
n = 20

# データは重複した値もあるので、分位の値が同じ値を持つ可能性がある。
# 例えば、100個のデータがあって、0が20個あれば下から10番目の値も20番目の値も0。
# uniqueしておく。
breaks = unique(quantile(as.matrix(demo2), probs = seq(0,1,length.out=n)))

# uniqueで数が減っている可能性があるので、作成したbreaksの長さを取り直しておく。
n2 = length(breaks)

pheatmap(demo2, show_rownames = F,
         breaks = breaks,
         color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))(n2)
         )

legend内の色の偏りはすごいが、ヒートマップ内の表現力は改善する。