r - 固定された外部階層クラスターを使用してヒートマップを作成する方法

Question

マトリックスデータがあり、ヒートマップで視覚化したいと思います。行は種であるため、行の横にある系統樹を視覚化し、ツリーに従ってヒートマップの行を並べ替えます。Rの関数で階層的クラスタリングヒートマップを作成できることは知っていheatmapますが、プロットでデフォルトで作成された距離クラスタリングの代わりに、系統的クラスタリングを使用するにはどうすればよいですか？

Answer 1

まず、パッケージを使用してデータをオブジェクトapeとして読み込む必要があります。phylo

library(ape)
dat <- read.tree(file="your/newick/file")
#or
dat <- read.tree(text="((A:4.2,B:4.2):3.1,C:7.3);")

以下は、ツリーが超距離である場合にのみ機能します。

次のステップは、系統樹をクラスに変換することdendrogramです。

次に例を示します。

data(bird.orders) #This is already a phylo object
hc <- as.hclust(bird.orders) #Compulsory step as as.dendrogram doesn't have a method for phylo objects.
dend <- as.dendrogram(hc)
plot(dend, horiz=TRUE)

mat <- matrix(rnorm(23*23),nrow=23, dimnames=list(sample(bird.orders$tip, 23), sample(bird.orders$tip, 23))) #Some random data to plot

まず、系統樹の順序に従ってマトリックスを順序付ける必要があります。

ord.mat <- mat[bird.orders$tip,bird.orders$tip]

次に、それをに入力しますheatmap：

heatmap(ord.mat, Rowv=dend, Colv=dend)

編集：これは、超距離および非超距離ツリーを処理するための関数です。

heatmap.phylo <- function(x, Rowp, Colp, ...){
    # x numeric matrix
    # Rowp: phylogenetic tree (class phylo) to be used in rows
    # Colp: phylogenetic tree (class phylo) to be used in columns
    # ... additional arguments to be passed to image function
    x <- x[Rowp$tip, Colp$tip]
    xl <- c(0.5, ncol(x)+0.5)
    yl <- c(0.5, nrow(x)+0.5)
    layout(matrix(c(0,1,0,2,3,4,0,5,0),nrow=3, byrow=TRUE),
                  width=c(1,3,1), height=c(1,3,1))
    par(mar=rep(0,4))
    plot(Colp, direction="downwards", show.tip.label=FALSE,
               xlab="",ylab="", xaxs="i", x.lim=xl)
    par(mar=rep(0,4))
    plot(Rowp, direction="rightwards", show.tip.label=FALSE, 
               xlab="",ylab="", yaxs="i", y.lim=yl)
    par(mar=rep(0,4), xpd=TRUE)
    image((1:nrow(x))-0.5, (1:ncol(x))-0.5, x, 
           xaxs="i", yaxs="i", axes=FALSE, xlab="",ylab="", ...)
    par(mar=rep(0,4))
    plot(NA, axes=FALSE, ylab="", xlab="", yaxs="i", xlim=c(0,2), ylim=yl)
    text(rep(0,nrow(x)),1:nrow(x),Rowp$tip, pos=4)
    par(mar=rep(0,4))
    plot(NA, axes=FALSE, ylab="", xlab="", xaxs="i", ylim=c(0,2), xlim=xl)
    text(1:ncol(x),rep(2,ncol(x)),Colp$tip, srt=90, pos=2)
    }

これが前の（超距離）例です：

heatmap.phylo(mat, bird.orders, bird.orders)

そして非超距離で：

