У меня есть матричные данные, и я хочу визуализировать ее с помощью тепловой карты. Строки - это виды, поэтому я хочу визуализировать филогенетическое дерево в стороне от строк и переупорядочить строки тепловой карты в соответствии с деревом. Я знаю, что функция heatmap
в R может создать тепловую карту иерархической кластеризации, но как я могу использовать свою филогенетическую кластеризацию вместо созданной по умолчанию дистанционной кластеризации на графике?
как создать тепловую карту с фиксированным внешним иерархическим кластером
Ответы (6)
Сначала вам нужно использовать пакет 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)
heatmap.phylo
функция! новый подход, независимый от концепции депрограммы! Я почти уверен, что перенесу это в мир сетки! +10! Я почти уверен, что могу передать его в пакет сетки (решетка и сетка не уверены в ggplot2)
- person agstudy; 01.03.2013
Error in image.default((1:ncol(x)) - 0.5, (1:nrow(x)) - 0.5, x, xaxs = "i", : dimensions of z are not length(x)(-1) times length(y)(-1)
работает heatmap.phylo(c, d1, d2)
на моих собственных данных. Я проверил, что размер матрицы и длина вершин двух деревьев определенно совпадают. Вы знаете, в чем может быть проблема? Спасибо.
- person RNA; 02.03.2013
heatmap.phylo()
ncol и nrow были переключены в функции image
. Я поправил.
- person RNA; 02.03.2013
heatmap.phylo()
функция. Однако, когда я попытался последовать вашему рабочему примеру, я получил следующее сообщение об ошибке: Error in plot.default(0, type = "n", xlim = x.lim, ylim = y.lim, xlab = "", : formal argument "xlab" matched by multiple actual arguments
. Я понятия не имею, откуда взялась ошибка, потому что я просто копирую / вставляю сюда ваш пример для не ультраметрического дерева ... Любые идеи?
- person Antonio Canepa; 14.09.2020
ape
с тех пор изменил свой код для построения филогении. Просто избавьтесь от xlab='', ylab=''
в строках 13 и 16 кода heatmap.phylo
, и он должен работать.
- person plannapus; 14.09.2020
Сначала я создаю воспроизводимый пример. Без данных мы можем просто угадать, чего вы хотите. Поэтому, пожалуйста, постарайтесь сделать лучше в следующий раз (особенно если вы являетесь подтвержденным пользователем). Например, вы можете сделать это, чтобы создать свое дерево в формате 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)
Затем я использую dendrogramGrob из 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)
}
)
Я адаптировал ответ 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)
tree
, выполнив: 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="")
, у меня появилось это сообщение об ошибке Error in if (z[i]) { : missing value where TRUE/FALSE needed
. Любые идеи?
- person Antonio Canepa; 14.09.2020
Это точное применение тепловой карты уже реализовано в функции plot_heatmap
(на основе ggplot2) в пакет phyloseq, который открыто / свободно разрабатывается на GitHub. Примеры с полным кодом и результатами включены здесь:
http://joey711.github.io/phyloseq/plot_heatmap-examples
Одно предостережение, а не то, о чем вы здесь явно просите, но phyloseq::plot_heatmap
не перекрывает иерархическое дерево ни для одной из осей. Есть веская причина не основывать порядок осей на иерархической кластеризации - и это связано с тем, что индексы в конце длинных ветвей могут по-прежнему располагаться рядом друг с другом произвольно, в зависимости от того, как ветки расположены. повернуты в узлах. Этот момент и альтернатива, основанная на неметрическом многомерном масштабировании, объясняется далее в статье о Пакет NeatMap, который также написан для R и использует ggplot2. Такой подход к упорядочиванию индексов на тепловой карте с уменьшением размеров (ординацией) адаптирован для филогенетических данных о численности в phyloseq::plot_heatmap
.
plot_heatmap
создаст тепловую карту без иерархического дерева кластеризации рядом с ней, но не может (по запросу OP) кластеризовать по филогении (или поместить филогенетическое дерево рядом с графиком, чтобы указать эту филогению). Звучит правильно или я что-то упускаю?
- person ohnoplus; 23.12.2017
phyloseq::plot_heatmap
может упорядочивать таксоны на тепловой карте в соответствии с их порядком в дереве. Это достигается с помощью команды taxa.order
, которая может принимать либо таксономический ранг для кластеризации индексов, либо произвольный порядок самих индексов. rdocumentation.org/packages/phyloseq/versions/1.16.2/ темы / github.com/joey711/phyloseq/issues/230
- person Paul McMurdie; 08.01.2018
Хотя мое предложение для phlyoseq::plot_heatmap
поможет вам в этом, мощный пакет "ggtree" может сделать это или даже больше, если представление данных на деревьях - это действительно то, что вы собираетесь делать.
Некоторые примеры показаны в верхней части следующей страницы документации ggtree:
http://www.bioconductor.org/packages/3.7/bioc/vignettes/ggtree/inst/doc/advanceTreeAnnotation.html
Обратите внимание, что я вообще не связан с ggtree dev. Просто фанат проекта и того, что он уже умеет.
После общения с @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)
}
reorderfun
в тепловой карте ... - person Roman Luštrik   schedule 01.03.2013dput(head(mymatrixdata))
позволит людям легко восстановить часть ваших данных и упростит им помощь. - person Simon O'Hanlon   schedule 01.03.2013(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);
- person RNA   schedule 01.03.2013