R/plotMIDAS.R

In palmerito0/kboolnet: kboolnet: A Toolkit for the Verification, Validation, and Visualization of rxncon Models

#############################################################
# plotMIDASlist.R
# Adrian C
#
# Plots the output of readMIDASExcel. Loosely based on 
# plotCNOlist from the CellNOptR package.
#
# Arguments:
#   - MIDASlist: MIDAS list generated by readMIDASExcel
#   - dataWidth: Relative width of the data compared to the cues, optional
#
# Depndencies: tidyr, dplyr, ggplot2, egg
#############################################################

no_zero <- function(x) {
  y <- sprintf('%.1f',x)
  y[x > 0 & x < 1] <- sprintf('.%s',x[x > 0 & x < 1]*10)
  y[x == 0] <- '0'
  y[x > -1 & x < 0] <- sprintf('-.%s',x[x > -1 & x < 0]*-10)
  y
}

plotMIDAS <- function(MIDASlist, dataWidth = NA) {
  # First we must do some processing to put the MIDASlist data into a format that can be plotted (key-value-ish)
  data <- data.frame(cueNum=numeric(), time=numeric(), signal=character(), value=numeric(), variance=numeric())
  for (i in 1:length(MIDASlist$timeSignals)) {
    # Get the values an gather them
    timeValues <- as.data.frame(MIDASlist$valueSignals[,,i]) # Get data values
    timeValues <- cbind(timeValues, 1:nrow(timeValues)) # Add column to keep track of the cue combination
    colnames(timeValues) <- c(MIDASlist$namesSignals, "cueNum") # Set column names
    timeValues <- tidyr::gather(timeValues, "signal", "value", -cueNum)
    
    # Do the same as above for variances
    timeVariances <- as.data.frame(MIDASlist$valueVariances[,,i]) # Get data variances
    timeVariances <- cbind(timeVariances, 1:nrow(timeVariances)) # Add column to keep track of the cue combination
    colnames(timeVariances) <- c(MIDASlist$namesSignals, "cueNum") # Set column names
    timeVariances <- tidyr::gather(timeVariances, "signal", "variance", -cueNum)
    
    # Merge timeVariances and timeValues
    merged <- merge(timeVariances, timeValues)
    
    # Set the time column
    merged$time <- MIDASlist$timeSignals[i]
    
    # Append data to the original df
    data <- rbind(data, merged)
  }
  
  # Remove NA values - this causes weird behavior, don't do it
  # data <- data[!is.na(data$value),]
  
  # Also process the cues data
  cues <- as.data.frame(MIDASlist$valueCues)
  colnames(cues) <- MIDASlist$namesCues
  cues <- cues[,colSums(cues) != 0]
  cues$cueNum <- 1:nrow(cues)
  cuesGather <- tidyr::gather(cues, "cue", "value", -cueNum)
  cuesGather$cue <- factor(cuesGather$cue, levels=MIDASlist$treatmentDefs$name)
  cuesGather$value[cuesGather$value == 0] <- "None" # Set all zeroes to None
  for (cue in MIDASlist$namesCues) { # Set cues to their types
    cuesGather$value[cuesGather$cue == cue & cuesGather$value == 1] <- MIDASlist$treatmentDefs$type[MIDASlist$treatmentDefs$name == cue]
  }
  cuesGather$value <- factor(cuesGather$value)
  
  # Check which panels will have no data in them and make a new df with that info
  allPanels <- unique(data[,c("cueNum","signal")]) # Get all the panels
  dataPanels <- unique(data[!is.na(data$value),][,c("cueNum","signal")]) # Get only panels with data
  emptyPanels <- rbind(allPanels, dataPanels)
  emptyPanels <- emptyPanels[!(duplicated(emptyPanels) | duplicated(emptyPanels, fromLast = T)),] # Remove duplicates, leaving only panels without data
  
  # Remove empty values from data
  data <- data[!(is.na(data$value)),]
  
  # If there are empty panels, add empty values to them so they will still be drawn
  if (nrow(emptyPanels) > 0) {
    emptyPanels[,c("variance", "value", "time")] <- NA
    data <- rbind(data, emptyPanels)
    
    # Data plot
    xlimits <- c(min(data$time, na.rm = T), max(data$time, na.rm = T))
    ylimits <- c(min(data$value, na.rm = T), max(data$value, na.rm =T))
    panelHeight <- (ylimits[2] - ylimits[1]) * 1.3
    panelWidth <- (xlimits[2] - xlimits[1]) * 1.3
    dataPlot <- ggplot2::ggplot() + 
      ggplot2::geom_tile(data=emptyPanels, ggplot2::aes(x=mean(xlimits), y=mean(ylimits), height=panelHeight, width = panelWidth), fill = "grey") +
      ggplot2::geom_line(data=data, ggplot2::aes(x=time, y=value)) +
      ggplot2::geom_point(data=data, ggplot2::aes(x=time, y=value)) +
      ggplot2::geom_errorbar(data=data, ggplot2::aes(x=time, y=value, ymin=value-sqrt(variance), ymax=value+sqrt(variance)), width=1.5) + # Error bars
      ggplot2::coord_cartesian(xlim = xlimits, ylim = ylimits) + # Set graph limits
      ggplot2::scale_y_continuous(breaks=c(ylimits[1], mean(ylimits), ylimits[2]), labels = no_zero) +
      ggplot2::scale_x_continuous(breaks=c(xlimits[1], mean(xlimits), xlimits[2]), labels = no_zero) +
      ggplot2::facet_grid(cols=vars(signal), rows=vars(cueNum)) + # This splits the graph into times and cues
      ggplot2::xlab("Time (min)") + ggplot2::ylab("Signal") + ggplot2::labs(fill="MSE") + # Add labels
      ggplot2::theme(strip.text.y = ggplot2::element_blank(), # Removes the facet labels on the y axis
            panel.border = ggplot2::element_rect(size=1, fill=NA, colour="black"), # Add borders
            strip.background = ggplot2::element_blank(), strip.text.x = ggplot2::element_text(angle=90, hjust=0),
            plot.margin = ggplot2::unit(c(5.5,2.75,5.5,5.5), "pt"))
  } else {
    # Data plot
    xlimits <- c(min(data$time, na.rm = T), max(data$time, na.rm = T))
    ylimits <- c(min(data$value, na.rm = T), max(data$value, na.rm =T))
    panelHeight <- (ylimits[2] - ylimits[1]) * 1.3
    panelWidth <- (xlimits[2] - xlimits[1]) * 1.3
    dataPlot <- ggplot2::ggplot() + 
      ggplot2::geom_line(data=data, ggplot2::aes(x=time, y=value)) +
      ggplot2::geom_point(data=data, ggplot2::aes(x=time, y=value)) +
      ggplot2::geom_errorbar(data=data, ggplot2::aes(x=time, y=value, ymin=value-sqrt(variance), ymax=value+sqrt(variance)), width=1.5) + # Error bars
      ggplot2::coord_cartesian(xlim = xlimits, ylim = ylimits) + # Set graph limits
      ggplot2::scale_y_continuous(breaks=c(ylimits[1], mean(ylimits), ylimits[2]), labels = no_zero) +
      ggplot2::scale_x_continuous(breaks=c(xlimits[1], mean(xlimits), xlimits[2]), labels = no_zero) +
      ggplot2::facet_grid(cols=vars(signal), rows=vars(cueNum)) + # This splits the graph into times and cues
      ggplot2::xlab("Time (min)") + ggplot2::ylab("Signal") + ggplot2::labs(fill="MSE") + # Add labels
      ggplot2::theme(strip.text.y = ggplot2::element_blank(), # Removes the facet labels on the y axis
            strip.background = ggplot2::element_blank(), strip.text.x = ggplot2::element_text(angle=90, hjust=0),
            panel.border = ggplot2::element_rect(size=1, fill=NA, colour="black"), # Add borders
            plot.margin = ggplot2::unit(c(5.5,2.75,5.5,5.5), "pt"))
  }

    
  # Cues plot
  cueColors <- c("None"="white", "Inhibitor"="red2", "Stimulus"="green3", "KO"="darkblue", "Mutant"="magenta4")
  cuesPlot <- ggplot2::ggplot(cuesGather, ggplot2::aes(x=.5, y=.5)) +
    ggplot2::facet_grid(cols=vars(cue), rows=vars(cueNum)) +
    ggplot2::geom_tile(ggplot2::aes(fill=value), colour="black", size=1) + # Create the tiles (colour and size affect borders)
    ggplot2::scale_fill_manual(values=cueColors, drop=FALSE) + # Color the tiles appropriately
    ggplot2::scale_x_continuous(limits=c(0,1), expand=c(0,0)) + ggplot2::scale_y_reverse(limits=c(1,0), expand=c(0,0)) + # Reverse the ordering and set proper scales
    ggplot2::theme(axis.title = ggplot2::element_blank(), axis.ticks = ggplot2::element_blank(), axis.text = ggplot2::element_blank(), # Remove axis labels
          legend.position = "none", panel.background = ggplot2::element_blank(), strip.text.y = ggplot2::element_blank(), # Remove legend, background, and y labels
          strip.background = ggplot2::element_blank(), strip.text.x = ggplot2::element_text(angle=90, hjust=0), # Turn cue labels 90 degrees
          plot.margin = ggplot2::unit(c(0,5.5,0,0), "pt"), plot.background = ggplot2::element_blank()) # Remove margins, these are set by data plot
  
  # If width wasn't provided, make a best guess for it
  if (is.na(dataWidth)) {
    widths = c(length(MIDASlist$namesSignals), length(MIDASlist$namesCues)/4)
  } else {
    widths = c(dataWidth, 1)
  }
  
  # Put them together and plot
  return(suppressWarnings(suppressMessages(egg::ggarrange(dataPlot, cuesPlot, widths = widths))))
}