京橋のバイオインフォマティシャンの日常

南国のビーチパラソルの下で、Rプログラムを打ってる日常を求めて、、Daily Life of Bioinformatician in Kyobashi of Osaka

R言語でヒートマップ入門:インタラクティブなヒートマップでデータ分析をもっと楽しくしよう!!

はじめに

過去に、heatmap.2パッケージと階層的クラスタリングを使って、ヒートマップの描く方を紹介しました。 こんなヒートマップでした。いま思うと、懐かしいような、古めかしい感じがしますね。

 

今回は、R言語でのヒートマップの特集ということで、 色々なパッケージのヒートマップ表現を調査して、取り上げてみました。

ヒートマップとは、データを可視化するため、 「数値データの高い・低いの表」を「色や濃淡」セルのグリットに変換して表現した視覚化法、あるいはデータ解析法の1つです。 数値の一覧を目視するよりも、データ間の意味や差異がより直感的に導けますし、 「色や濃淡」の曖昧な表現から、データ全体のニュアンスを掴む行為とも言えます。

また、項目間の近縁関係のパターンを強調するために並べられ、しばしばデンドログラムが添えられます。 こうすることで、ヒートマップの見え方が大きく変わり、データからふとした新しい発見にも繋がるかもで、 それは、このヒートマップ解析の興味深いところです。 ヒートマップは、観測値、相関関係、欠損値パターンなどを可視化するために、多くの分野で使用されています。

この記事では、色々なパッケージでのヒートマップの作り方に加えて、 インタラクティブにヒートマップを作成・可視化できるパッケージを紹介します。 全体を通じて、ツール間の比較ができ、良し悪しも感じてもらえるのではと思います。

実行環境

実行環境は以下の通りです。

macOS Big Sur (バージョン11.5.2)
MacBook Air (M1, 2020)
チップ Apple M1
メモリ 16GB

Rパッケージの読み込み

R環境を立ち上げて、必要なパッケージをインストールしておきます。

#パッケージ・インストール
pack <- c("htmlwidgets", "qtlcharts", "qtl", "iheatmapr", "heatmap3", "heatmaply", "remotes", "corrplot")
install.packages(pack[!(pack %in% unique(rownames(installed.packages())))])

#ロード
for(n in 1:length(pack)){ eval(parse(text = paste0("library(", pack[n], ")"))) }; rm("n", "pack")

heatmap3パッケージを用いた、ハイカラなヒートマップ表現

heatmap2というパッケージがありますが、最近、heatmap3というのも出ているようです。

heatmap3は、改良されたヒートマップパッケージで、Rのオリジナル関数'heatmap'と 完全な互換性で、便利な機能を提供すると説明されています。

それでは、早速、heatmap3のサンプルコードを実行してみます。

heatmap3でのヒートマップ作成(1)シンプルなヒートマップ

この実行例では、乱数を生成して、サンプルデータを作成・利用して、シンプルなヒートマップを作成します。

#乱数でサンプルデータ作成
rnormData <- matrix(rnorm(1000), 40, 25)
#head(rnormData)

#微調整
rnormData[1:15, seq(6, 25, 2)] <- rnormData[1:15, seq(6, 25, 2)] + 2
rnormData[16:40, seq(7, 25, 2)] = rnormData[16:40, seq(7, 25, 2)] + 4

#列名を与える
colnames(rnormData) <- c(paste("Control", 1:5, sep = ""),
                         paste(c("TrtA", "TrtB"), rep(1:10,each=2), sep = ""))
#行名を与える
rownames(rnormData) <- paste("Probe", 1:40, sep = "")

#グループのカラーリング
ColSideColors <- cbind(Group1=c(rep("steelblue2",5), rep(c("brown1", "mediumpurple2"),10)),
                       Group2=sample(c("steelblue2","brown1", "mediumpurple2"),25,replace=TRUE))

#特定セルのカラーリング
colorCell <- data.frame(row=c(1,3,5),col=c(2,4,6),
                      color=c("green4","black","orange2"),stringsAsFactors=FALSE)

