加入树状图和热图

我有heatmap(来自一组样本的基因表达):

set.seed(10)
mat <- matrix(rnorm(24*10,mean=1,sd=2),nrow=24,ncol=10,dimnames=list(paste("g",1:24,sep=""),paste("sample",1:10,sep="")))
dend <- as.dendrogram(hclust(dist(mat)))
row.ord <- order.dendrogram(dend)
mat <- matrix(mat[row.ord,],nrow=24,ncol=10,dimnames=list(rownames(mat)[row.ord],colnames(mat)))
mat.df <- reshape2::melt(mat,value.name="expr",varnames=c("gene","sample"))

require(ggplot2)
map1.plot <- ggplot(mat.df,aes(x=sample,y=gene))+geom_tile(aes(fill=expr))+scale_fill_gradient2("expr",high="darkred",low="darkblue")+scale_y_discrete(position="right")+
  theme_bw()+theme(plot.margin=unit(c(1,1,1,-1),"cm"),legend.key=element_blank(),legend.position="right",axis.text.y=element_blank(),axis.ticks.y=element_blank(),panel.border=element_blank(),strip.background=element_blank(),axis.text.x=element_text(angle=45,hjust=1,vjust=1),legend.text=element_text(size=5),legend.title=element_text(size=8),legend.key.size=unit(0.4,"cm"))

(由于plot.margin我正在使用的论据,左侧被切断,但我需要这个,如下所示).

然后,我prune行dendrogram根据深度截止值来获得较少的集群(即,只有深深的分裂),并做一些编辑所产生dendrogram有它绘制他们的方式,我希望它:

depth.cutoff <- 11
dend <- cut(dend,h=depth.cutoff)$upper
require(dendextend)
gg.dend <- as.ggdend(dend)
leaf.heights <- dplyr::filter(gg.dend$nodes,!is.na(leaf))$height
leaf.seqments.idx <- which(gg.dend$segments$yend %in% leaf.heights)
gg.dend$segments$yend[leaf.seqments.idx] <- max(gg.dend$segments$yend[leaf.seqments.idx])
gg.dend$segments$col[leaf.seqments.idx] <- "black"
gg.dend$labels$label <- 1:nrow(gg.dend$labels)
gg.dend$labels$y <- max(gg.dend$segments$yend[leaf.seqments.idx])
gg.dend$labels$x <- gg.dend$segments$x[leaf.seqments.idx]
gg.dend$labels$col <- "black"
dend1.plot <- ggplot(gg.dend,labels=F)+scale_y_reverse()+coord_flip()+theme(plot.margin=unit(c(1,-3,1,1),"cm"))+annotate("text",size=5,hjust=0,x=gg.dend$label$x,y=gg.dend$label$y,label=gg.dend$label$label,colour=gg.dend$label$col)

我用cowplot's将它们绘制在一起plot_grid:

require(cowplot)
plot_grid(dend1.plot,map1.plot,align='h',rel_widths=c(0.5,1))

虽然align='h'它的工作并不完美.

使用dendrogramwith 绘制未切割图表说明了这一点:map1.plotplot_grid

dend0.plot <- ggplot(as.ggdend(dend))+scale_y_reverse()+coord_flip()+theme(plot.margin=unit(c(1,-1,1,1),"cm"))
plot_grid(dend0.plot,map1.plot,align='h',rel_widths=c(1,1))

顶部和底部的树枝dendrogram似乎被挤向中心.使用它scale似乎是一种调整它的方式,但是比例值似乎是特定于数字的,所以我想知道是否有任何方式以更有原则的方式做到这一点.

接下来,我对我的每个群集进行一些术语丰富分析heatmap.假设这个分析给了我这个data.frame:

enrichment.df <- data.frame(term=rep(paste("t",1:10,sep=""),nrow(gg.dend$labels)),
                          cluster=c(sapply(1:nrow(gg.dend$labels),function(i) rep(i,5))),
                          score=rgamma(10*nrow(gg.dend$labels),0.2,0.7),
                          stringsAsFactors = F)

我想做的是将其绘制data.frame为a heatmap并将切口dendrogram放在它下面(类似于它放在表达式左边的方式heatmap).

所以我plot_grid再次尝试认为这align='v'可以在这里工作:

首先重新生成树状图,使其面朝上:

dend2.plot <- ggplot(gg.dend,labels=F)+scale_y_reverse()+theme(plot.margin=unit(c(-3,1,1,1),"cm"))

现在尝试将它们一起绘制:

plot_grid(map2.plot,dend2.plot,align='v')

plot_grid 似乎无法对齐它们,如图所示和它抛出的警告消息:

In align_plots(plotlist = plots, align = align) :
  Graphs cannot be vertically aligned. Placing graphs unaligned.

似乎接近的是:

