R语言代码-Nature主刊中相关系系数热图绘制

需求描述

计算多个基因和免疫浸润的相关性，画出类似于文章中这种热图，用颜色展示相关系数，用**标出显著的。

example.png

出自https://www.nature.com/articles/s41586-019-1032-7

应用场景

评价基因对免疫浸润的影响，甚至可以批量筛选出影响免疫浸润的候选基因。

其本质是计算两种特征之间的相关性，此处的免疫浸润可以换成其他变量，例如临床表型等等。

如果要同时展示两个基因之间的相关性，及其与免疫浸润之间的关系，请参考“FigureYa96R2”。

输入文件

要求前两个文件的样本名一致

• ssGSEA_output.csv，免疫细胞矩阵，列是免疫细胞，行是样本，由FigureYa71ssGSEA产生
• easy_input_expr.csv，基因表达矩阵，列是样本，行是基因。行也可以是转录本，甚至是临床信息。
• genelist.txt，基因名，将要计算这些基因跟免疫浸润的相关性。要求跟easy_input_expr.csv里的基因名一致。

### 读入ssGSEA的结果
tcga_gsva <- read.csv("ssGSEA_output.csv",row.names = 1)
rownames(tcga_gsva) <- gsub("\\.","-",rownames(tcga_gsva))
head(tcga_gsva)[,1:3]

### 读入表达量数据
tcga_expr <- data.table::fread("easy_input_expr.csv", data.table = F, )
## 第一列转为行名
rownames(tcga_expr) <- tcga_expr[,1]
tcga_expr <- tcga_expr[,-1]
## 调整免疫矩阵中的样本顺序
tcga_gsva <- tcga_gsva[colnames(tcga_expr),]
tcga_expr[1:3,1:4]

### 读入感兴趣的基因
genelist <- read.table("genelist.txt")
genelist <- genelist$V1

批量计算相关性

先写一个函数
输入一个基因，即可返回跟免疫基因的相关性、p值

gene <- genelist
immuscore <- function(gene){
  y <- as.numeric(tcga_expr[gene,])
  colnames <- colnames(tcga_gsva)
  do.call(rbind,lapply(colnames, function(x){
    dd  <- cor.test(as.numeric(tcga_gsva[,x]),y,type="spearman")
    data.frame(gene=gene,immune_cells=x,cor=dd$estimate,p.value=dd$p.value )
  }))
}

#以FOXP3为例，测试一下函数
immuscore("FOXP3")

批量计算genelist跟免疫浸润相关性的结果

data <- do.call(rbind,lapply(genelist,immuscore))
head(data)

#保存到文件
write.csv(data, "correlation.csv", quote = F, row.names = F)

增加一列，区分p值的大小

使用两个ifelse实现三分类

data$pstar <- ifelse(data$p.value < 0.05,
                     ifelse(data$p.value < 0.01,"**","*"),
                     "")
data$pstar[1:20]
[1] ""   ""   "**" "**" "**" ""   "**" "*"  ""   ""   "**" "**" "**" "**" ""   "**" "**" "**" "*"  "**"

开始画图

使用ggplot2画图，主要用到的是geom_tile函数：

• 相关性用颜色的不同来表示，相关性的大小用颜色的深浅来反映；

• 有差异的把*号打印在热图上

library(ggplot2)
library(dplyr)
ggplot(data, aes(immune_cells, gene)) + 
  geom_tile(aes(fill = desc(cor)), colour = "white",size=1)+
  scale_fill_gradient2(low = "#2b8cbe",mid = "white",high = "#e41a1c")+
  geom_text(aes(label=pstar),col ="black",size = 5)+
  theme_minimal()+# 不要背景
  theme(axis.title.x=element_blank(),#不要title
        axis.ticks.x=element_blank(),#不要x轴
        axis.title.y=element_blank(),#不要y轴
        axis.text.x = element_text(angle = 45, hjust = 1),# 调整x轴文字
        axis.text.y = element_text(size = 8))+#调整y轴文字
  #调整legen
  labs(fill =paste0(" * p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"))

#保存到文件
ggsave("correlation.pdf", width = 8, height = 4)

correlation.png

geom_tile函数中的colour可以改变颜色，例如black

library(ggplot2)
library(dplyr)
ggplot(data, aes(immune_cells, gene)) + 
  geom_tile(aes(fill = desc(cor)), colour = "black",size=1)+
  scale_fill_gradient2(low = "#2b8cbe",mid = "white",high = "#e41a1c")+
  geom_text(aes(label=pstar),col ="black",size = 5)+
  theme_minimal()+# 不要背景
  theme(axis.title.x=element_blank(),#不要title
        axis.ticks.x=element_blank(),#不要x轴
        axis.title.y=element_blank(),#不要y轴
        axis.text.x = element_text(angle = 45, hjust = 1),# 调整x轴文字
        axis.text.y = element_text(size = 8))+#调整y轴文字
  #调整legen
  labs(fill =paste0(" * p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"))

correlation2.png

geom_tile函数中的colour可以改变颜色，例如black

如果colour缺失，就没有边框的颜色

library(ggplot2)
library(dplyr)
ggplot(data, aes(immune_cells, gene)) + 
  geom_tile(aes(fill = desc(cor)),size=1)+
  scale_fill_gradient2(low = "#2b8cbe",mid = "white",high = "#e41a1c")+
  geom_text(aes(label=pstar),col ="black",size = 5)+
  theme_minimal()+# 不要背景
  theme(axis.title.x=element_blank(),#不要title
        axis.ticks.x=element_blank(),#不要x轴
        axis.title.y=element_blank(),#不要y轴
        axis.text.x = element_text(angle = 45, hjust = 1),# 调整x轴文字
        axis.text.y = element_text(size = 8))+#调整y轴文字
  #调整legend
  labs(fill =paste0(" * p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"))

correlation3.png

如果把其他变量映射到size参数，可以实现热图小方块的大小改变
本次不推荐，我觉得不好看

library(ggplot2)
library(dplyr)
ggplot(data, aes(immune_cells, gene)) + 
  geom_tile(aes(fill = desc(cor),size=desc(cor)), colour = "white")+
  scale_fill_gradient2(low = "#2b8cbe",mid = "white",high = "#e41a1c")+
  scale_size_continuous(range =c(0,6))+
  geom_text(aes(label=pstar),col ="black",size = 5)+
  theme_minimal()+
  theme(axis.title.x=element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y = element_text(size = 8))+
  labs(fill =paste0("* p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"),size="Correlation")

correlation4.png