#特定セルのハイライト
highlightCell <- data.frame(row=c(2,4,6),col=c(1,3,5),
                            color=c("black", "green4","orange2"),lwd=1:3,stringsAsFactors=FALSE)

#ヒートマップ実行
heatmap3(rnormData,
         ColSideColors=ColSideColors, 
         showRowDendro=FALSE,
         colorCell=colorCell,
         highlightCell=highlightCell)

#注釈
#ColSideColors: 水平サイドバーの色名を含む長さncol(x)の文字ベクトル
#showRowDendro: 行のデンドログラムをプロットするかどうかを指定する論理値
#colorCell: 3つの列を持つ data.frame で、どのセルが特定の色で着色されるか
#highlightCell: 3列または4列の data.frame で、どのセルが特定の色の矩形でハイライトされるか

 

このヒートマップはインタラクティブではない分、見栄えが良いような気がしますね。

heatmap3でのヒートマップ作成(2)少し複雑なヒートマップ

次には、少し複雑なヒートマップの作成を実行します。

#乱数でサンプルデータ作成
ColSideAnn <- data.frame(Information=rnorm(25),
                         Group=c(rep("Control",5),rep(c("TrtA","TrtB"),10)),
                         stringsAsFactors=TRUE)

#行名を与える(上記の実行から使用した変数転用)
row.names(ColSideAnn) <- colnames(rnormData)

#カラーパレット作成
RowSideColors <- colorRampPalette(c("chartreuse4", "white", "firebrick"))(40)

result <- heatmap3(rnormData,
                   ColSideCut=1.2,
                   ColSideAnn=ColSideAnn,
                   ColSideFun=function(x) {showAnn(x)},
                   ColSideWidth=0.8,
                   RowSideColors=RowSideColors,
                   col=colorRampPalette(c("green","black", "red"))(1024),
                   RowAxisColors=1,
                   legendfun=function() {showLegend(legend=c("Low","High"), col=c("chartreuse4","firebrick"))},
                   verbose=TRUE)
#注釈
#ColSideCut: coloumデンドログラムをカットする際に使用される数値
#ColSideAnn: 連続変数と因子変数をアノテーション情報として持つデータフレーム。
#ColSideFun: 列方向に注釈やラベリング図を生成するために使用される関数
#ColSideWidth: 数値列の側面領域の高さ
#RowSideColors: (オプション) nrow(x) の長さの文字ベクトル
#RowAxisColors: RowSideColors のどの色を行軸のラベルの色として使用するかを示す整数値
#legendfun: 関数で、図の左上に凡例を生成する。

 

実際、ヒートマップの調整は結構ややこしいです。

heatmaplyパッケージを使った、ヒートマップ作図

次に、heatmaplyパッケージは、ggplot2とplotly.jsのエンジンをベースにして開発されています。 d3heatmapと同様のヒートマップを生成しますが、plotly.jsはより大きなサイズの行列を扱えて、 デンドログラムからのズーム機能の点で優れているようです。

デフォルトのplot_method = "ggplot"の代わりに、 plot_method = "plotly"を選択することもできます。

heatmaplyによる通常のヒートマップ

まずは、mtcarsのデータセットで、heatmaply関数のシンプル実行をしてみます。

#データ表示
head(mtcars)
#                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
#Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
#Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
#Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
#Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
#Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
#Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

#実行
heatmaply(mtcars)

https://kumes.github.io/Blog/R_heatmap/04_heatmaply.html

heatmaply_corによる相関ヒートマップ

heatmaply_corは、相関行列に適したヒートマップ関数です。 具体的には、限界値は -1 から 1 に設定され、 カラーパレットはRdBuになっています。

#実行
heatmaply_cor(
  x = cor(mtcars),
  xlab = "Features",
  ylab = "Features",
  k_col = 2,
  k_row = 2
)

https://kumes.github.io/Blog/R_heatmap/05_heatmaply_cor.html

heatmaplyによるpercentizeヒートマップ