plot_grid(map2.plot,dend2.plot,rel_heights=c(1,0.5),nrow=2,ncol=1,scale=c(1,0.675))

这是在scale玩完参数后实现的,虽然情节太宽了.如此反复,我想知道是否有办法解决它,或在某种程度上预先确定什么是正确scale的的任何给定的名单dendrogram和heatmap,也许它们的尺寸.

前段时间我遇到了同样的问题.我使用的基本技巧是在给出树形图的结果的情况下直接指定基因的位置.为简单起见,首先是绘制完整树状图的情况:

# For the full dendrogram
library(plyr)
library(reshape2)
library(dplyr)
library(ggplot2)
library(ggdendro)
library(gridExtra)
library(dendextend)

set.seed(10)

# The source data
mat <- matrix(rnorm(24 * 10, mean = 1, sd = 2), 
              nrow = 24, ncol = 10, 
              dimnames = list(paste("g", 1:24, sep = ""), 
                              paste("sample", 1:10, sep = "")))

sample_names <- colnames(mat)

# Obtain the dendrogram
dend <- as.dendrogram(hclust(dist(mat)))
dend_data <- dendro_data(dend)

# Setup the data, so that the layout is inverted (this is more 
# "clear" than simply using coord_flip())
segment_data <- with(
    segment(dend_data), 
    data.frame(x = y, y = x, xend = yend, yend = xend))
# Use the dendrogram label data to position the gene labels
gene_pos_table <- with(
    dend_data$labels, 
    data.frame(y_center = x, gene = as.character(label), height = 1))

# Table to position the samples
sample_pos_table <- data.frame(sample = sample_names) %>%
    mutate(x_center = (1:n()), 
           width = 1)

# Neglecting the gap parameters
heatmap_data <- mat %>% 
    reshape2::melt(value.name = "expr", varnames = c("gene", "sample")) %>%
    left_join(gene_pos_table) %>%
    left_join(sample_pos_table)

# Limits for the vertical axes
gene_axis_limits <- with(
    gene_pos_table, 
    c(min(y_center - 0.5 * height), max(y_center + 0.5 * height))
) + 
    0.1 * c(-1, 1) # extra spacing: 0.1

# Heatmap plot
plt_hmap <- ggplot(heatmap_data, 
                   aes(x = x_center, y = y_center, fill = expr, 
                       height = height, width = width)) + 
    geom_tile() +
    scale_fill_gradient2("expr", high = "darkred", low = "darkblue") +
    scale_x_continuous(breaks = sample_pos_table$x_center, 
                       labels = sample_pos_table$sample, 
                       expand = c(0, 0)) + 
    # For the y axis, alternatively set the labels as: gene_position_table$gene
    scale_y_continuous(breaks = gene_pos_table[, "y_center"], 
                       labels = rep("", nrow(gene_pos_table)),
                       limits = gene_axis_limits, 
                       expand = c(0, 0)) + 
    labs(x = "Sample", y = "") +
    theme_bw() +
    theme(axis.text.x = element_text(size = rel(1), hjust = 1, angle = 45), 
          # margin: top, right, bottom, and left
          plot.margin = unit(c(1, 0.2, 0.2, -0.7), "cm"), 
          panel.grid.minor = element_blank())

# Dendrogram plot
plt_dendr <- ggplot(segment_data) + 
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) + 
    scale_x_reverse(expand = c(0, 0.5)) + 
    scale_y_continuous(breaks = gene_pos_table$y_center, 
                       labels = gene_pos_table$gene, 
                       limits = gene_axis_limits, 
                       expand = c(0, 0)) + 
    labs(x = "Distance", y = "", colour = "", size = "") +
    theme_bw() + 
    theme(panel.grid.minor = element_blank())

library(cowplot)
plot_grid(plt_dendr, plt_hmap, align = 'h', rel_widths = c(1, 1))

请注意,我在热图图表的左侧保持y轴刻度,只是为了显示树状图和刻度线完全匹配.

现在,对于剪切树状图的情况,应该记住,树状图的叶子将不再以对应于给定聚类中的基因的确切位置结束.为了获得基因和簇的位置,需要从切割完整的树形图中的两个树形图中提取数据.总的来说,为了澄清簇中的基因,我添加了划分簇的矩形.

# For the cut dendrogram
library(plyr)
library(reshape2)
library(dplyr)
library(ggplot2)
library(ggdendro)
library(gridExtra)
library(dendextend)

set.seed(10)

# The source data
mat <- matrix(rnorm(24 * 10, mean = 1, sd = 2), 
              nrow = 24, ncol = 10, 
              dimnames = list(paste("g", 1:24, sep = ""), 
                              paste("sample", 1:10, sep = "")))

sample_names <- colnames(mat)

# Obtain the dendrogram
full_dend <- as.dendrogram(hclust(dist(mat)))

