R/CovariatePlot.R

In katrikorpela/mare: Microbiota Analysis in R Easily

CovariatePlot <- function(meta, taxonomic.table, covariate, taxa=NULL, smooth.method = "loess", 
    readcount.cutoff = 0, select.by = NULL, select = NULL, pdf = F, logtrans = F, nozeros = F, relative = T) {

   if(Sys.info()[['sysname']] == "Linux") {
  quartz <- function() {X11()}
  }
    if(Sys.info()[['sysname']] == "Windows") {
  quartz <- function() {X11()}
}
   
    meta <- read.delim(meta)
    if(!relative) meta[,"ReadCount"]<-1
    taxatable <- read.delim(taxonomic.table)
   
    names(taxatable) <- sapply(names(taxatable), function(x) gsub("_NA", ".", x))
    names(taxatable) <- sapply(names(taxatable), function(x) strsplit(x, split = "_")[[1]][length(strsplit(x, split = "_")[[1]])])
    if(length(taxa)==0) taxa <- colnames(taxatable)
    taxa <- sapply(taxa, function(x) gsub("_NA", ".", x))
    taxa <- sapply(taxa, function(x) strsplit(x, split = "_")[[1]][length(strsplit(x, split = "_")[[1]])])
    
    if(nozeros) taxatable[taxatable==0]<-NA
    
    if(logtrans) { dataset <- data.frame(meta, log(((taxatable+1)/meta$ReadCount)))
    } else dataset <- data.frame(meta, (taxatable/meta$ReadCount) * 100)

    if (length(select.by) != 0) {
        dataset$selection <- dataset[, select.by]
        dataset <- dataset[dataset$selection == select, ]
    }
    dataset <- dataset[dataset$ReadCount > readcount.cutoff, ]
    
   if(relative) lab <- "Relative abundance (%)" else lab <- "Abundance" 
    
    df = na.omit(reshape2::melt(dataset[, c(covariate, taxa)], id = c(covariate)))
    names(df) <- c("x", "variable", "value")
    p <- ggplot2::ggplot(df, ggplot2::aes(y = value, x = x), environment = environment()) + 
          ggplot2::stat_smooth(method = smooth.method, formula = y ~ x, se = T,color="gray30",  fill = "cornflowerblue") + 
          ggplot2::geom_point(pch=20,color="gray30")+
          ggplot2::facet_wrap(~variable, ncol = floor(sqrt(length(taxa))), scales = "free") + 
          ggplot2::theme_bw() + 
          ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
                         panel.border = ggplot2::element_blank(),
                         panel.grid.major = ggplot2::element_blank(), 
                         panel.grid.minor = ggplot2::element_blank(),
                         strip.background =  ggplot2::element_rect(color = "white",fill="white"))+
          ggplot2::xlab(covariate) + 
          ggplot2::ylab(lab) 
         
 quartz()
    plot(p)
    
    if (pdf) {
        pdf(paste(covariate,"Plot.pdf",sep=""))
       plot(p)
        dev.off()
    }
}