percentize関数を使った、経験的パーセンタイル変換は、変数のランクを取ることに似た挙動です。 違いは、変換された変数を比較したり解釈したりするのがよりシンプルになることです。

これは、ヒートマップ(例:heatmaply)において、複数の変数を比較するのに便利な側面があります。

#実行
heatmaply(
  percentize(mtcars),
  xlab = "Features",
  ylab = "Cars", 
  main = "Data transformation using 'percentize'"
)

https://kumes.github.io/Blog/R_heatmap/06_heatmaply_percentile.html

qtlchartsパッケージのiplotCorrを利用した、相関行列のヒートマップ表現

qtlchartsパッケージは、D3.jsを使った、ウェブベースのインタラクティブチャートを提供しています。 パッケージの主な用途は、遺伝子座を実験的に同定するための解析ツールのようです。 その中の1つの関数である、iplotCorrを使用します。

#データ読み込み
data(geneExpr)

#設定
chartOpts <- list(cortitle="Correlation matrix", scattitle="Scatterplot")

#ヒートマップ実行
iplotCorr(mat=geneExpr$expr, group=geneExpr$genotype, 
          reorder=TRUE, chartOpts=chartOpts)

https://kumes.github.io/Blog/R_heatmap/07_iplotCorr.html

次に、Spearmanの相関行列を使う方法です。

# Spearmanの相関行列計算
corr <- cor(geneExpr$expr, method="spearman", use="pairwise.complete.obs")

#設定
chartOpts <- list(cortitle="Spearman matrix", scattitle="Scatterplot")

# 階層的クラスタリングによる並び替え
o <- hclust(as.dist(1-corr))$order

#ヒートマップ実行
iplotCorr(mat=geneExpr$expr[,o], group=geneExpr$genotype, 
          corr=corr[o,o], chartOpts=chartOpts)

https://kumes.github.io/Blog/R_heatmap/08_iplotCorr.html

iheatmaprパッケージを用いたヒートマップ作成

iheatmaprパッケージは、複雑かつインタラクティブなヒートマップを構築するためのRパッケージです。

#データ読み込み
data(measles, package = "iheatmapr")

#ヘッド表示
head(measles)
#                1930  1931  1932  1933  1934
#ALABAMA         4389  8934   270  1735 15849
#ALASKA             0     0     0     0     0
#AMERICAN.SAMOA     0     0     0     0     0
#ARIZONA         2107  2135    86  1261  1022
#ARKANSAS         996   849    99  5438  7222
#CALIFORNIA     43585 27807 12618 26551 25650

#実行:メインのみ
main_heatmap(measles, name = "Measles<br>Cases", x_categorical = FALSE,
             layout = list(font = list(size = 8)))

まずは、main_heatmapのみの結果です。

https://kumes.github.io/Blog/R_heatmap/09_iheatmapr_main.html

#実行
main_heatmap(measles, name = "Measles<br>Cases", x_categorical = FALSE,
             layout = list(font = list(size = 8))) %>%
  add_col_groups(ifelse(1930:2001 < 1961,"No","Yes"),
                  side = "bottom", name = "Vaccine<br>Introduced?",
                  title = "Vaccine?",
                  colors = c("lightgray","blue")) %>%
  add_col_labels(ticktext = seq(1930,2000,10),font = list(size = 8)) %>%
  add_row_labels(size = 0.3,font = list(size = 6)) %>% 
  add_col_summary(layout = list(title = "Average<br>across<br>states"),
                  yname = "summary")  %>%                 
  add_col_title("Measles Cases from 1930 to 2001", side= "top") %>%
  add_row_summary(groups = TRUE, type = "bar",
                  layout = list(title = "Average<br>per<br>year",
                                font = list(size = 8)))

main_heatmapの結果に、6つのレイアーを付与して、最終的に複雑なヒートマップ図を作成します。

https://kumes.github.io/Blog/R_heatmap/10_iheatmapr.html

corrplotパッケージを用いた、相関ヒートマップ解析

