复现GMM文章（五）：图5代码和数据

介绍

复现GMM文章的的Fig5图。

加载R包

library(tidyr)
library(tidyverse)
library(dplyr)
library(ggsci)
library(ggpubr)
library(psych)
library(ggcorrplot)
library(showtext)
library(ggrepel)

导入数据

所有的数据可以通过下列链接下载：

百度网盘链接: https://pan.baidu.com/s/1isKEK1G5I6X90KYqLufmWw

提取码: t9ca

图5A

• 数据

load('01_data/plot_data/F5A.RData')

head(auc_marker)

cor <- corr.test(auc_marker[,1:2],method = 'spearman')

• 画图

p1 <- ggscatter(data=auc_marker,x='dis',y='auc',shape='level',color='disease',palette = "igv",alpha=0.8,size=2)+
  geom_smooth(mapping=aes(x=dis,y=auc),data=auc_marker,formula =y~x, method=lm, inherit.aes = F,size=0.5)+
  geom_text_repel(data=auc_marker,aes(x=dis,y=auc,label = disease),inherit.aes = F,size=3)+
  xlab("MSI (median)") + ylab("External AUC (median)")+
  labs(color = "disease")+
  theme(legend.position = 'none')+
  annotate('text',x=0.7,y=0.84,label=paste0('r=',round(cor$r[1,2],2),'  ','p','=',format(cor$p[1,2],scientific=T,digits=3)),size=3)+
  border()

p2 <- ggdensity(auc_marker, 'dis', fill="disease_type2",palette = "npg")+
  clean_theme()+
  theme(legend.position = 'none')
p3 <- ggdensity(auc_marker, 'auc', fill="disease_type2", palette = "npg")+
  clean_theme()+
  theme(legend.position = 'none')+
  ggpubr::rotate()

p <- ggarrange(p2, NULL, p1, p3, ncol = 2, nrow = 2, align = "hv", widths = c(3, 1), 
            heights = c(1, 3))

p

图5B

• 数据

load('01_data/plot_data/F5B.RData')

head(self.e)

stat.test <- compare_means(
  dis_value~group1,data = self.e, group.by = "level",
  method = "wilcox.test") %>% 
  mutate(y.position = c(seq(from=4, to=7.5,length.out=10),rep(seq(from=4, to=6,length.out=6), times=2)))
x <- stat.test$p.adj
stat.test$p.adj.signif<-ifelse(x<0.05, ifelse(x<0.01, ifelse(x<0.001, ifelse(x<=0.0001, '****','***'),'**'),'*'),'ns')

• 画图

p <- ggboxplot(self.e, x = "group1", y = "dis_value",fill = "group1", palette = "jco", 
               facet.by = "level",width = 0.5,outlier.shape = NA)+
  stat_compare_means(aes(label = paste0('p = ',..p.format..)),label.y = 6.9,)+
  theme(axis.text.x = element_blank(),axis.ticks=element_blank(),legend.position="none")+
  xlab("Disease category") + ylab("MSI")+
  theme(text = element_text(size=13,face = 'plain',family ='',colour = 'black')) + 
  stat_pvalue_manual(stat.test,label = "p.adj.signif")

p

图5C

• 数据

load('01_data/plot_data/F5C.RData')

head(self.e)

stat.test <- compare_means(
  dis_value~level,data = self.e, group.by = "subtype",
  method = "wilcox.test") %>% 
  mutate(y.position = rep(seq(from=4.5, to=5.5,length.out=3), times=4))
x <- stat.test$p.adj
stat.test$p.adj.signif <- ifelse(x<0.05, ifelse(x<0.01, ifelse(x<0.001, ifelse(x<=0.0001, '****','***'),'**'),'*'),'ns')

• 画图

p <- ggboxplot(self.e, x = "level", y = "dis_value",fill = "level", 
               palette = c('#774ec7','#bd93cc','#a2c4b1'),
               facet.by = "subtype",nrow = 1,width = 0.5,outlier.shape = NA)+
  stat_compare_means(aes(label = paste0('p = ',..p.format..)),label.y = 6.9,)+
  theme(axis.text.x = element_blank(),axis.ticks=element_blank(),legend.position="top")+ 
  xlab("Data type") + ylab("MSI")+
  theme(text = element_text(size=13,face = 'plain',family ='',colour = 'black'),legend.position = 'none') + 
  stat_pvalue_manual(stat.test,label = "p.adj.signif")

p

图5D

• 数据

load('01_data/plot_data/F5D.RData')

head(auc_marker_lodo)

• 画图

p <- ggpaired(auc_marker_lodo, x = 'method', y = 'dis',color = 'method', palette = "npg", 
              line.color = "gray", line.size = 0.4,
              short.panel.labs = FALSE,width = 0.4)+
  stat_compare_means(paired = T,label.y=4.6)+
  theme(legend.position = 'none')+
  labs(y='MSI (median)',x='')+
  theme(text = element_text(family = '',size = 13,face = 'plain'))

p

图5E

左侧

• 数据

load('01_data/plot_data/F5E_left.RData')

head(p_r_test)

for (i in unique(p_r_test$num)){
  dis_median <- median(subset(p_r_test,num==i)$dis)
  if (i==1){
    p_r_test_median <- data.frame(dis=dis_median,method='test',num=i)
  }else{
    p_r_test_median <- rbind(p_r_test_median,data.frame(dis=dis_median,method='test',num=i))
  } 
}

• 画图

p1 <- ggboxplot( p_r_test,x='num',y='dis',color=pal_npg("nrc")(10)[1],add = "jitter",width = 0.3)+
  xlab('No. of training cohorts')+
  ylab('MSI (median)')+
  labs(title = 'ASD')+
  geom_line(data = p_r_test_median,aes(x=num,y=dis),color=pal_npg("nrc")(10)[2])+
  geom_point(data = p_r_test_median,aes(x=num,y=dis),color=pal_npg("nrc")(10)[2])+
  # geom_hline(yintercept =0.5,color=pal_d3('category10')(10)[4])+
  theme(text = element_text(size=13,face = 'plain',family ='',colour = 'black'))+  
  theme(panel.background = element_blank(),axis.line = element_line(colour = "black")) +
  ylim(0,2.5)

p1

右侧

• 数据

load('01_data/plot_data/F5E_right.RData')

for (i in unique(p_r_test$num)){
  dis_median <- median(subset(p_r_test,num==i)$dis)
  if (i==1){
    p_r_test_median <- data.frame(dis=dis_median,method='test',num=i)
  }else{
    p_r_test_median <- rbind(p_r_test_median,data.frame(dis=dis_median,method='test',num=i))
  } 
}

• 画图

p2 <- ggboxplot( p_r_test,x='num',y='dis',color=pal_npg("nrc")(10)[1],add = "jitter",width = 0.3)+
  xlab('No. of training cohorts')+
  ylab('MSI (median)')+
  labs(title = 'PD')+
  geom_line(data = p_r_test_median,aes(x=num,y=dis),color=pal_npg("nrc")(10)[2])+
  geom_point(data = p_r_test_median,aes(x=num,y=dis),color=pal_npg("nrc")(10)[2])+
  theme(text = element_text(size=13,face = 'plain',family ='',colour = 'black'))+  
  theme(panel.background = element_blank(),axis.line = element_line(colour = "black"))

p2

合并

pl <- ggarrange(p1,p2)
pl