R/epiMA.R

In TheJacksonLaboratory/epihet: Determining Epigenetic Heterogeneity from Bisulfite Sequencing Data

#' @title Make MA Plot
#' @description
#' Creates an MA plot from the differential
#' heterogeneity data calculated from the diffHet()
#' function. For each loci, graphs the average of both
#' group means on the x-axis and the heterogeneity difference
#' on the y-axis. Graphs coordinates with significant adjusted
#' p-values in red.
#' @param pval.matrix The data frame returned from the
#' diffHet() function that contains means, p-values,
#' adjusted p-values, and heterogeneity difference
#' @param padjust.cutoff The adjusted p-value cutoff
#' to confirm a significant value. (default: 0.05)
#' @param pch The plotting character to be used in the
#' MA plot (default: '.')
#' @param sve A boolean to save the plot (default: FALSE)
#' @param pointsize A numeric value to adjust point size (default:1.5)
#' @return A figure
#' @examples
#' diff.het.matrix<-data.frame(chromosome=c(rep("1",10)),
#' loci=paste("loci",seq_len(10),sep="-"),subtype.mean=c(0.21,0.23,0.37,0.26,
#' 0.29,0.31,0.29,0.13,0.12,0.093),Normal.mean=c(0.01,0.01,0.01,0.02,
#' 0.02,0.01,0.01,0.79,0.73,0.79),het.dif=c(0.20,0.220,0.360,0.240,0.270,
#' 0.300,0.280,-0.660,-0.610,-0.697),p.value=c(3.08e-03,1.43e-02,9.27e-03,
#' 3.45e-02,2.99e-02,3.68e-02, 4.60e-02, 5.65e-10, 9.18e-10,
#' 9.98e-11),p.adjust=c(8.84e-03,2.76e-02, 2.04e-02, 5.01e-02,
#' 4.56e-02, 5.24e-02, 6.08e-02, 3.74e-08, 5.22e-08,
#' 1.06e-08),type=rep("pdr",10))
#'
#' epiMA(pval.matrix = diff.het.matrix, padjust.cutoff = 0.05,
#' pch = ".", sve = TRUE)
#' @export
epiMA <- function(pval.matrix, padjust.cutoff = 0.05, pch = ".", sve = FALSE,pointsize=1.5) {
    if (!is(pval.matrix[,3],"numeric") | !is(pval.matrix[,4],"numeric")){
      stop("Error: The third and forth column should be the mean value of epigenetic heterogeneity in cancer or normal samples")
    }
    sig.values <- pval.matrix[pval.matrix$p.adjust <= padjust.cutoff, ]
    sig.means <- (sig.values[, 3] + sig.values[, 4])/2
    means <- (pval.matrix[, 3] + pval.matrix[, 4])/2
    max.het.dif <- max(abs(pval.matrix$het.dif))
    columns <- colnames(pval.matrix)
    group1 <- strsplit(columns[3], ".mean")[[1]]
    group2 <- strsplit(columns[4], ".mean")[[1]]
    value <- pval.matrix$type[1]
    pdftitle <- paste0(group1, "_", group2, "_", value, "_MAplot.pdf")
    title <- paste0("MA Plot of ", group1, " & ", group2,
        " by ", value, " Values")
    x.label <- paste0("Mean Epigenetic Heterogeneity by ",value)
    if (sve) {
        pdf(pdftitle)
        plot(means, pval.matrix$het.dif, col = "grey",
            xlab = x.label, ylab = "Heterogeneity Difference",
            pch = pch, ylim = c(-max.het.dif, max.het.dif),
            main = title,cex=pointsize,cex.lab=1.5,cex.axis=1.5)
        points(sig.means, sig.values$het.dif, col = "red",pch = pch,cex=pointsize)
        dev.off()
    } else {
        plot(means, pval.matrix$het.dif, col = "grey",
            xlab = x.label, ylab = "Heterogeneity Difference",
            pch = pch, ylim = c(-max.het.dif, max.het.dif),
            main = title,cex=pointsize,cex.lab=1.5,cex.axis=1.5)
        points(sig.means, sig.values$het.dif, col = "red",pch = pch,cex=pointsize)
    }
}