R可视化：组间微生物代谢通路

前言

组间微生物富集代谢通路常用reporter score表示，在得到reporter score结果后，我们常需要可视化其结果。更多知识分享请到 https://zouhua.top/。

Importing data

library(dplyr)
library(tibble)
library(ggplot2)
library(xlsx)

grp <- c("ET_B", "ET_P")
grp.col <- c("#FF6000", "#0080C0")
stg <- c("Before", "After")
stg.col <- c("#6288BA", "#C25D97")

phen <- read.csv("phenotype.csv") 
stg.path <- read.xlsx("reporter_score.20181114.xlsx", sheetIndex = 1)
bp.path <- read.xlsx("reporter_score.20181114.xlsx", sheetIndex = 3)

Function

ScoreVisuStage <- function(x, cfg, num=1){
  idx <- c("Pathway", "Level2", "Description")
  
  if (num==1) {
    level <- c("Amino_Acid_Metabolism","Carbohydrate_Metabolism",
               "Energy_Metabolism", 
               "Membrane_Transport","Xenobiotics_Biodegradation_and_Metabolism",
               "Metabolism_of_Cofactors_and_Vitamins")
  } else {
    level <- c("Amino_Acid_Metabolism","Carbohydrate_Metabolism",
               "PTS","Xenobiotics_Biodegradation_and_Metabolism",
               "Membrane_Transport",
               "Metabolism_of_Cofactors_and_Vitamins", "Energy_Metabolism")
  }
  
  if (cfg=="Before"){
    tmp <- x[ ,c(1:4)] %>% setNames(c(idx, "Reportscore"))
    
  } else {
    tmp <- x[ ,c(1:3, 5)] %>% setNames(c(idx, "Reportscore"))
  }
  
  dat <- tmp %>% filter(Level2 %in% level) %>% 
               filter(Reportscore < -1.96 | Reportscore > 1.96)
  dat.z <- dat %>% select(c(idx, grep("score", colnames(dat)))) %>%
           setNames(c(idx, "ET_B vs ET_P")) %>%
           tidyr::gather(Group, Reporter, -idx) %>%
           mutate(newdir=ifelse(Reporter < 0, grp[2], grp[1])) %>%
           mutate(Group=factor(Group, 
             levels = c("ET_B vs ET_P"))) 
  
  ord_aa <- dat.z %>% filter(Group %in% "ET_B vs ET_P") %>%
            arrange(Reporter, Level2) 
  
  dat.z$Description <- factor(dat.z$Description, levels = ord_aa$Description)
  dat.z$Level2 <- factor(dat.z$Level2, levels = level)
  dat.z$newdir <- factor(dat.z$newdir, levels=grp)  
  return(
    ggplot(dat.z, aes(y=Description, x=Reporter, color=newdir))+ 
    geom_segment(xend=0, aes(yend=Description), size=3)+
    geom_vline(xintercept = 0)+
    geom_vline(xintercept = c(-1.96,1.96), linetype="dashed", color="grey50")+
    facet_grid(Level2 ~ Group, scales = "free", space = "free_y")+
    theme_bw()+
    scale_color_manual(values = rev(grp.col),
                       labels = grp)+
    labs(x="Reporter score", y="")+
    scale_x_continuous(breaks=seq(-6, 6, 2),
                       limits=c(-6.1, 6.1))+
    theme(panel.grid.major = element_blank(),
          panel.grid.minor.x = element_blank(),
          legend.position = "NA",
          strip.text.y = element_blank(),
          axis.title = element_text(size=9, color="black", face="bold"),
          strip.text = element_text(size=10, color="black", face="bold"),
          axis.text = element_text(size=8, color="black"))
  )
}

# groups 
ScoreVisuBP <- function(x, cfg, num=1){
  idx <- c("Pathway", "Level2", "Description")
  
  if (num==1) {
    level <- c("Amino_Acid_Metabolism","Carbohydrate_Metabolism",
               "Energy_Metabolism", 
               "Membrane_Transport","Xenobiotics_Biodegradation_and_Metabolism",
               "Metabolism_of_Cofactors_and_Vitamins")
  } else {
    level <- c("Amino_Acid_Metabolism","Carbohydrate_Metabolism",
               "PTS","Xenobiotics_Biodegradation_and_Metabolism",
               "Membrane_Transport",
               "Metabolism_of_Cofactors_and_Vitamins", "Energy_Metabolism")
  }
  
  if (cfg=="B"){
    tmp <- x[ ,c(1:4)] %>% setNames(c(idx, "Reportscore"))
    
  } else {
    tmp <- x[ ,c(1:3, 5)] %>% setNames(c(idx, "Reportscore"))
  }
  
  dat <- tmp %>% filter(Level2 %in% level) %>% 
               filter(Reportscore < -1.96 | Reportscore > 1.96)
  dat.z <- dat %>% select(c(idx, grep("score", colnames(dat)))) %>%
           setNames(c(idx, "Baseline vs After")) %>%
           tidyr::gather(Group, Reporter, -idx) %>%
           mutate(newdir=ifelse(Reporter < 0, stg[1], stg[2])) %>%
           mutate(Group=factor(Group, 
             levels = c("Baseline vs After"))) 
  
  ord_aa <- dat.z %>% filter(Group %in% "Baseline vs After") %>%
            arrange(Reporter, Level2) 
  
  dat.z$Description <- factor(dat.z$Description, levels = ord_aa$Description)
  dat.z$Level2 <- factor(dat.z$Level2, levels = level)
  dat.z$newdir <- factor(dat.z$newdir, levels=stg)
  return(
    ggplot(dat.z, aes(y=Description, x=Reporter, color=newdir))+ 
    geom_segment(aes(yend=Description),xend=0, size=3)+
    geom_vline(xintercept = 0)+
    geom_vline(xintercept = c(-1.96,1.96), linetype="dashed", color="grey50")+
    facet_grid(Level2 ~ Group, scales = "free", space = "free_y")+
    theme_bw()+
    labs(x="Reporter score",
         y="")+
    scale_x_continuous(breaks=seq(-6, 6, 2),
                       limits=c(-6.1, 6.1))+
    scale_color_manual(values = stg.col,
                       labels = stg)+  
    theme(panel.grid.major = element_blank(),
          panel.grid.minor.x = element_blank(),
          legend.position = "NA",
          strip.text.y = element_blank(),
          axis.title = element_text(size=9, color="black", face="bold"),
          strip.text = element_text(size=10, color="black", face="bold"),
          axis.text = element_text(size=8, color="black"))
  )
}

Run

• Pathway in Stage

ScoreVisuStage(stg.path, "Before", 1)

• Pathway in Group

ScoreVisuBP(bp.path, "B", 1)

参考

参考文章如引起任何侵权问题，可以与我联系，谢谢。