ggtreeで描図した系統樹にgheatmapでアノテーションを付与する
はじめに
R言語のggtreeで描図した系統樹に、メタデータを読み込んでアノテーションを付与(ヒートマップとして表示)する方法の備忘録です。
動作環境 (Hardware)
-MacBook Pro
-チップ Apple M1 (Rosetta2 インストール済み)
-macOS Sequoia 15.6.1
-メモリ 16GB
動作環境 (Software)
-R version 4.5.1
-ggtree 3.16.3
-ggplot2 3.5.2
-ggnewscale 0.5.2
-ape 5.8.1
-RStudio
手順
1.パッケージの読み込み
Rを起動してパッケージを読み込みます。
#インストールが済んでいない場合は先にパッケージをインストール
#install.packages("ggnewscale")
#install.packages("ape")
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("ggtree")
#パッケージの読み込み
library(ggplot2)
library(ggnewscale)
library(ape)
library(ggtree)
packageVersion関数を使用すると各パッケージのバージョンを確認できます。
#packageVersionの例
packageVersion("ggplot2")

2.treeオブジェクトの作成と描図
系統樹のファイルをread.tree関数を使って読み込み、treeオブジェクトに代入して描図します。
tree <- read.tree('/path_to/treefile')
#描図
p <- ggtree(tree) + geom_tiplab(size = 4) + geom_treescale(width = 0.0005)
plot(p)
geom_tiplab関数でラベル(株名)を表示、
geom_treescale関数でscale barを表示します(ラベルの大きさやscale barの値は適宜調整)。
plotの結果です。

scale barの表示位置をずらすには、geom_treescale関数で xとyに座標の値を指定します(それぞれの値は適宜)。
p2 <- ggtree(tree) + geom_tiplab(size = 4) +
geom_treescale(width = 0.0005, x = -0.0001, y = -0.8)
plot(p2)
xとyの両方にマイナスの値を指定したところ、scale barが先ほどより左下側に移動しました。

3.メタデータの読み込み
以下のようなメタデータが csvファイルで作成されている(metadata.csv)とします。
系統樹のラベル名と、メタデータのファイルの株名の列のデータは揃えておきます。

data_4の列のように、ある特性を持っている(+)か否かで空欄にしている2値データの列がある場合、negative(ある特性を持っていない)データには事前にダミーの値(例では「0」)を入力して埋めておきます。

data_1のような、数値データのみの列(例えば、サンプルの分離年など)はそのままread.csv関数で読み込むと数値型のデータと判別されるので、colClassesを使用して文字列型(chr)として読み込まれるようにします。
#メタデータの読み込み
#colClassesのベクトルに、文字列型に変更させたいデータの列名 = "character" を指定する
meta_data <- read.csv('/path_to/metadata.csv', colClasses = c(data_1 = "character"))
#確認
meta_data

data_1の列も「chr」型で読み込まれています。
4.アノテーションの付与
メタデータの各列からデータフレームを作成して、gheatmap関数で系統樹にそれらのデータをヒートマップとして次々重ねて表示していきます。
最初のメタデータからデータフレームを作成します。
data_1 <- data.frame("data_1" = meta_data[, c("data_1")])
rownames(data_1) <- meta_data$name
上の1行目、「data_1」が3回出てきてややこしいですが
一つ目の「data_1」が変数(作成するデータフレーム)の名前
二つ目の「data_1」が作成するデータフレームの列名につける名前
三つ目の「data_1」がメタデータを代入した meta_data オブジェクトから抽出する列の名前
です。
2行目で、作成したデータフレームの行名をメタデータの「name」列の値にすることで系統樹のラベルと、それぞれの株のメタデータとを紐付けています。
データフレームを確認します。
#確認
data_1

次に、gheatmap関数とscale_fill_manual関数で表示の設定をします。
h1 <- gheatmap(p = p2,
data = data_1,
offset = 0.001,
width = 0.07,
color = "white",
colnames = TRUE,
colnames_position = "top",
font.size = 2) +
scale_fill_manual(
name = "data_1",
values = c("2000" = "#8b4513", "2001" = "#0000FF", "2002" = "green", "2003" = "red", "2004" = "yellow"),
breaks = c("2000", "2001", "2002", "2003", "2004"),
guide = guide_legend(position = "right",
keywidth = 1,
keyheight = 1,
order = 1)
)
gheatmapでは以下を指定します。
p:系統樹のオブジェクト
data:表示させるデータフレーム
offset:系統樹とヒートマップの間の距離
width:ヒートマップの表示幅
color:ヒートマップの各セルの区切り(グリッド)の色
colnames:TRUEでヒートマップの列名を表示させる
colnames_position:列名の表示位置
font.size:列名のフォントサイズ
scale_fill_manualでは以下を指定します。
name:凡例のタイトル
values:データフレーム内の各データとその表示色(RGB値または色名で指定)のベクトル
上記のコードのように、「#rrggbb表記」と「色の名前」を混ぜて指定可能
breaks:凡例に表示させるデータのベクトル
guide_legend:
position:系統樹に対しての凡例の表示位置
keywidthおよびkeyheight:凡例の表示サイズ
order:凡例の並び順
#描図
plot(h1)

