ggplot2实现分半小提琴图绘制基因表达谱和免疫得分

生信宝典

发布于 2021-04-15 10:20:43

3.5K0

发布于 2021-04-15 10:20:43

文章被收录于专栏：生信宝典

最近看到很多人问下面这个图怎么绘制，看着确实不错。于是我查了一些资料，这个图叫split violin或者half violin，本质上是一种小提琴图。参考代码在https://gist.github.com/Karel-Kroeze/746685f5613e01ba820a31e57f87ec87

这里利用上期处理好的TCGA HNSCC的配对数据进行练习，数据包含43个肿瘤样本和43个癌旁样本。

除了基因表达量绘制的结果展示，最后还附带一个ESTIMATE计算免疫评分的例子。此外，计算免疫浸润主流的方法还有Cibersort、ssGSEA等算法，在之后的推文里我会做一些教程介绍。

具体代码如下:

remove(list = ls()) #一键清空
#加载包
library(ggplot2)
library(reshape2)
library(plyr)

suppressMessages(library(ggpubr))
suppressMessages(library(dplyr))

读入Deseq2标准化后的表达数据

# 1.1 表达数据
data <- read.csv("./Rawdata/TCGA_HNSCpaired_Norexpr_data_paired.csv",
                 header = T,row.names = 1, check.names=F)
data[1:3,1:4]

##       TCGA.CV.6943.01 TCGA.CV.6959.01 TCGA.CV.7438.11 TCGA.CV.7242.11
## AAAS        10.795406       11.198490       10.864833        10.89324
## AACS        12.033001       11.427610       12.173660        11.37317
## AADAC        4.857712        4.740191        8.864625        11.14658

保存为Tab键分割的格式，供estimate包使用。

# 安装ImageGP 包
# library(devtools)
# install_git("https://gitee.com/ct5869/ImageGP.git")
library(ImageGP)
sp_writeTable(data, file="./Rawdata/TCGA_HNSCpaired_Norexpr_data_paired.tsv")

筛选要绘制的基因

selected_gene <- c("S100A9","MT-CO2","MT-CO3","MT-CO1","KRT16","S100A8","KRT14","MT-CYB")

data <- data[selected_gene,]

读入分组数据

# 1.2 分组数据
phenotype <- read.csv("./Rawdata/TCGA_HNSCpaired.metadata.csv",header = T,row.names = 1,check.names=F, stringsAsFactors = T)
phenotype <- phenotype[,c("group","subject"),drop=F]
head(phenotype)

##                    group      subject
## TCGA.CV.6943.01    Tumor TCGA.CV.6943
## TCGA.CV.6959.01    Tumor TCGA.CV.6959
## TCGA.CV.7438.11 Nontumor TCGA.CV.7438
## TCGA.CV.7242.11 Nontumor TCGA.CV.7242
## TCGA.CV.7432.01    Tumor TCGA.CV.7432
## TCGA.CV.6939.11 Nontumor TCGA.CV.6939

绘图数据的格式准备

## 
# 2.1 把分组信息加进去
data_new <- data.frame(t(data))
data_new$sample = row.names(data_new)
data_new <- merge(data_new, phenotype,by.x = "sample",by.y = 0)

# 2.2 融合数据
data_new = melt(data_new)

## Using sample, group, subject as id variables

colnames(data_new) = c("sample","group","subject", "gene","expression")

head(data_new)

##            sample    group      subject   gene expression
## 1 TCGA.CV.6933.01    Tumor TCGA.CV.6933 S100A9   17.30560
## 2 TCGA.CV.6933.11 Nontumor TCGA.CV.6933 S100A9   10.83050
## 3 TCGA.CV.6934.01    Tumor TCGA.CV.6934 S100A9   13.57274
## 4 TCGA.CV.6934.11 Nontumor TCGA.CV.6934 S100A9   16.36757
## 5 TCGA.CV.6935.01    Tumor TCGA.CV.6935 S100A9   12.29116
## 6 TCGA.CV.6935.11 Nontumor TCGA.CV.6935 S100A9   20.22816

加载绘图函数

## 3. 这里加载包装好的2个函数,用于后面的统计和绘图
source("./assist/Function_for_violin_plot.R")

计算均值和误差

## 4. 绘图
# 4.1 这里注意到原图用的是误差线，这里用步骤三加载的函数，计算一下误差信息
Data_summary <- summarySE(data_new, measurevar="expression", groupvars=c("group","gene"))
head(Data_summary)

##      group   gene  N expression        sd        se        ci
## 1 Nontumor S100A9 43   18.02898 2.8049172 0.4277459 0.8632261
## 2 Nontumor MT.CO2 43   17.76159 0.8913531 0.1359301 0.2743180
## 3 Nontumor MT.CO3 43   18.18671 0.9300174 0.1418263 0.2862171
## 4 Nontumor MT.CO1 43   19.27698 1.0075457 0.1536493 0.3100768
## 5 Nontumor  KRT16 43   15.84541 3.5569291 0.5424266 1.0946612
## 6 Nontumor S100A8 43   17.17403 3.4525920 0.5265153 1.0625510

