R/DifferentialGeneExpression.R

In sagrudd/cDNA_tutorial: Tutorial For Nanopore Sequence Based Gene Expression Profiling

#' Prepare a volcanoPlot of gene expression data using ggplot2
#'
#' Volcano plot presents a scatter plot of the -log10(p-value) / log2 gene
#' expression from a differential gene expression study
#'
#' @param deSeqObj DESeq2 differential expression object
#' @param lfcThreshold a log2 threshold required for definition of differential expression
#' @param adjPValueThreshold minimal threshold of FDR corrected p-value for differential expression
#'
#' @return None
#'
#' @author Oxford Nanopore Bioinformatics, \email{support@@nanoporetech.com}
#' @references \url{https://en.wikipedia.org/wiki/Volcano_plot_(statistics)}
#' @keywords volcano
#'
#' @export
volcanoPlot <- function(deSeqObj, lfcThreshold, adjPValueThreshold) {
  volcanoSubstr <- results(deSeqObj)
  logUp <- which(volcanoSubstr$log2FoldChange >= lfcThreshold)
  logDown <- which(volcanoSubstr$log2FoldChange <= -lfcThreshold)
  withStat <- which(volcanoSubstr$padj <= adjPValueThreshold)
  colours <- c(noDifference="gray", upRegulated="red", downRegulated="green")
  gene <- rep("noDifference", nrow(volcanoSubstr))
  gene[logUp[logUp %in% withStat]] <- "upRegulated"
  gene[logDown[logDown %in% withStat]] <- "downRegulated"

  plot <- ggplot(data.frame(volcanoSubstr), aes(x=log2FoldChange, y=-log10(padj))) +
    geom_point(size=1.2) +
    geom_hline(yintercept=-log10(adjPValueThreshold), color="orange") +
    geom_vline(xintercept=-lfcThreshold, color="green") +
    geom_vline(xintercept=lfcThreshold, color="red") +
    aes(colour=gene) +
    scale_colour_manual(values=colours) +
    ggtitle("Volcano plot showing distribution of lfc vs adjp for expressed genes")

  print(plot)
}