cat("owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);",
    file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
mat2 <- matrix(rnorm(4*4),nrow=4, 
             dimnames=list(sample(tree.owls$tip,4),sample(tree.owls$tip,4)))
is.ultrametric(tree.owls)
[1] FALSE
heatmap.phylo(mat2,tree.owls,tree.owls)

Answer 2

まず、再現可能な例を作成します。データがなければ、私たちはあなたが望むものを推測することができます。ですので、次回はもっと上手くやってみてください（特にあなたは確認済みのユーザーです）。たとえば、これを実行して、newick形式でツリーを作成できます。

tree.text='(((XXX:4.2,ZZZ:4.2):3.1,HHH:7.3):6.3,AAA:13.6);'

@plannpusのように、私はapeこのツリーをhclustクラスに変換するために使用しています。残念ながら、変換は超距離ツリーに対してのみ実行できるようです。ルートから各先端までの距離は同じです。

library(ape)
tree <- read.tree(text='(((XXX:4.2,ZZZ:4.2):3.1,HHH:7.3):6.3,AAA:13.6);')
is.ultrametric(tree)
hc <- as.hclust.phylo(tree)

次に、dendrogramGrobfromを使用しlatticeExtraてツリーをプロットしています。ヒートマップを描画するためにlevelplot。lattice

library(latticeExtra)
dd.col <- as.dendrogram(hc)
col.ord <- order.dendrogram(dd.col)
mat <- matrix(rnorm(4*4),nrow=4)
colnames(mat) <- tree$tip.label
rownames(mat) <- tree$tip.label
levelplot(mat[tree$tip,tree$tip],type=c('g','p'),
          aspect = "fill",
          colorkey = list(space = "left"),
          legend =
            list(right =
                   list(fun = dendrogramGrob,
                        args =
                          list(x = dd.col, 
                               side = "right",
                               size = 10))),
          panel=function(...){
            panel.fill('black',alpha=0.2)
            panel.levelplot.points(...,cex=12,pch=23)
          }
)

Answer 3

私はplannapusの答えを複数のツリーを処理するように適応させました（また、プロセスで必要のないいくつかのオプションを切り取りました）：

library(ape)

heatmap.phylo <- function(x, Rowp, Colp, breaks, col, denscol="cyan", respect=F, ...){
    # x numeric matrix
    # Rowp: phylogenetic tree (class phylo) to be used in rows
    # Colp: phylogenetic tree (class phylo) to be used in columns
    # ... additional arguments to be passed to image function

    scale01 <- function(x, low = min(x), high = max(x)) {
        x <- (x - low)/(high - low)
        x
    }

    col.tip <- Colp$tip
    n.col <- 1
    if (is.null(col.tip)) {
        n.col <- length(Colp)
        col.tip <- unlist(lapply(Colp, function(t) t$tip))
        col.lengths <- unlist(lapply(Colp, function(t) length(t$tip)))
        col.fraction <- col.lengths / sum(col.lengths)
        col.heights <- unlist(lapply(Colp, function(t) max(node.depth.edgelength(t))))
        col.max_height <- max(col.heights)
    }

    row.tip <- Rowp$tip
    n.row <- 1
    if (is.null(row.tip)) {
        n.row <- length(Rowp)
        row.tip <- unlist(lapply(Rowp, function(t) t$tip))
        row.lengths <- unlist(lapply(Rowp, function(t) length(t$tip)))
        row.fraction <- row.lengths / sum(row.lengths)
        row.heights <- unlist(lapply(Rowp, function(t) max(node.depth.edgelength(t))))
        row.max_height <- max(row.heights)
    }

    cexRow <- min(1, 0.2 + 1/log10(n.row))
    cexCol <- min(1, 0.2 + 1/log10(n.col))

    x <- x[row.tip, col.tip]
    xl <- c(0.5, ncol(x)+0.5)
    yl <- c(0.5, nrow(x)+0.5)

    screen_matrix <- matrix( c(
        0,1,4,5,
        1,4,4,5,
        0,1,1,4,
        1,4,1,4,
        1,4,0,1,
        4,5,1,4
    ) / 5, byrow=T, ncol=4 )

    if (respect) {
        r <- grconvertX(1, from = "inches", to = "ndc") / grconvertY(1, from = "inches", to = "ndc")
        if (r < 1) {
            screen_matrix <- screen_matrix * matrix( c(r,r,1,1), nrow=6, ncol=4, byrow=T)
        } else {
            screen_matrix <- screen_matrix * matrix( c(1,1,1/r,1/r), nrow=6, ncol=4, byrow=T)
        }
    }


    split.screen( screen_matrix )

    screen(2)
    par(mar=rep(0,4))

    if (n.col == 1) {
        plot(Colp, direction="downwards", show.tip.label=FALSE,xaxs="i", x.lim=xl)
    } else {
        screens <- split.screen( as.matrix(data.frame( left=cumsum(col.fraction)-col.fraction, right=cumsum(col.fraction), bottom=0, top=1)))
        for (i in 1:n.col) {
            screen(screens[i])
            plot(Colp[[i]], direction="downwards", show.tip.label=FALSE,xaxs="i", x.lim=c(0.5,0.5+col.lengths[i]), y.lim=-col.max_height+col.heights[i]+c(0,col.max_height))
        }
    }

    screen(3)
    par(mar=rep(0,4))

    if (n.col == 1) {
        plot(Rowp, direction="rightwards", show.tip.label=FALSE,yaxs="i", y.lim=yl)
    } else {
        screens <- split.screen( as.matrix(data.frame( left=0, right=1, bottom=cumsum(row.fraction)-row.fraction, top=cumsum(row.fraction))) )
        for (i in 1:n.col) {
            screen(screens[i])
            plot(Rowp[[i]], direction="rightwards", show.tip.label=FALSE,yaxs="i", x.lim=c(0,row.max_height), y.lim=c(0.5,0.5+row.lengths[i]))
        }
    }


    screen(4)
    par(mar=rep(0,4), xpd=TRUE)
    image((1:nrow(x))-0.5, (1:ncol(x))-0.5, x, xaxs="i", yaxs="i", axes=FALSE, xlab="",ylab="", breaks=breaks, col=col, ...)

    screen(6)
    par(mar=rep(0,4))
    plot(NA, axes=FALSE, ylab="", xlab="", yaxs="i", xlim=c(0,2), ylim=yl)
    text(rep(0,nrow(x)),1:nrow(x),row.tip, pos=4, cex=cexCol)

    screen(5)
    par(mar=rep(0,4))
    plot(NA, axes=FALSE, ylab="", xlab="", xaxs="i", ylim=c(0,2), xlim=xl)
    text(1:ncol(x),rep(2,ncol(x)),col.tip, srt=90, adj=c(1,0.5), cex=cexRow)

    screen(1)
    par(mar = c(2, 2, 1, 1), cex = 0.75)

    symkey <- T
    tmpbreaks <- breaks
    if (symkey) {
        max.raw <- max(abs(c(x, breaks)), na.rm = TRUE)
        min.raw <- -max.raw
        tmpbreaks[1] <- -max(abs(x), na.rm = TRUE)
        tmpbreaks[length(tmpbreaks)] <- max(abs(x), na.rm = TRUE)
    } else {
        min.raw <- min(x, na.rm = TRUE)
        max.raw <- max(x, na.rm = TRUE)
    }
    z <- seq(min.raw, max.raw, length = length(col))

    image(z = matrix(z, ncol = 1), col = col, breaks = tmpbreaks, 
          xaxt = "n", yaxt = "n")
    par(usr = c(0, 1, 0, 1))
    lv <- pretty(breaks)
    xv <- scale01(as.numeric(lv), min.raw, max.raw)
    axis(1, at = xv, labels = lv)

    h <- hist(x, plot = FALSE, breaks = breaks)
    hx <- scale01(breaks, min.raw, max.raw)
    hy <- c(h$counts, h$counts[length(h$counts)])
    lines(hx, hy/max(hy) * 0.95, lwd = 1, type = "s", 
          col = denscol)
    axis(2, at = pretty(hy)/max(hy) * 0.95, pretty(hy))
    par(cex = 0.5)
    mtext(side = 2, "Count", line = 2)

    close.screen(all.screens = T)

}

tree <- read.tree(text = "(A:1,B:1);((C:1,D:2):2,E:1);((F:1,G:1,H:2):5,((I:1,J:2):2,K:1):1);", comment.char="")
N <- sum(unlist(lapply(tree, function(t) length(t$tip))))

set.seed(42)
m <- cor(matrix(rnorm(N*N), nrow=N))
rownames(m) <- colnames(m) <- LETTERS[1:N]
heatmap.phylo(m, tree, tree, col=bluered(10), breaks=seq(-1,1,length.out=11), respect=T) 

Answer 4

ヒートマップのこの正確なアプリケーションは、 GitHubでオープン/フリーに開発されたphyloseqパッケージのplot_heatmap関数（ggplot2に基づく）にすでに実装されています。完全なコードと結果の例はここに含まれています：

http://joey711.github.io/phyloseq/plot_heatmap-examples

1つの注意点であり、ここで明示的に求めているものではありませんがphyloseq::plot_heatmap、どちらの軸の階層ツリーもオーバーレイしません。階層的クラスタリングに基づいて軸の順序を決定しないのには十分な理由があります。これは、ノードでのブランチの回転方法に応じて、長いブランチの最後のインデックスを任意に並べることができるためです。この点、および非計量的多次元尺度構成法に基づく代替案は、NeatMapパッケージに関する記事でさらに説明されています。これもR用に記述されており、ggplot2を使用しています。ヒートマップ内のインデックスを順序付けるためのこの次元削減（調整）アプローチは、の系統発生的存在量データに適合していphyloseq::plot_heatmapます。

Answer 5

私の提案はphlyoseq::plot_heatmapそこへの道の一部になりますが、ツリーでデータを表現することが本当に目的である場合は、強力な「ggtree」パッケージでこれを実行できます。

いくつかの例は、次のggtreeドキュメントページの上部に示されています。

http://www.bioconductor.org/packages/3.7/bioc/vignettes/ggtree/inst/doc/advanceTreeAnnotation.html

私はggtreedevとはまったく関係がないことに注意してください。プロジェクトのファンであり、すでに何ができるのか。

Answer 6

@plannapusと通信した後xlab=""、上記のコードに関するいくつかの余分な情報を削除するために、コードを（ほんの少しだけ）変更しました。ここにコードがあります。コメントされた行に余分なコードがあり、新しい行がそれらを消去しているのがわかります。これが私のような新しいユーザーに役立つことを願っています！:)