绘制分半小提琴图

# 4.2. 出图
  # 这个是我自己写的一个ggplot2的主题，可以自定义修改其中的参数
if(T){
  mytheme <- theme(plot.title = element_text(size = 12,color="black",hjust = 0.5),
                       axis.title = element_text(size = 12,color ="black"), 
                       axis.text = element_text(size= 12,color = "black"),
                       panel.grid.minor.y = element_blank(),
                       panel.grid.minor.x = element_blank(),
                       axis.text.x = element_text(angle = 45, hjust = 1 ),
                       panel.grid=element_blank(),
                       legend.position = "top",
                       legend.text = element_text(size= 12),
                       legend.title= element_text(size= 12)
) 
}

  # 自行调整下面的参数
gene_split_violin <- ggplot(data_new,aes(x= gene,y= expression,fill= group))+
  geom_split_violin(trim= F,color="white",scale = "area") + #绘制分半的小提琴图
  geom_point(data = Data_summary,aes(x= gene, y= expression),pch=19,
             position=position_dodge(0.5),size= 1)+ #绘制均值为点图
  geom_errorbar(data = Data_summary,aes(ymin = expression-ci, ymax= expression+ci), 
                width= 0.05, 
                position= position_dodge(0.5), 
                color="black",
                alpha = 0.8,
                size= 0.5) +
  scale_fill_manual(values = c("#56B4E9", "#E69F00"))+ 
  labs(y=("Log2 expression"),x=NULL,title = "Split violin") + 
  theme_bw()+ mytheme +
  stat_compare_means(aes(group = group),
                     label = "p.signif",
                     method = "anova",
                     label.y = max(data_new$expression),
                      hide.ns = T)
gene_split_violin;ggsave(gene_split_violin,
                         filename = "./Output/gene_split_violin.pdf",
                         height = 10,width = 16,units = "cm")

用ESTIMATE算法计算免疫得分

library(estimate)

# eestimate 包安装
# library(utils)
# rforge <- "http://r-forge.r-project.org"
# install.packages("estimate", repos=rforge, dependencies=TRUE)

# 5. ESTIMATE计算免疫得分
# 5.1 输入txt格式的表达矩阵，输出ESIMATE计算结果
filterCommonGenes(input.f= "./Rawdata/TCGA_HNSCpaired_Norexpr_data_paired.tsv", 
                  output.f="./Output/TCGA_estimate.gct", id="GeneSymbol")

## [1] "Merged dataset includes 9219 genes (1193 mismatched)."

estimateScore(input.ds = "./Output/TCGA_estimate.gct",
              output.ds = "./Output/TCGA_estimate_score.gct", 
              platform="illumina")

## [1] "1 gene set: StromalSignature  overlap= 135"
## [1] "2 gene set: ImmuneSignature  overlap= 139"

ESTI_score <- read.table("./Output/TCGA_estimate_score.gct",skip = 2,header = T,row.names = 1)
ESTI_score <- as.data.frame(t(ESTI_score[2:ncol(ESTI_score)]))
head(ESTI_score)

##                 StromalScore ImmuneScore ESTIMATEScore
## TCGA.CV.6943.01     906.2923  1649.01369    2555.30599
## TCGA.CV.6959.01    -352.6656   318.22117     -34.44448
## TCGA.CV.7438.11   -1183.4705   276.89782    -906.57268
## TCGA.CV.7242.11   -1067.1461   -47.83809   -1114.98415
## TCGA.CV.7432.01   -1234.5253  -449.47317   -1683.99846
## TCGA.CV.6939.11     424.8381  -674.84391    -250.00579

# 5.2 融合数据
table(row.names(ESTI_score) == rownames(phenotype))

## 
## TRUE 
##   86

ESTI_score$group <- phenotype$group
ESTI_score$sample <- rownames(ESTI_score)

ESTI_score_New =  melt(ESTI_score)

## Using group, sample as id variables

colnames(ESTI_score_New)=c("group","sample","status","score")  #设置行名
head(ESTI_score_New)

##      group          sample       status      score
## 1    Tumor TCGA.CV.6943.01 StromalScore   906.2923
## 2    Tumor TCGA.CV.6959.01 StromalScore  -352.6656
## 3 Nontumor TCGA.CV.7438.11 StromalScore -1183.4705
## 4 Nontumor TCGA.CV.7242.11 StromalScore -1067.1461
## 5    Tumor TCGA.CV.7432.01 StromalScore -1234.5253
## 6 Nontumor TCGA.CV.6939.11 StromalScore   424.8381

# 5.3 计算误差线
ESTI_Data_summary <- summarySE(ESTI_score_New, measurevar="score", groupvars=c("group","status"))
head(ESTI_Data_summary)