(参考1)offsetの値を0.001から0.005に変更してみます。値を大きくすることで系統樹とヒートマップとを離すことができます。
参考1

(参考2)keywidthおよびkeyheightの値をそれぞれ1から0.5に変更してみます。凡例の表示が小さくなりました。
参考2

二つ目のメタデータからデータフレームを作成します。
data_2 <- data.frame("data_2" = meta_data[, c("data_2")])
rownames(data_2) <- meta_data$name
new_scale_fillで先の系統樹のオブジェクト(h1)に新しいデータの追加を宣言したのち、gheatmapのoffsetの値をずらしていくことで、二つ目のメタデータのヒートマップを重ねます。
h2 <- h1 + new_scale_fill()
h3 <- gheatmap(p = h2,
data = data_2,
offset = 0.00135, #0.001 + 0.00035
width = 0.07,
color = "white" ,
colnames = TRUE,
colnames_position = "top",
font.size = 2) +
scale_fill_manual(
name = "data_2",
values = c("A" = "#FFC0CB", "B" = "#40E0D0", "C" = "#006400", "D" = "#FFE4B5"),
breaks = c("A", "B", "C", "D"),
guide = guide_legend(keywidth = 1, keyheight = 1, order = 2)
)
#描図
plot(h3)

(参考3)offsetの増分を0.00035 から 0.0005にしてみます。
参考3

今回の例では offsetの増分を0.00035にすると、ヒートマップを隣接して描くことができるのがわかりました。
三つ目のメタデータからデータフレームを作成します。
data_3 <- data.frame("data_3" = meta_data[, c("data_3")])
rownames(data_3) <- meta_data$name
先ほどと同様に、new_scale_fillでこれまで作成したオブジェクト(h3)に新しいデータの追加を宣言し、gheatmapのoffsetの値をさらに今回の増分(0.00035)ずらすことで三つ目のメタデータのヒートマップを重ねます。
h4 <- h3 + new_scale_fill()
h5 <- gheatmap(p = h4,
data = data_3,
offset = 0.0017, #0.00135 + 0.00035
width = 0.07,
color = "white" ,
colnames = TRUE,
colnames_position = "top",
font.size = 2) +
scale_fill_manual(
name = "data_3",
values = c("a" = "#FF1493", "b" = "#00BFFF", "c" = "#F0E68C"),
breaks = c("a", "b", "c"),
guide = guide_legend(keywidth = 1, keyheight = 1, order = 3)
)
#描図
plot(h5)

(参考4)data_2のguide_legendのorderを3、data_3のorderを2にして図を作成してみます。凡例の並び順が変わっています。
参考4

四つ目のメタデータ(data_4)のような、2値データを追加します。
scale_fill_manualのvaluesの指定で、「0」を入力していた negativeの(ある特性を持っていない)データの色指定を、図のバックグラウンドと同じ色(白)に設定します。
また、guides(fill="none")と設定することで、凡例には非表示にすることができます。
data_4 <- data.frame("data_4" = meta_data[, c("data_4")])
rownames(data_4) <- meta_data$name
h6 <- h5 + new_scale_fill()
h7 <- gheatmap(p = h6,
data = data_4,
offset = 0.00205, #0.0017 + 0.00035
width = 0.07,
color = "white",
colnames = TRUE,
colnames_position = "top",
font.size = 2) +
scale_fill_manual(
name = "data_4",
values = c("+" = "orange", "0" = "white")
) +
guides(fill = "none")
#描図
plot(h7)

(参考5)valuesの値の記載を「values = c("+" = "orange")」と、"+"のデータだけ指定して図を作成してみます。他の、色を指定していないデータが自動で灰色に塗られました。
参考5

5.画像の保存
ggsave関数を利用して画像ファイルとして書き出します。plotに、最終的に出力する「系統樹+ヒートマップ」のオブジェクトを指定します。解像度やサイズは一例です。
ggsave("plot.png", plot = h7, dpi = 400, width = 10, height = 8)
plot.png

書き出した画像が、RStudioやJupyter Notebookなどで作業中にplotした画像とずれがあることがあるので、必要に応じて数値を微調整して修正します。
おわりに
-
系統樹ごとに値が異なるので試行錯誤の必要そうなパラメータ:
・geom_treescaleの width
・scale barの表示位置をずらす場合の xとy
・gheatmapの offset(最初の値と増分の値) および width -
ヒートマップの列名の表示は 表示する列数が多い・文字数が長いなどの場合、オプションで操作しても文字が被ったり画像から見切れたりすることがあります。
「colnames = FALSE」で表示しないようにして図を作成して、別途画像編集ソフトなどでテキストを追加する方が視認性がよい図になるケースもあると思います。
参考
gheatmap: gheatmap in ggtree: an R package for visualization of tree and annotation data
Create your own discrete scale — scale_manual • ggplot2
Discussion