R/summaryPlots.R

In rhondabacher/methylscaper: Visualization of Methylation Data

#' Calculates the percentage of methylated cells/molecules per site
#'
#' @param orderObject An object of class \code{orderObject}
#' @param makePlot Logical, indicates whether to generate the percentage plot
#' @param ... Additional parameters used by the \code{plot} function.
#' @return The percent of molecules or cells methylated (endogenous (yellow)
#'      or accessible) at each site. Output is a list with names "red"
#'      and "yellow". Red represents the endogenous methylation and
#'      yellow represents the accessibility. Within each list object is a
#'      vector of the percent of cells/molecules methylated. The location
#'      of the site is also represent in the form CXX, where XX is the
#'      position of the site within the defined region.
#' @importFrom graphics hist lines plot points
#' @export
#' @examples
#'
#' data(singlemolecule_example)
#'
#' orderObj <- initialOrder(singlemolecule_example, Method = "PCA")
#' methyl_percent_sites(orderObj, makePlot = TRUE)
methyl_percent_sites <- function(orderObject, makePlot = TRUE, ...) {
    dat <- orderObject$toClust

    red_sites <- c(which(dat == 4, arr.ind = TRUE)[, 2], which(dat == 1, arr.ind = TRUE)[, 2])
    red_sites <- sort(unique(red_sites))
    yellow_sites <- c(which(dat == -4, arr.ind = TRUE)[, 2], which(dat == -1, arr.ind = TRUE)[, 2])
    yellow_sites <- sort(unique(yellow_sites))

    c_red <- vapply(red_sites, function(i) {
        sum(dat[, i] == 4) / nrow(dat)
    }, numeric(1))
    c_yellow <- vapply(yellow_sites, function(i) {
        sum(dat[, i] == -4) / nrow(dat)
    }, numeric(1))
    if (makePlot) {
        plot(
            x = red_sites - ncol(dat) / 2, y = c_red,
            col = "brown1", pch = 19, ylim = c(0, 1),
            xlab = "Position in the genomic region", ylab = "% methylated",
            bty = "n", cex.lab = 1.3, ...
        )
        lines(x = red_sites - ncol(dat) / 2, y = c_red, col = "brown1")
        points(x = yellow_sites, y = c_yellow, col = "gold2", pch = 19)
        lines(x = yellow_sites, y = c_yellow, col = "gold2")
        legend("topright", c("Endogenous", "Accessibility"),
            pch = 16,
            col = c("brown1", "gold2"), bty = "n", title = "Methylation Type"
        )
        n_sites <- length(union(red_sites, yellow_sites))
        labs <- union(red_sites, yellow_sites)[seq(1, n_sites, by = n_sites / 12)]
    }
    final <- list(c_red, c_yellow)
    names(final) <- c("meth", "acc")
    return(final)
}

#' Calculate the proportion of methylated bases for each cell/molecule
#'
#' @param orderObject An object of class \code{orderObject}
#' @param type Indicates which data set to compute proportions for.
#'      This should be 'met' or 'hcg' or 'red' for endogenous methylation; 'acc' or
#'      'gch' or 'yellow' for accessibility.
#' @param makePlot Indicates whether to plot a histogram of the proportions
#'  across all cells/molcules.
#' @param ... Additional parameters used by the \code{hist} function.
#'
#' @return The proportion of methylated (endogenous (yellow)
#'      or accessible) bases for each cell/molecule. Output is vector
#'      with length the numbner of cells/molecules and contains a proportion.
#' @importFrom graphics hist
#' @export
#' @examples
#'
#' data(singlemolecule_example)
#'
#' orderObj <- initialOrder(singlemolecule_example, Method = "PCA")
#' methyl_proportion(orderObj, makePlot = TRUE)
methyl_proportion <- function(orderObject, type = "yellow",
    makePlot = TRUE, ...) {
    type <- tolower(type)
    if (type == "gch") type <- "yellow"
    if (type == "accessibility methylation") type <- "yellow"
    if (type == "acc") type <- "yellow"
    color_indicator <- ifelse(type == "yellow", -1, 1)
    Proportion <- apply(orderObject$toClust, 1, function(x) {
        sum(x == color_indicator * 3 | x == color_indicator * 4) / (length(x) / 2)
    })
    if (makePlot) {
        opar <- par(lwd = 4)
        H <- hist(Proportion, plot = FALSE, breaks = 15)
        plot(H,
            xlim = c(0, 1), border = ifelse(type == "yellow", "gold2", "brown1"),
            col = "gray75", cex.lab = 1.3,
            lwd = 2, ...
        )
        par(opar)
    }
    Proportion <- data.frame(methylationProportion = Proportion)
    return(Proportion)
}

#' Calculate the average methylation/accessibility status across all cells/molecules.
#'
#' @param orderObject An object of class \code{orderObject}
#' @param window_length Length of the window to be used to compute a
#'          moving average. Default is 20.
#' @param makePlot Logical, indicates whether to generate a line plot of
#'          average status.
#' @param ... Addition parameters used by the \code{plot} function.
#'
#' @return The proportion of methylated bases for each cell/molecule
#'      within a defined moving window. Output is a list with elements
#'      "meth_avg" and "acc_avg", indicating endogenous
#'      or accessible methylation respectively.
#' @importFrom stats filter
#' @importFrom graphics legend
#' @export
#' @examples
#'
#' data(singlemolecule_example)
#'
#' orderObj <- initialOrder(singlemolecule_example, Method = "PCA")
#' methyl_average_status(orderObj, makePlot = TRUE)
#'
methyl_average_status <- function(orderObject, window_length = 20,
    makePlot = TRUE, ...) {
    colLength <- ncol(orderObject$toClust)
    gch_num <- orderObject$toClust[, seq(1, (colLength / 2))]
    hcg_num <- orderObject$toClust[, seq((colLength / 2 + 1), colLength)]

    acc_sum <- colSums(gch_num <= -3)
    meth_sum <- colSums(hcg_num >= 3)
    acc_denom <- colSums(gch_num != 0)
    meth_denom <- colSums(hcg_num != 0)
    acc_avg <- acc_sum / acc_denom
    meth_avg <- meth_sum / meth_denom

    width <- window_length
    moving_acc_avg <- filter(x = acc_avg, filter = rep(1, width)) / width
    moving_meth_avg <- filter(x = meth_avg, filter = rep(1, width)) / width

    if (makePlot) {
        plot(
            x = seq(1, length(moving_acc_avg)), y = moving_acc_avg,
            type = "l", col = "gold2", lwd = 2, bty = "n",
            xlab = "Position in the genomic region",
            ylab = "Averaged % methylation status",
            cex.lab = 1.3, ylim = c(0, 1), ...
        )
        lines(moving_meth_avg, col = "brown1", lwd = 2)
        legend("topright", c("Endogenous", "Accessibility"),
            lwd = 2,
            col = c("brown1", "gold2"), bty = "n", title = "Methylation Type"
        )
    }
    return(list(meth_avg = moving_meth_avg, acc_avg = moving_acc_avg))
}