##      group        status  N       score        sd        se       ci
## 1 Nontumor  StromalScore 43  -918.91119  837.0492 127.64881 257.6057
## 2 Nontumor   ImmuneScore 43  -210.02121  504.0363  76.86481 155.1195
## 3 Nontumor ESTIMATEScore 43 -1128.93240 1138.1268 173.56271 350.2637
## 4    Tumor  StromalScore 43  -517.35577  659.5859 100.58590 202.9906
## 5    Tumor   ImmuneScore 43   -67.37634  638.1790  97.32139 196.4025
## 6    Tumor ESTIMATEScore 43  -584.73211 1138.8398 173.67145 350.4832

ESTI_split_violin <- ggplot(ESTI_score_New,aes(x= status,y= score,fill= group))+
  geom_split_violin(trim= F,color="white",scale = "area") + #绘制分半的小提琴图
  geom_point(data = ESTI_Data_summary,aes(x= status, y= score),pch=19,
             position=position_dodge(0.4),size= 1)+ #绘制均值为点图
  geom_errorbar(data = ESTI_Data_summary,aes(ymin = score-ci, ymax= score+ci), 
                width= 0.05, 
                position= position_dodge(0.4), 
                color="black",
                alpha = 0.8,
                size= 0.5) +
  scale_fill_manual(values = c("#56B4E9", "#E69F00"))+ 
  labs(y=("ESTIMATE score"),x=NULL,title = "Split violin") + 
  theme_bw()+ mytheme +
  scale_x_discrete(labels=c("Stromal","Immune","ESTIMATE")) +
  stat_compare_means(aes(group = group),
                     label = "p.signif",
                     method = "wilcox",
                     label.y = max(ESTI_score_New$score),
                      hide.ns = T)
ESTI_split_violin; ggsave(ESTI_split_violin,filename = "./Output/ESTIMATE_plot.pdf", height = 10,width = 10,units = "cm")

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 14393)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.936  LC_CTYPE=Chinese (Simplified)_China.936   
## [3] LC_MONETARY=Chinese (Simplified)_China.936 LC_NUMERIC=C                              
## [5] LC_TIME=Chinese (Simplified)_China.936    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ImageGP_0.1.0   devtools_2.3.0  usethis_1.6.1   dplyr_1.0.0     ggpubr_0.4.0   
## [6] estimate_1.0.13 plyr_1.8.6      reshape2_1.4.4  ggplot2_3.3.2  
## 
## loaded via a namespace (and not attached):
##  [1] matrixStats_0.56.0   fs_1.4.2             RColorBrewer_1.1-2   rprojroot_1.3-2     
##  [5] rstan_2.21.1         tools_4.0.2          backports_1.1.7      R6_2.4.1            
##  [9] colorspace_1.4-1     withr_2.2.0          tidyselect_1.1.0     gridExtra_2.3       
## [13] prettyunits_1.1.1    processx_3.4.3       curl_4.3             compiler_4.0.2      
## [17] git2r_0.27.1         cli_2.0.2            desc_1.2.0           labeling_0.3        
## [21] scales_1.1.1         callr_3.4.3          stringr_1.4.0        digest_0.6.25       
## [25] StanHeaders_2.21.0-5 foreign_0.8-80       rmarkdown_2.3        rio_0.5.16          
## [29] htmltools_0.5.1.1    pkgconfig_2.0.3      sessioninfo_1.1.1    rlang_0.4.6         
## [33] readxl_1.3.1         rstudioapi_0.11      farver_2.0.3         generics_0.1.0      
## [37] jsonlite_1.7.0       zip_2.1.1            car_3.0-8            inline_0.3.15       
## [41] magrittr_1.5         loo_2.3.1            Rcpp_1.0.5           munsell_0.5.0       
## [45] fansi_0.4.1          abind_1.4-5          lifecycle_0.2.0      stringi_1.4.6       
## [49] yaml_2.2.1           carData_3.0-4        pkgbuild_1.1.0       grid_4.0.2          
## [53] parallel_4.0.2       forcats_0.5.0        crayon_1.3.4         haven_2.3.1         
## [57] hms_0.5.3            knitr_1.29           ps_1.3.3             pillar_1.4.6        
## [61] ggsignif_0.6.0       codetools_0.2-16     stats4_4.0.2         pkgload_1.1.0       
## [65] glue_1.4.1           evaluate_0.14        V8_3.2.0             data.table_1.14.0   
## [69] remotes_2.1.1        BiocManager_1.30.10  RcppParallel_5.0.2   vctrs_0.3.1         
## [73] testthat_2.3.2       cellranger_1.1.0     gtable_0.3.0         purrr_0.3.4         
## [77] tidyr_1.1.0          assertthat_0.2.1     xfun_0.15            openxlsx_4.1.5      
## [81] broom_0.7.0          rstatix_0.6.0        tibble_3.0.2         pheatmap_1.0.12     
## [85] memoise_1.1.0        ellipsis_0.3.1