美文网首页基因组数据绘图
跟着Nature Methods学画图:R语言画热图(pheat

跟着Nature Methods学画图:R语言画热图(pheat

作者: 小明的数据分析笔记本 | 来源:发表于2021-02-20 19:52 被阅读0次

在简书 土豆学生信 分享的内容看到这篇论文 简书的链接是 https://www.jianshu.com/p/bbf9cb13b41a

论文是


image.png

论文对应的代码是公开的 https://github.com/ajwilk/2020_Wilk_COVID

image.png

在学习他这个代码的时候发现其中自定义了一个函数可以操作热图的文字标签,可以让热图上只显示我们感兴趣的文字标签。

我在运行这个代码的时候遇到了报错,没有把代码完全运行完,但是已经获得和NK.markers这个表达量文件,部分内容如下

image.png

我们用这个表达量文件先做一个简单的热图

读入数据
df<-read.csv("NM/NK_markers_1.csv",header=T,row.names = 1)
head(df)
最简单的热图
library(pheatmap)
pdf(file = "NM/hp-1.pdf",width = 4,height = 10)
pheatmap(df,fontsize = 3)
dev.off()
image.png

我们可以看到上图右侧所有的基因名都显示出来了,如果我们想只显示自己感兴趣的,那该如何实现呢?可以用开头提到的自定义函数

add.flag <- function(pheatmap,
                     kept.labels,
                     repel.degree) {
  
  # repel.degree = number within [0, 1], which controls how much 
  #                space to allocate for repelling labels.
  ## repel.degree = 0: spread out labels over existing range of kept labels
  ## repel.degree = 1: spread out labels over the full y-axis
  
  heatmap <- pheatmap$gtable
  
  new.label <- heatmap$grobs[[which(heatmap$layout$name == "row_names")]] 
  
  # keep only labels in kept.labels, replace the rest with ""
  new.label$label <- ifelse(new.label$label %in% kept.labels, 
                            new.label$label, "")
  
  # calculate evenly spaced out y-axis positions
  repelled.y <- function(d, d.select, k = repel.degree){
    # d = vector of distances for labels
    # d.select = vector of T/F for which labels are significant
    
    # recursive function to get current label positions
    # (note the unit is "npc" for all components of each distance)
    strip.npc <- function(dd){
      if(!"unit.arithmetic" %in% class(dd)) {
        return(as.numeric(dd))
      }
      
      d1 <- strip.npc(dd$arg1)
      d2 <- strip.npc(dd$arg2)
      fn <- dd$fname
      return(lazyeval::lazy_eval(paste(d1, fn, d2)))
    }
    
    full.range <- sapply(seq_along(d), function(i) strip.npc(d[i]))
    selected.range <- sapply(seq_along(d[d.select]), function(i) strip.npc(d[d.select][i]))
    
    return(unit(seq(from = max(selected.range) + k*(max(full.range) - max(selected.range)),
                    to = min(selected.range) - k*(min(selected.range) - min(full.range)), 
                    length.out = sum(d.select)), 
                "npc"))
  }
  new.y.positions <- repelled.y(new.label$y,
                                d.select = new.label$label != "")
  new.flag <- segmentsGrob(x0 = new.label$x,
                           x1 = new.label$x + unit(0.15, "npc"),
                           y0 = new.label$y[new.label$label != ""],
                           y1 = new.y.positions)
  
  # shift position for selected labels
  new.label$x <- new.label$x + unit(0.2, "npc")
  new.label$y[new.label$label != ""] <- new.y.positions
  
  # add flag to heatmap
  heatmap <- gtable::gtable_add_grob(x = heatmap,
                                     grobs = new.flag,
                                     t = 4, 
                                     l = 4
  )
  
  # replace label positions in heatmap
  heatmap$grobs[[which(heatmap$layout$name == "row_names")]] <- new.label
  
  # plot result
  grid.newpage()
  grid.draw(heatmap)
  
  # return a copy of the heatmap invisibly
  invisible(heatmap)
}

将以上函数放到文本文件里,通过source()加载这个函数

source("useful_R_function/add_flag.r")

选择感兴趣的基因名,我这里就随机选取几个了

gene_name<-sample(rownames(df),10)

画图

source("useful_R_function/add_flag.r")
library(grid)
gene_name<-sample(rownames(df),10)
p1<-pheatmap(df)
add.flag(p1,
         kept.labels = gene_name,
         repel.degree = 0.2)

结果就变成了如下


image.png
接下来是简单的美化

代码

source("useful_R_function/add_flag.r")
df<-read.csv("NM/NK_markers_1.csv",header=T,row.names = 1)
head(df)
library(pheatmap)
library(grid)
gene_name<-sample(rownames(df),10)
paletteLength <- 100
mycolor<-colorRampPalette(c("blue","white","red"))(100)
mycolor
myBreaks <- unique(c(seq(min(df), 0, length.out=ceiling(paletteLength/2) + 1), 
                     seq(max(df)/paletteLength, max(df),
                         length.out=floor(paletteLength/2))))
p1<-pheatmap(df,color = mycolor,breaks = myBreaks)
pdf(file = "NM/hp-2.pdf",width = 4,height = 8)
add.flag(p1,
         kept.labels = gene_name,
         repel.degree = 0.2)
dev.off()
image.png

这个图和开头提到的论文里的Figure3f就有几分相似了,但是还没有添加分组信息

欢迎大家关注我的公众号
小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享:1、R语言和python做数据分析和数据可视化的简单小例子;2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记;3、生物信息学入门学习资料及自己的学习笔记!

相关文章

网友评论

    本文标题:跟着Nature Methods学画图:R语言画热图(pheat

    本文链接:https://www.haomeiwen.com/subject/xldxfltx.html