当前位置：首页 > news >正文

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

news 2025/7/22 22:37:06

介绍

复现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'))+ #AUCgeom_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())