最近看到很多人問下面這個圖怎么繪制，看著確實不錯。于是我查了一些資料，這個圖叫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 variablescolnames(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 ## 86ESTI_score$group <- phenotype$group ESTI_score$sample <- rownames(ESTI_score)ESTI_score_New = melt(ESTI_score)## Using group, sample as id variablescolnames(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.4832ESTI_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

參考資料：

《數據可視化——R語言ggplot2包繪制精美的小提琴圖》

數據和代碼下載：

https://gitee.com/ct5869/shengxin-baodian/tree/master/TCGA

作者：趙法明

編輯：生信寶典

往期精品(點擊圖片直達文字對應教程)

機器學習

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的ggplot2实现分半小提琴图绘制基因表达谱和免疫得分的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。