复现GMM文章（一）：图1代码和数据

介绍

复现GMM文章的的Fig1图。

加载R包

  library(tidyr)
  library(tidyverse)
  library(dplyr)
  library(ggsci)
  library(ggpubr)

导入数据

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

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

提取码: t9ca

图1B

• 数据

load("01_data/plot_data/F1B.RData")

head(temp)

• 画图

temp %>%
  mutate(type=factor(type, levels=c("Intestinal" ,"Metabolic" , "Mental"  ,  "Autoimmune", "Liver"))) %>%
  arrange(desc(count)) %>% 
  mutate(disease=factor(disease, levels=unique(disease))) %>%
  ggplot(aes(x=disease, y=count,group=data_type)) +

  geom_bar(stat="identity",position='stack', aes(fill=data_type)) +

  geom_text(aes(label=count),position=position_stack(vjust = 0.5),size=7)+
  facet_grid(~type, scales="free", space="free") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=45, hjust=1,face = 'bold',size=12),
        axis.text.y=element_text(face = 'bold',size=12),
        plot.title=element_text(hjust=0.5)) +
  ylab("No. of project") +
  xlab('disease') +
  coord_cartesian(ylim=c(0,11),expand=FALSE) +
  scale_y_continuous(breaks=seq(0, 12, 2))+
  theme(panel.border = element_blank(), axis.line = element_line())+
  scale_fill_d3(alpha = 0.5)+
  theme(text = element_text(size=16,face = 'plain',family ='',colour = 'black'))

图1C

• 数据

load("01_data/plot_data/F1C.RData")

head(project_stat0)

project_stat0 <- gather(project_stat0,phenotype,num,c('case','control'))
project_stat0$phenotype <- factor(project_stat0$phenotype,levels = c('control','case'))

• 画图

ggdensity(project_stat0, 'num', color="phenotype",palette = "aaas",add = "median",alpha = 0.1,size=1,fill ="phenotype",rug = TRUE)+
  labs(x = 'No. of samples in each cohort',y='Density')+
  annotate("text", label = paste0("Median: ",median(subset(project_stat0,phenotype=='case')$num)), x = 150, y = 0.015, size = 4, colour = pal_aaas("default", alpha = 0.6)(10)[2])+
  annotate("text", label = paste0("Median: ",median(subset(project_stat0,phenotype=='control')$num)), x = 150, y = 0.013, size = 4, colour = pal_aaas("default", alpha = 0.6)(10)[1])

图1D

• 数据

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

head(auc_self)

stat.test <- compare_means(
  auc~group1,data = auc_self, 
  # group.by = "level",
  method = "wilcox.test") %>% 
  mutate(y.position = seq(from=1.05, to=1.65,length.out=10))
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')

• 画图

ggboxplot(auc_self, x = "group1", y = "auc", fill = "group1",
               palette = "jco",width = 0.2)+ 
  geom_hline(yintercept =0.5,color='#dbdcdc')+
  geom_hline(yintercept =0.6,color='#ffd09a')+
  geom_hline(yintercept =0.7,color='#ffcbd8')+
  geom_hline(yintercept =0.8,color='#7b77ff')+
  geom_hline(yintercept =0.9,color='#e60020')+
  # stat_compare_means()+
  ylim(0.05,1.68)+
  theme(legend.position="none")+    
  ylab("Internal AUC")+xlab('')+
  ggtitle('Disease category')+
  theme(axis.text.x=element_text(angle=20, hjust=0.8,face = 'plain',size=13),
        text = element_text(size=13,face = 'plain',family ='',colour = 'black')) +
  stat_pvalue_manual(stat.test,label = "p.adj.signif")

图1E

• 数据

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

head(self.e)

stat.test <- compare_means(
  auc~level,data = self.e, 
  # group.by = "level",
  method = "wilcox.test") %>% 
  mutate(y.position = seq(from=1.2, to=1.65,length.out=3))
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')

• 画图

ggboxplot(self.e, x = "level", y = "auc", fill = "level",
               width = 0.2,palette = c('#774ec7','#bd93cc','#a2c4b1'))+
  geom_hline(yintercept =0.5,color='#dbdcdc')+
  geom_hline(yintercept =0.6,color='#ffd09a')+
  geom_hline(yintercept =0.7,color='#ffcbd8')+
  geom_hline(yintercept =0.8,color='#7b77ff')+
  geom_hline(yintercept =0.9,color='#e60020')+
  ylim(0.05,1.68)+
  # stat_compare_means()+
  theme(legend.position="none")+    
  ylab("Internal AUC")+xlab('')+
  ggtitle('Data type')+
  theme(axis.text.x=element_text(angle=20, hjust=0.8,face = 'plain',size=13),
        text = element_text(size=13,face = 'plain',family ='',colour = 'black')) + 
  stat_pvalue_manual(stat.test,label = "p.adj.signif")

图1F-G

• 数据

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

head(a_all)

• 画图

ggboxplot(a_all, x = "method", y = "auc", fill = "method",
               palette = c('#1fb8b4','#ff7f0e'),width = 0.15)+ 
  geom_hline(yintercept =0.5,color='#dbdcdc')+
  geom_hline(yintercept =0.6,color='#ffd09a')+
  geom_hline(yintercept =0.7,color='#ffcbd8')+
  geom_hline(yintercept =0.8,color='#7b77ff')+
  geom_hline(yintercept =0.9,color='#e60020')+
  facet_wrap(~group1,nrow = 1)+
  # annotate('text',x=1:2,y=0.15,label=c('0.765','0.638'))+ #AUC
  geom_signif(comparisons =list(c('internal','external')),
              y_position = c(1.12, 1.32),
              test = 'wilcox.test',
              map_signif_level = function(x){ifelse(x<0.05, 
                                                    ifelse(x<0.01, 
                                                           ifelse(x<0.001, 
                                                                  ifelse(x<=0.0001, '****','***'),'**'),'*'),'ns')})+
  ylim(0.05,1.32)+
  theme(legend.position="top")+    
  xlab("") + ylab("AUC")+
  labs(fill = "AUC type")+
  theme(text = element_text(size=13,face = 'plain',family ='',colour = 'black'),
        axis.text.x = element_blank(),
        axis.ticks=element_blank())