R/plotstab.R

In SAMPLE: Estimate Sampling Effort for Species Occurrence Rates

#' Plot
#'
#'This function generates the prevalence rate plots for all the variables in question (e.g. species, habitats).
#' @name plotstab
#' @usage plotstab(data, info, outputName = "Results", outputDir = tempdir())
#' @importFrom magrittr %>%
#' @importFrom dplyr if_else
#' @import ggplot2
#' @importFrom utils read.table write.table
#' @importFrom grDevices pdf dev.off
#' @param data A dataframe generated by the function RunPerm() (a dataframe object).
#' @param info A dataframe generated by the function stability() (a dataframe object).
#' @param outputName Prefix used for the output (a character; default output_N="Results").
#' @param outputDir Output directory (a character). If omitted, files are generated in the R session's temporary directory.
#' @return A text file (.txt) with the output values of the analysis, and a PDF file with the generated plots from the analysis.
#'
#' @examples
#' data("coral_symbionts")
#' set.seed(812)
#' perm = RunPerm(input = coral_symbionts,replicates = 50)
#' stable = stability(data = perm,stability_thresh = 5 ,success_points = 5,diff = 2 )
#' plotstab(data = perm, info = stable, outputName = "Stability_example")
#' unlink(file.path(tempdir(), "Stability_example.pdf"))
#' unlink(file.path(tempdir(), "Stability_example.txt"))
#' @export

# Declare global variables
utils::globalVariables(c("Substract", "Prevalence", "avg", "Taxa", "lci", "uci", "thres", "Host_sp", "colonies"))

