需求描述

計(jì)算多個(gè)基因和免疫浸潤(rùn)的相關(guān)性，畫(huà)出類似于文章中這種熱圖，用顏色展示相關(guān)系數(shù)，用**標(biāo)出顯著的。

example.png

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

應(yīng)用場(chǎng)景

評(píng)價(jià)基因?qū)γ庖呓?rùn)的影響，甚至可以批量篩選出影響免疫浸潤(rùn)的候選基因。

其本質(zhì)是計(jì)算兩種特征之間的相關(guān)性，此處的免疫浸潤(rùn)可以換成其他變量，例如臨床表型等等。

如果要同時(shí)展示兩個(gè)基因之間的相關(guān)性，及其與免疫浸潤(rùn)之間的關(guān)系，請(qǐng)參考“FigureYa96R2”。

輸入文件

要求前兩個(gè)文件的樣本名一致

• ssGSEA_output.csv，免疫細(xì)胞矩陣，列是免疫細(xì)胞，行是樣本，由FigureYa71ssGSEA產(chǎn)生
• easy_input_expr.csv，基因表達(dá)矩陣，列是樣本，行是基因。行也可以是轉(zhuǎn)錄本，甚至是臨床信息。
• genelist.txt，基因名，將要計(jì)算這些基因跟免疫浸潤(rùn)的相關(guān)性。要求跟easy_input_expr.csv里的基因名一致。

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

### 讀入表達(dá)量數(shù)據(jù)
tcga_expr <- data.table::fread("easy_input_expr.csv", data.table = F, )
## 第一列轉(zhuǎn)為行名
rownames(tcga_expr) <- tcga_expr[,1]
tcga_expr <- tcga_expr[,-1]
## 調(diào)整免疫矩陣中的樣本順序
tcga_gsva <- tcga_gsva[colnames(tcga_expr),]
tcga_expr[1:3,1:4]

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

批量計(jì)算相關(guān)性

先寫(xiě)一個(gè)函數(shù)
輸入一個(gè)基因，即可返回跟免疫基因的相關(guān)性、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為例，測(cè)試一下函數(shù)
immuscore("FOXP3")

批量計(jì)算genelist跟免疫浸潤(rùn)相關(guān)性的結(jié)果

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

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

增加一列，區(qū)分p值的大小

使用兩個(gè)ifelse實(shí)現(xiàn)三分類

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

開(kāi)始畫(huà)圖

使用ggplot2畫(huà)圖，主要用到的是geom_tile函數(shù)：

• 相關(guān)性用顏色的不同來(lái)表示，相關(guān)性的大小用顏色的深淺來(lái)反映；

• 有差異的把*號(hào)打印在熱圖上

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),# 調(diào)整x軸文字
        axis.text.y = element_text(size = 8))+#調(diào)整y軸文字
  #調(diào)整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函數(shù)中的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),# 調(diào)整x軸文字
        axis.text.y = element_text(size = 8))+#調(diào)整y軸文字
  #調(diào)整legen
  labs(fill =paste0(" * p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"))

correlation2.png

geom_tile函數(shù)中的colour可以改變顏色，例如black

如果colour缺失，就沒(méi)有邊框的顏色

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),# 調(diào)整x軸文字
        axis.text.y = element_text(size = 8))+#調(diào)整y軸文字
  #調(diào)整legend
  labs(fill =paste0(" * p < 0.05","\n\n","** p < 0.01","\n\n","Correlation"))

correlation3.png

如果把其他變量映射到size參數(shù)，可以實(shí)現(xiàn)熱圖小方塊的大小改變
本次不推薦，我覺(jué)得不好看

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