复现GMM文章（八）：附图7代码和数据

介绍

附图7的代码

library(psych)
library(ggcorrplot)
library(showtext)
library(ggpubr)
library(ggrepel)
library(tidyverse)
library(patchwork)
library(UpSetR)        
library(VennDiagram) 

Fig S7A

load('01_data/plot_data/FS7A.RData')
cor <- corr.test(auc_marker[,c('dis','auc')],method = 'spearman')
p1_1 <- 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=15,y=1,label=paste0('r=',round(cor$r[1,2],2),'  ','p','=',format(cor$p[1,2],scientific=T,digits=3)),size=4)+
  border() +ggtitle('ALDEx2')


p2_1 <- 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=40,y=1,label=paste0('r=',round(cor$r[1,2],2),'  ','p','=',format(cor$p[1,2],scientific=T,digits=3)),size=4)+
  border() +ggtitle('MaAslin2')


pA <- p1_1+p2_1

pA

Fig S7B

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

result_MSI$level <- paste0(result_MSI$data_type,'_',result_MSI$t)
result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'
result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Amplicon_genus'] <- '16S_genus'
result_MSI$level <- factor(result_MSI$level,levels=c("mNGS_species","mNGS_genus","16S_genus"))

result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'


stat.test <- compare_means(
  ALDEx2~group1,data = result_MSI, 
  # group.by = "level",
  method = "wilcox.test") %>%
  mutate(y.position = seq(from=60, to=100,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')

p3_1 <- ggboxplot(result_MSI, x = "group1", y = "ALDEx2", fill = "group1",
                  palette = "jco",width = 0.2,outlier.shape = NA)+ 
  theme_bw() +
  stat_compare_means()+
  ylim(0,110)+
  theme(legend.position="none")+    
  ylab("MSI (ALDEx2)")+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")
p3_1

result_MSI <- get(load('01_data/result_MSI_DA.RData'))
result_MSI <- subset(result_MSI,Group=='cutoff')
result_MSI$level <- paste0(result_MSI$data_type,'_',result_MSI$t)

result_MSI <- subset(result_MSI,disease %in% names(which(table(unique(result_MSI[,c('disease','data_type')])$disease)>1)))
result_MSI$group1_new <- factor(result_MSI$group1 , levels = c("Intestinal","Metabolic","Mental","Autoimmune"))

result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'
result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Amplicon_genus'] <- '16S_genus'
result_MSI$level <- factor(result_MSI$level,levels=c("mNGS_species","mNGS_genus","16S_genus"))

stat.test <- compare_means(
  ALDEx2~level,data = result_MSI, 
  # group.by = "level",
  method = "wilcox.test") %>%
  mutate(y.position = seq(from=70, to=90,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')

p4_1 <- ggboxplot(result_MSI, x = "level", y = "ALDEx2", fill = "level",
                  width = 0.2,palette = c('#774ec7','#bd93cc','#a2c4b1'),outlier.shape = NA)+
  ylim(0,110)+
  theme_bw() +
  stat_compare_means()+
  theme(legend.position="none")+    
  ylab("MSI (ALDEx2)")+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")
p4_1

pB <- ggarrange(p3_1,p4_1,
               ncol = 2, nrow = 1,
               widths = c(4,3))
pB

Fig S7C

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

result_MSI$level <- paste0(result_MSI$data_type,'_',result_MSI$t)
result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'
result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Amplicon_genus'] <- '16S_genus'
result_MSI$level <- factor(result_MSI$level,levels=c("mNGS_species","mNGS_genus","16S_genus"))

result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'

stat.test <- compare_means(
  ALDEx2~group1,data = result_MSI, 
  # group.by = "level",
  method = "wilcox.test") %>%
  mutate(y.position = seq(from=350, to=550,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')

p5_1 <- ggboxplot(result_MSI, x = "group1", y = "ALDEx2", fill = "group1",
                  palette = "jco",width = 0.2,outlier.shape = NA)+ 
  theme_bw() +
  stat_compare_means()+
  ylim(0,560)+
  theme(legend.position="none")+    
  ylab("MSI (MaAslin2)")+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")

p5_1


result_MSI <- get(load('01_data/result_MSI_DA.RData'))
result_MSI$ALDEx2 <- result_MSI$MaAslin2
result_MSI <- subset(result_MSI,Group=='cutoff')
result_MSI$level <- paste0(result_MSI$data_type,'_',result_MSI$t)

result_MSI <- subset(result_MSI,disease %in% names(which(table(unique(result_MSI[,c('disease','data_type')])$disease)>1)))
result_MSI$group1_new <- factor(result_MSI$group1 , levels = c("Intestinal","Metabolic","Mental","Autoimmune"))
result_MSI$level[result_MSI$level=='Metagenomics_species'] <- 'mNGS_species'
result_MSI$level[result_MSI$level=='Metagenomics_genus'] <- 'mNGS_genus'
result_MSI$level[result_MSI$level=='Amplicon_genus'] <- '16S_genus'
result_MSI$level <- factor(result_MSI$level,levels=c("mNGS_species","mNGS_genus","16S_genus"))

stat.test <- compare_means(
  ALDEx2~level,data = result_MSI, 
  # group.by = "level",
  method = "wilcox.test") %>%
  mutate(y.position = seq(from=350, to=500,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')

p6_1 <- ggboxplot(result_MSI, x = "level", y = "ALDEx2", fill = "level",
                  width = 0.2,palette = c('#774ec7','#bd93cc','#a2c4b1'),outlier.shape = NA)+
  ylim(0,560)+
  theme_bw() +
  stat_compare_means()+
  theme(legend.position="none")+    
  ylab("MSI (MaAslin2)")+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")

pC <- ggarrange(p5_1,p6_1,
               ncol = 2, nrow = 1,
               widths = c(4,3))
pC

Fig S7D

result_MSI <- get(load('01_data/result_MSI_DA.RData'))
result_MSI <- subset(result_MSI,Group=='cutoff')

r <- paste("cor: ",round(cor(result_MSI$LEfSe,result_MSI$ALDEx2,method = 'spearman'),3), sep = "")
p_v <- paste("p: ",format(cor.test(result_MSI$LEfSe,result_MSI$ALDEx2,method = 'spearman')$p.value,scientific=T,digits=3), sep = "")
title1 <- paste(r," ",p_v, sep = "")

r <- paste("cor: ",round(cor(result_MSI$LEfSe,result_MSI$MaAslin2,method = 'spearman'),3), sep = "")
p_v <- paste("p: ",format(cor.test(result_MSI$LEfSe,result_MSI$MaAslin2,method = 'spearman')$p.value,scientific=T,digits=3), sep = "")
title2 <- paste(r," ",p_v, sep = "")

p1 <- ggplot(result_MSI,mapping = aes(x=LEfSe, y=ALDEx2),size=1)+
  geom_point(size=2)+
  labs(x='MSI (LEfSe)',y='MSI (ALDEx2)')+
  color_palette('igv')+
  stat_smooth(mapping =aes(x=LEfSe, y=ALDEx2) ,method=lm,inherit.aes = F,size=0.5)+
  theme_classic()+
  annotate('text',x=0.5,y=70,label=title1,size=4)+
  border()+ggtitle('ALDEx2')

p2 <- ggplot(result_MSI,mapping = aes(x=LEfSe, y=MaAslin2),size=1)+
  geom_point(size=2)+
  labs(x='MSI (LEfSe)',y='MSI (MaAslin2)')+
  color_palette('igv')+
  stat_smooth(mapping =aes(x=LEfSe, y=MaAslin2) ,method=lm,inherit.aes = F,size=0.5)+
  theme_classic()+
  annotate('text',x=1,y=400,label=title2,size=4)+
  border() +ggtitle('MaAslin2')

pB <- p1+p2
pB

Fig S7E

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

pD <- upset(fromList(upset_list),  
            matrix.color = "gray23", main.bar.color = "#3B4992", sets.bar.color = "#3B4992",
            nsets = 10000,     
            nintersects = 10000, 
            order.by = "freq", 
            keep.order = F, 
            mb.ratio = c(0.6,0.4),   
            text.scale = 1.5, 
            mainbar.y.label = "Marker intersections", sets.x.label = "No. of markers per method")
pD