corrplotは、相関行列を視覚的に探索するツールで、変数間の隠れたパターンを検出するのに役立ちます。また、自動的な変数の並べ替えをサポートしています。また、可視化方法、グラフィックレイアウト、色、凡例、テキストラベルなど、豊富なプロットオプションを提供しています。相関の統計的な有意性を判断するためのp値や信頼区間も提供しています。

#データ準備: mtcarsは自動車の燃費などに関するデータ
head(mtcars)
#                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
#Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
#Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
#Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
#Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
#Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
#Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

#相関の計算
M <- cor(mtcars)

#デフォルト相関行列
corrplot(M, method = 'circle')

#数字表示有りの相関行列
corrplot(M, method = 'number') # colorful number

corrmorantを使った、相関行列マップ

ヒートマップのパッケージとは言い難いのですが、、このパッケージも追記しておきます。

corrmorantパッケージは、通常のggplot2構文で相関行列のプロットを作図できます。 さらに、相関行列に基づく探索的データ解析のための大規模な可視化ツール群を提供しています。

#インストール
remotes::install_github("r-link/corrmorant")
library(corrmorant)

#データ読み込み
data(drosera)

#表示
head(drosera)
#   species variety petiole_length petiole_width blade_length blade_width
#1 capensis   rubra          54.00          1.95        38.30        3.95
#2 capensis   rubra          45.90          1.15        27.60        2.95
#3 capensis   rubra          21.20          0.85        21.85        2.65
#4 capensis   rubra          56.40          1.65        38.20        3.30
#5 capensis   rubra          28.60          1.00        15.80        2.80
#6 capensis   rubra          32.85          1.50        24.65        4.40

#詳細
str(drosera)
#'data.frame': 150 obs. of  6 variables:
# $ species       : Factor w/ 3 levels "capensis","madagascariensis",..: 1 1 1 1 1 1 1 1 1 1 ...
# $ variety       : Factor w/ 3 levels "alba","rubra",..: 2 2 2 2 2 2 2 2 2 2 ...
# $ petiole_length: num  54 45.9 21.2 56.4 28.6 ...
# $ petiole_width : num  1.95 1.15 0.85 1.65 1 1.5 1.85 1.65 2.4 1.5 ...
# $ blade_length  : num  38.3 27.6 21.9 38.2 15.8 ...
# $ blade_width   : num  3.95 2.95 2.65 3.3 2.8 4.4 4.8 3.4 3.9 4.4 ...

それでは、droseraのデータを使って、ggcorrmを実行します。 mappingで、colとfillにはspeciesの列を指定します。 2列目のvariety列が無視されるのは、Factorだからですかね。

#実行
ggcorrm(drosera, mapping = aes(col = species, fill = species)) +
  lotri(geom_smooth(alpha = 0.4, method = "lm")) +
  lotri(geom_point(alpha = 0.4)) +
  utri_corrtext(nrow = 2, squeeze = 0.6) +
  dia_names(y_pos = 0.15, size = 3) +
  dia_density(lower = 0.3, color = 1, alpha = 0.4)

#注釈
#lotri: プロットパネルのサブセットへのggplotレイヤーのマッピングを変更するために使用できるセレクタ関数
#utri_corrtext: プロットの非対角成分での二変量相関の強さをテキストラベルで表示するために使用されます。
#dia_names: プロットの変数名を対角成分の適切な位置にプロットします。
#dia_density: プロットの対角パネルに密度カーブを追加します。

以下のような、結果が得られます。

まとめ

R言語でのヒートマップ特集ということで、最近流行りのヒートマップ表現について調査してみました。 ヒートマップもインタラクティブな可視化が多くなってきている印象ですね。

1つのヒートマップ上に、色々な情報や解析結果を盛り込めるのは良いのですが、 実行コードが複雑で、覚えるのがやや面倒な感じがします。私だけでしょうか。。

参考資料

cran.r-project.org

kbroman.org

github.com

www.karada-good.net

github.com

cran.r-project.org

cran.r-project.org

github.com