heatmap.phylo <- function(x, Rowp, Colp, ...){
    # x numeric matrix
    # Rowp: phylogenetic tree (class phylo) to be used in rows
    # Colp: phylogenetic tree (class phylo) to be used in columns
    # ... additional arguments to be passed to image function
    x <- x[Rowp$tip, Colp$tip]
    xl <- c(0.5, ncol(x) + 0.5)
    yl <- c(0.5, nrow(x) + 0.5)
    layout(matrix(c(0,1,0,2,3,4,0,5,0),nrow = 3, byrow = TRUE),
                  width = c(1,3,1), height = c(1,3,1))
    par(mar = rep(0,4))
    # plot(Colp, direction = "downwards", show.tip.label = FALSE,
    #            xlab = "", ylab = "", xaxs = "i", x.lim = xl)
      plot(Colp, direction = "downwards", show.tip.label = FALSE,
               xaxs = "i", x.lim = xl)
    par(mar = rep(0,4))
    # plot(Rowp, direction = "rightwards", show.tip.label = FALSE, 
    #            xlab = "", ylab = "", yaxs = "i", y.lim = yl)
    plot(Rowp, direction = "rightwards", show.tip.label = FALSE, 
               yaxs = "i", y.lim = yl)
    par(mar = rep(0,4), xpd = TRUE)
    image((1:nrow(x)) - 0.5, (1:ncol(x)) - 0.5, x, 
           #xaxs = "i", yaxs = "i", axes = FALSE, xlab = "", ylab = "", ...)
           xaxs = "i", yaxs = "i", axes = FALSE, ...)
    par(mar = rep(0,4))
    plot(NA, axes = FALSE, ylab = "", xlab = "", yaxs = "i", xlim = c(0,2), ylim = yl)
    text(rep(0, nrow(x)), 1:nrow(x), Rowp$tip, pos = 4)
    par(mar = rep(0,4))
    plot(NA, axes = FALSE, ylab = "", xlab = "", xaxs = "i", ylim = c(0,2), xlim = xl)
    text(1:ncol(x), rep(2, ncol(x)), Colp$tip, srt = 90, pos = 2)
}