R/plotMIDASComp.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
#############################################################

gatherMIDASlist <- function(MIDASlist) {
  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)
  }
  
  return(data)
}

expandLimits <- function(limits, factor) {
  expand_by <- (limits[2] - limits[1]) * (1 + (factor/2))
  return(c(limits[1]-factor, limits[2]+factor))
  
}

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
}

plotMIDASComp <- function(expMIDASlist, simMIDASlist, errMat = NA, dataWidth = NA) {
  # First we must do some processing to put the MIDASlist data into a format that can be plotted (key-value-ish)
  expData <- gatherMIDASlist(expMIDASlist)
  simData <- gatherMIDASlist(simMIDASlist)
  allData <- rbind(expData, simData)
  
  # Do something similar with the error matrix
  errMIDASlist <- expMIDASlist
  errMIDASlist$valueSignals <- array(errMat, dim = c(nrow(expMIDASlist$valueCues), length(expMIDASlist$namesSignals),1))
  errMIDASlist$timeSignals <- c(1)
  errData <- gatherMIDASlist(errMIDASlist)
  
  # 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(expMIDASlist$valueCues)
  colnames(cues) <- expMIDASlist$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 expMIDASlist$namesCues) { # Set cues to their types
    cuesGather$value[cuesGather$cue == cue & cuesGather$value == 1] <- expMIDASlist$treatmentDefs$type[expMIDASlist$treatmentDefs$name == cue]
  }
  cuesGather$value <- factor(cuesGather$value)
  
  # Check which panels will have no experimental data in them and make a new df with that info
  allPanels <- unique(expData[,c("cueNum","signal")]) # Get all the panels
  dataPanels <- unique(expData[!is.na(expData$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
  expData <- expData[!(is.na(expData$value)),]
  simData <- simData[!(is.na(simData$value)),]

  # Data plot
  xlimits <- c(min(allData$time, na.rm = T), max(allData$time, na.rm = T))
  ylimits <- c(min(allData$value, na.rm = T), max(allData$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=errData, ggplot2::aes(x=mean(xlimits), y=mean(ylimits), height=panelHeight,
                               width = panelWidth, fill=value), alpha = 0.6) +
    ggplot2::scale_fill_gradient(low="green", high="red", limits=c(0,1), position="left") +
    ggplot2::geom_line(data=simData, ggplot2::aes(x=time, y=value), linetype = "twodash", color = "gray50") +
    ggplot2::geom_point(data=simData, ggplot2::aes(x=time, y=value), color = "gray50") +
    ggplot2::geom_line(data=expData, ggplot2::aes(x=time, y=value)) +
    ggplot2::geom_point(data=expData, ggplot2::aes(x=time, y=value)) +
    # ggplot2::geom_errorbar(data=data, ggplot2::aes(ymin=value-sqrt(variance), ymax=value+sqrt(variance)), width=10) + # 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::xlab("Time (min)") + ggplot2::ylab("Signal") + ggplot2::labs(fill="MSE") +
    ggplot2::facet_grid(cols=vars(signal), rows=vars(cueNum)) + # This splits the graph into times and cues
    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"), legend.position = "left")

  # 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
    expData <- rbind(expData, emptyPanels)
    daataPlot <- dataPlot + ggplot2::geom_tile(data=emptyPanels, ggplot2::aes(x=mean(xlimits), y=mean(ylimits), height=panelHeight,
                                                            width = panelWidth), fill = "grey", alpha = 0.6)
  }


  # 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(expMIDASlist$namesSignals), length(expMIDASlist$namesCues)/4)
  } else {
    widths = c(dataWidth, 1)
  }
  
  # Put them together and plot
  return(suppressWarnings(suppressMessages(egg::ggarrange(dataPlot, cuesPlot, widths = widths))))
}