【Plot】[pheatmap] アスタリスクをヒートマップ内に表示

pheatmapでは値を色で表現するだけでなく、文字列をplot内に埋め込むことができる。特定の値以上の箇所にアスタリスクなどを入れれば注目すべき箇所を強調できる。

ライブラリ

library(pheatmap)
library(psych)

【特定の値以上の箇所にアスタリスクをつける】

例としてmtcarsの各列間の全組み合わせで計算した相関係数を使ってみる。

相関係数を算出するcor()関数をdata.frameに対して行うと、各列間の相関係数が簡単に算出できる。

cor(mtcars)

各列間の相関係数

この相関係数matrixをpheatmap()で図示すると相関関係をわかりやすく表現できる。

cor_df <- cor(mtcars)
pheatmap(cor_df)

相関係数がある一定以上の値を持つ箇所にアスタリスクをつけてみる。

まずはpheatmapに渡したマトリクスと同じサイズの空のマトリクスを用意する。

空セルのmatrixを用意

marks <- matrix(data = "", 
                nrow = nrow(cor_df),
                ncol = ncol(cor_df))

そして元のマトリクスで特定の条件を満たす箇所と対応する空マトリクスに表示したい文字を埋めていく。
例えば次では相関係数が0.6より高い箇所の情報を使って、空マトリクスの対応箇所にアスタリスクを入れている。

marks[ cor_df > 0.6 ] <- "*"

この情報をpheatmap()機能のdisplay_numbers=引数に渡すと、目的のものが得られる。

pheatmap(cor_df,
         display_numbers = marks)

次のように多段階の閾値でアスタリスクを増やしていってもいいだろう。

marks[ abs(cor_df) > 0.3 ] <- "*"
marks[ abs(cor_df) > 0.6 ] <- "**"
marks[ abs(cor_df) > 0.9 ] <- "***"

pheatmap(cor_df,
         display_numbers = marks)

※ この時に条件が緩い➔厳しい順に代入していかなければならない。厳しい➔緩い順でやると緩い条件で入れたアスタリスクで更新されてしまう。

（応用例）【相関検定で有意水準を満たす箇所にアスタリスクをつける】

上記と同じmtcarsのデモデータを使って、相関検定のp値に基づいたアスタリスク付の例を紹介する。

data.frameの相関検定

Rにデフォルトで入っている相関検定のcor.test()関数は2つのベクトルにしか対応していない。

cor.test()にdata.frameを渡すとエラーが出る

data.frameの各列間の相関検定を行うにはpsychパッケージのcorr.test()関数が使用できる。
もしくはppcorパッケージのpcor()関数で偏相関係数の検定をしても良い。(ここでは解説しない。)

デフォルトではpearson相関検定をpairwiseに行い、多重検定補正はholm法で行う。

corr.testの引数確認

args(corr.test)

デフォルト設定で相関検定を行う。corr.test()の返り値はリストとなる。

corr.testの実行

hoge <- corr.test(mtcars)

rに相関係数、p.adjに調整済みp値が保存されている。

corr.test出力の中身

pheatmap用にp値のmatrixを成型

この調整済みp値を使ってヒートマップ上に有意水準の*を入れていきたいのだが、上三角行列の値のベクトルになっており、調整済みp値は (行数 * 列数 - 対角成分数)/2 の数しかない。

まずは調整済みp値を本来の行列サイズにしていく。

調整済みp値のmatrix作成

# 同じサイズのNAのmatrixを用意
tmp <- matrix(data = NA, nrow = nrow(hoge$p), ncol = ncol(hoge$p))
# 上三角行列に調整済みp値を代入
tmp[upper.tri(tmp)] <- hoge$p.adj

# 上三角行列を転置して、下三角行列に変換 (上三角行列の箇所はNA)
p <- t(tmp)
# 上三角行列のNAの箇所に調整済みp値を代入
p[upper.tri(p)] <- hoge$p.adj

完成した調整済みp値のmatrix

ちなみに、corr.test()の返り値pの行列は上三角が調整済みp値で下三角が調整前のp値なので、以下のようにしても同様の結果になる。

p値matrixの下三角に上三角の値を入れる

p <- t(hoge$p)
p[upper.tri(p)] <- hoge$p[upper.tri(hoge$p)]
# 対角成分の値が0になっているのでNAにしておく
diag(p) <- NA

次にp値の値に応じて*や空白セルに変換してゆく。

値を記号に変換

for(i in 1:nrow(p)){
    for(j in 1:ncol(p)){
        tmp <- as.numeric(p[i,j])
        if(is.na(tmp)){
            symbol <- ""
        } else if(tmp < 0.001){
            symbol <- "***"
        } else if(tmp < 0.01){
            symbol <- "**"
        } else if(tmp < 0.05){
            symbol <- "*"
        } else{
            symbol <- ""
        }
        p[i,j] <- symbol
    }
}

最後にpheatmap()のdisplay_numbers=引数に上記で作成したmatrixを指定して完成。

display_numbers引数

pheatmap(hoge$r, display_numbers = p)