# Cut the dendrogram
depth_cutoff <- 11
h_c_cut <- cut(full_dend, h = depth_cutoff)
dend_cut <- as.dendrogram(h_c_cut$upper)
dend_cut <- hang.dendrogram(dend_cut)
# Format to extend the branches (optional)
dend_cut <- hang.dendrogram(dend_cut, hang = -1) 
dend_data_cut <- dendro_data(dend_cut)

# Extract the names assigned to the clusters (e.g., "Branch 1", "Branch 2", ...)
cluster_names <- as.character(dend_data_cut$labels$label)
# Extract the names of the haplotypes that belong to each group (using
# the 'labels' function)
lst_genes_in_clusters <- h_c_cut$lower %>% 
    lapply(labels) %>% 
    setNames(cluster_names)

# Setup the data, so that the layout is inverted (this is more 
# "clear" than simply using coord_flip())
segment_data <- with(
    segment(dend_data_cut), 
    data.frame(x = y, y = x, xend = yend, yend = xend))

# Extract the positions of the clusters (by getting the positions of the 
# leafs); data is already in the same order as in the cluster name
cluster_positions <- segment_data[segment_data$xend == 0, "y"]
cluster_pos_table <- data.frame(y_position = cluster_positions, 
                                cluster = cluster_names)

# Specify the positions for the genes, accounting for the clusters
gene_pos_table <- lst_genes_in_clusters %>%
    ldply(function(ss) data.frame(gene = ss), .id = "cluster") %>%
    mutate(y_center = 1:nrow(.), 
           height = 1)
# > head(gene_pos_table, 3)
#    cluster gene y_center height
# 1 Branch 1  g11        1      1
# 2 Branch 1  g20        2      1
# 3 Branch 1  g12        3      1

# Table to position the samples
sample_pos_table <- data.frame(sample = sample_names) %>%
    mutate(x_center = 1:nrow(.), 
           width = 1)

# Coordinates for plotting rectangles delimiting the clusters: aggregate
# over the positions of the genes in each cluster
cluster_delim_table <- gene_pos_table %>%
    group_by(cluster) %>%
    summarize(y_min = min(y_center - 0.5 * height), 
              y_max = max(y_center + 0.5 * height)) %>%
    as.data.frame() %>%
    mutate(x_min = with(sample_pos_table, min(x_center - 0.5 * width)), 
           x_max = with(sample_pos_table, max(x_center + 0.5 * width)))

# Neglecting the gap parameters
heatmap_data <- mat %>% 
    reshape2::melt(value.name = "expr", varnames = c("gene", "sample")) %>%
    left_join(gene_pos_table) %>%
    left_join(sample_pos_table)

# Limits for the vertical axes (genes / clusters)
gene_axis_limits <- with(
    gene_pos_table, 
    c(min(y_center - 0.5 * height), max(y_center + 0.5 * height))
) + 
    0.1 * c(-1, 1) # extra spacing: 0.1

# Heatmap plot
plt_hmap <- ggplot(heatmap_data, 
                   aes(x = x_center, y = y_center, fill = expr, 
                       height = height, width = width)) + 
    geom_tile() +
    geom_rect(data = cluster_delim_table, 
              aes(xmin = x_min, xmax = x_max, ymin = y_min, ymax = y_max), 
              fill = NA, colour = "black", inherit.aes = FALSE) + 
    scale_fill_gradient2("expr", high = "darkred", low = "darkblue") +
    scale_x_continuous(breaks = sample_pos_table$x_center, 
                       labels = sample_pos_table$sample, 
                       expand = c(0.01, 0.01)) + 
    scale_y_continuous(breaks = gene_pos_table$y_center, 
                       labels = gene_pos_table$gene, 
                       limits = gene_axis_limits, 
                       expand = c(0, 0), 
                       position = "right") + 
    labs(x = "Sample", y = "") +
    theme_bw() +
    theme(axis.text.x = element_text(size = rel(1), hjust = 1, angle = 45), 
          # margin: top, right, bottom, and left
          plot.margin = unit(c(1, 0.2, 0.2, -0.1), "cm"), 
          panel.grid.minor = element_blank())

# Dendrogram plot
plt_dendr <- ggplot(segment_data) + 
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend)) + 
    scale_x_reverse(expand = c(0, 0.5)) + 
    scale_y_continuous(breaks = cluster_pos_table$y_position, 
                       labels = cluster_pos_table$cluster, 
                       limits = gene_axis_limits, 
                       expand = c(0, 0)) + 
    labs(x = "Distance", y = "", colour = "", size = "") +
    theme_bw() + 
    theme(panel.grid.minor = element_blank())

library(cowplot)
plot_grid(plt_dendr, plt_hmap, align = 'h', rel_widths = c(1, 1.8))