plotstab <- function(data,info,outputName="Results",outputDir=tempdir()){

  ############################################
  # Set output directory defined by the user #
  ############################################
  # Ensure outputDir is valid
  outputDir <- ifelse(missing(outputDir) || outputDir == "", tempdir(), outputDir)
  dir.create(outputDir, showWarnings = FALSE, recursive = TRUE)  # Ensure directory exists
  
  #######################################
  # Set output name defined by the user #
  #######################################
  output_name = as.character(outputName)

  #####################
  # Define color skim #
  #####################
  n <- length(unique(data$Taxa)) # Count the number of total parasites/symbionts
  list = list()
  counter = 1

  for (i in unique(info$Host_sp)){
    species_number = info[info$Host_sp==i,]
    length(species_number$type_species)
    list[[counter]] = length(species_number$type_species)
    counter = counter + 1
  }
  maximum_species = max(do.call(c,list))

  message("\n")    
  message(paste("Number of hosts: ",length(unique(info$Host_sp)),sep=""))   
  message(paste("Number of symbionts: ", n,sep=""))   
  message(paste("Maximum number of symbionts per host: ",maximum_species,sep=""))   
  message("\n")     


  #if(maximum_species < 3){ # If it smaller than 6 is you color code below
  #  couleur = viridis(begin = 0, end = .75, n)
  #}else if((maximum_species > 2) && (maximum_species < 7)){
  #couleur = c(viridis(begin = 0, end = .75, n-1),"#ff8c00")
  # print(couleur)
  #}else{ # If it large use a random color code assigned by the brewer palette
  #  print("Warning: Too many host species the plot is going to be very cluttered and hard to read")
  #  qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
  #col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
  #couleur = col_vector[1:n]
  #}

  if(n < 8){
    couleur = c("#7B3014","#D04A07","#F98C40","black","#5AA5CD","#236CA7","#26456E")
  }else{ # If it large use a random color code assigned by the brewer palette
    warning("Too many host species the plot may be cluttered and hard to read.")
    
    qual_col_pals = RColorBrewer::brewer.pal.info[RColorBrewer::brewer.pal.info$category == 'qual',]
    col_vector = unlist(mapply(RColorBrewer::brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
    couleur = col_vector[1:n]
  }
  
  message("Color scheme defined")   
  message("Start plotting")   
  
  ###################
  # Start plotting  #
  ###################
  # The plotting works by adding one taxa (parasites/symbionts) at a time as a new layer.
  count = 1 # Counter that will by incremented for each taxa. This counter is used to choice the colors and figure out the spatial coordinate of the geom_text for the stability threshold.
  data$Taxa = gsub(pattern = "[_-]", replacement = " ", x=data$Taxa, perl = TRUE)
  data$Host_sp = gsub(pattern = "[_-]", replacement = " ", x=data$Host_sp, perl = TRUE)
  info$type_species = gsub(pattern = "[_-]", replacement = " ", x=info$type_species, perl = TRUE)
  info$Host_sp = gsub(pattern = "[_-]", replacement = " ", x=info$Host_sp, perl = TRUE)

  suppressMessages(for (taxonomy_names in sort(unique(data$Taxa))){ # Counter that will count the number of
    if(count==1){ # If count = 1 make the initial plotting with the ggplot function
      plot_data = data[data$Taxa==taxonomy_names,] #  Loop over taxa
      plot_data = plot_data %>%
        dplyr::group_by(Host_sp,Substract,colonies) %>%
        dplyr::summarise(avg = mean(Prevalence), lci = Rmisc::CI(Prevalence,ci = 0.95)[3], uci = Rmisc::CI(Prevalence,ci = 0.95)[1])
      plot_data$colonies = as.numeric(plot_data$colonies)
      plot_data$Taxa = as.factor(taxonomy_names)

      p = ggplot2::ggplot(plot_data, ggplot2::aes(x=colonies, y=avg,color=Taxa)) + ggplot2::geom_point(size=0.5) + ggplot2::scale_color_manual(values = couleur)
      p = p +  ggplot2::theme_bw() + ggplot2::facet_wrap(Host_sp~., scales = "free",  ncol = 3) +  ggplot2::xlab("Number of samples") + ggplot2::ylab("Prevalence (%)")
      p = p + ggplot2::theme(strip.background =element_rect(fill="wheat1"),strip.text.x = ggplot2::element_text(size = 14,face = "bold.italic"),legend.text=ggplot2::element_text(size=15),legend.title=ggplot2::element_text(size=16),axis.title=ggplot2::element_text(size=15),axis.text = ggplot2::element_text(size = 15))
      p = p + ggplot2::geom_ribbon(data=plot_data,aes(ymin=lci,ymax=uci),fill="grey", color="grey",alpha =0.5)
      p = p + ggplot2::geom_vline(data  = info[info$type_species==taxonomy_names,], aes(xintercept = thres),color = couleur[count] , linetype="dotted")
      p = p + ggplot2::geom_text(y = Inf, aes(x = Inf, label = ifelse(thres < 100, gsub(pattern = "(.+)", " \\1 ", thres), thres)),data = info[info$type_species==taxonomy_names,], color = couleur[count],hjust =1.2, vjust = 2,size = 4.5)
      p = p + ggplot2::geom_text(y = Inf, aes(x = Inf, label = gsub(pattern = "(.+)", " \\1 ", missing)), data = info[info$type_species==taxonomy_names,], color = couleur[count],hjust =1.5, vjust = 2,size = 4.5)
      p = p + ggplot2::geom_text(aes(x = Inf, y = Prevalence, label = round(Prevalence, digits = 1)), data =  info[info$type_species==taxonomy_names,], color = couleur[count],hjust = 1.5, vjust = 1.5,size = 4.5)
      p = p + ggplot2::guides(colour = guide_legend(override.aes = list(size=2)))
      count = count + 1
    }
    else{ #  if count > 1 add the successive layers (taxa) in each loop
      plot_data = data[data$Taxa==taxonomy_names,]
      plot_data = subset(data, data$Taxa==taxonomy_names,drop = FALSE)
      plot_data = plot_data %>%
        dplyr::group_by(Host_sp,Substract,colonies) %>%
        dplyr::summarise(avg = mean(Prevalence), lci = Rmisc::CI(Prevalence,ci = 0.95)[3], uci = Rmisc::CI(Prevalence,ci = 0.95)[1])
      plot_data$colonies = as.numeric(plot_data$colonies)
      plot_data$Taxa = as.factor(taxonomy_names)

      p = p + ggplot2::geom_point(data = plot_data, ggplot2::aes(x=colonies, y=avg, color=Taxa) ,size=0.5) +  ggplot2::scale_color_manual(values = couleur)
      p = p + ggplot2::theme_bw() + ggplot2::facet_wrap(Host_sp~., scales = "free",  ncol = 3)
      p = p + ggplot2::theme(strip.background =element_rect(fill="wheat1")) + ggplot2::theme(legend.position="bottom")
      p = p + ggplot2::theme(strip.text = element_text(face = "bold.italic",size = 14),legend.text=element_text(size=15,face="italic"),legend.title=element_text(size=16),axis.title=element_text(size=15),axis.text = element_text(size = 15))
      p = p + ggplot2::geom_ribbon(data=plot_data,aes(ymin=lci,ymax=uci),fill="grey", color="grey",alpha =0.5)
      p = p + ggplot2::geom_vline(data  = info[info$type_species==taxonomy_names,], aes(xintercept = thres),color = couleur[count] , linetype="dotted")
      p = p + ggplot2::geom_text(y = Inf, aes(x = Inf, label = ifelse((thres < 100) | (is.na(thres)), gsub(pattern = "(.+)", " \\1 ", thres), thres)), data = info[info$type_species==taxonomy_names,], color = couleur[count],hjust =  count * 1.2, vjust = 2,size = 4.5)
      p = p + ggplot2::geom_text(y = Inf, aes(x = Inf, label = gsub(pattern = "(.+)", " \\1 ", missing)), data = info[info$type_species==taxonomy_names,], color = couleur[count],hjust = count * 1.5, vjust = 2,size = 4.5)
      p = p + ggplot2::geom_text(aes(x = Inf, y = Prevalence, label = round(Prevalence, digits = 1)), data =  info[info$type_species==taxonomy_names,], color = couleur[count],hjust = 1, vjust = 1.5,size = 4.5)
      p = p + ggplot2::guides(colour = guide_legend(override.aes = list(size=2)))
      count = count + 1
    }
  })

  ###################
  # PLotting per se #
  ###################
  nb_sp = length(unique(data$Host_sp)) # Check the number of species. Import to know how large the pdf plot has to be.

  if(nb_sp > 1 && nb_sp < 4){ # If there are less than 4 species
    pdf(file.path(outputDir, paste0(output_name, ".pdf")), length(unique(data$Host_sp)) * 5, 6) # this define the width and the length of the pdf if there are less than 4 species to plot
    plot(p)
    dev.off()
  } else if (nb_sp > 4){ # If there are more than 4 species
    pdf(file.path(outputDir, paste0(output_name, ".pdf")), 12, round(length(unique(data$Host_sp)) / 4) * 4) # this define the width and the length of the pdf if there are more than 4 species to plot
    plot(p)
    dev.off()
  } else {  # If there are exactly 4 species
    pdf(file.path(outputDir, paste0(output_name, ".pdf")), 8, 5)   # this define the width and the length of the pdf if there are four species to plot
    plot(p)
    dev.off()
  }

  ###############################
  # Save results as a text file #
  ###############################
  textfile = info[,c(1,4,3,2)] # Create an object textfile coutaning the same information as info but reshape in a more meaningful order.
  colnames(textfile) = c("Host_species","Taxa","Prevalence","thres_stability")
  write.table(file = file.path(outputDir, paste0(output_name, ".txt")), 
              x = textfile, quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)
}