R/MethCP.R

In boyinggong/MethCP: Differential methylation anlsysis for bisulfite sequencing data

setClassUnion("characterORnumeric", c("character", "numeric"))


#' @title Class \code{MethCP}
#'
#' @description
#' A class for performing DMR (differentially methylated region) detection
#' analysis using method \code{MethCP} on whole genome bisulfite sequencing
#' data.
#'
#' @slot test a character string of the name of the per-cytosine
#' statistcis.
#' @slot group1 a character vector containing the sample names of
#' the treatment group.
#' @slot group2 a character vector containing the sample names of
#' the control group.
#' @slot stat a \code{GRangesList} object containing the results of
#' per-cytosine tests.
#' @slot segmentation a \code{GRanges} object containing the segments
#' and their infomation such as region-based statistics, coverages, etc.
#'
#' @importFrom stats lm na.omit cov qnorm
#' @importFrom bsseq getCoverage
#' @importFrom methods is
#'
#' @name MethCP-class
#' @rdname MethCP-class
#' @exportClass MethCP
MethCP <- setClass(
    "MethCP",
    representation(
        test = "character",
        stat = "GRangesList",
        group1 = "characterORnumeric",
        group2 = "characterORnumeric",
        segmentation = "GRanges"),
    prototype(
        test = NA_character_,
        group1 = NA,
        group2 = NA,
        stat = GenomicRanges::GRangesList(list()),
        segmentation = GenomicRanges::GRanges()))

setGeneric("test", function(x) standardGeneric("test"))
setMethod("test", "MethCP", function(x) x@test)

setGeneric("group1", function(x) standardGeneric("group1"))
setMethod("group1", "MethCP", function(x) x@group1)

setGeneric("group2", function(x) standardGeneric("group2"))
setMethod("group2", "MethCP", function(x) x@group2)

setGeneric("stat", function(x) standardGeneric("stat"))
setMethod("stat", "MethCP", function(x) x@stat)

setGeneric("stat<-", function(x, value) standardGeneric("stat<-"))
setMethod("stat<-", "MethCP", function(x, value) {
    x@stat <- value
    x
})

setGeneric("segmentation", function(x) standardGeneric("segmentation"))
setMethod("segmentation", "MethCP", function(x) x@segmentation)

setGeneric("segmentation<-", function(x, value) 
    standardGeneric("segmentation<-"))
setMethod("segmentation<-", "MethCP", function(x, value) {
    x@segmentation <- value
    x
})

#' @title The constructor function for \code{MethCP} objects.
#'
#' @description
#' The constructor function for \code{MethCP} objects.
#'
#' @usage
#' MethCP(
#'     test = NA_character_, group1 = NA, group2 = NA, chr, pos,
#'     pvals, effect.size)
#'
#' @details
#' If not specified by function \code{calcLociStatTimeCourse},
#' \code{calcLociStat}, the parameter \code{test} can be set to any
#' user-specified string indicating the name of the test performed.
#'
#' In the cases where the goal is not to compare between treatment and control
#' groups, parameter \code{group1} and \code{group2} can be set to \code{NA}.
#'
#' If generated by \code{calcLociStat}, parameter \code{stat} will be a
#' \code{GRangesList} object where each element in the list contains statistics
#' for each of the chromosome in the dataset.
#'
#' @param test a character string of the name of the per-cytosine
#' statistcis.
#' @param group1 a character vector containing the sample names of
#' the treatment group.
#' @param group2 a character vector containing the sample names of
#' the control group.
#' @param chr a character vector containing the cytosine chromosome
#' infomation.
#' @param pos a numeric vector containing the cytosine positions.
#' @param pvals a numeric vector containing the \code{p}-values for each
#' cytosine.
#' @param effect.size a numeric vector containing the effect sizes
#' for each cytosine.
#'
#' @return A \code{MethCP} objects.
#'
#' @examples
#' obj <- MethCP(
#'     test = "myTest",
#'     group1 = paste0("Treatment", 1:3),
#'     group2 = paste0("Control", 1:3),
#'     chr = rep("Chr01", 5),
#'     pvals = c(0, 0.1, 0.9, NA, 0.02),
#'     pos = c(2, 5, 9, 10, 18),
#'     effect.size = c(1, -1, NA, 9, Inf))
#' # MethCP will omit the NAs and infinite values.
#' obj
#'
#' @importFrom IRanges IRanges
#' @importFrom methods new
#' @importFrom stats qnorm na.omit
#'
#' @name MethCP

#' @rdname MethCP-class
#' @export
MethCP <- function(
        test = NA_character_,
        group1 = NA,
        group2 = NA,
        chr, pos, pvals, effect.size)
    {
        if (!(class(chr) %in% c("integer", "character"))){
            stop("chr must be integter or character.")
        }
        if (!(class(pos) %in% c("integer", "numeric"))){
            stop("pos must be integter or numeric.")
        }
        if (!(class(pvals) %in% c("numeric"))){
            stop("pvals must be numeric.")
        }
        if (!all(na.omit(pvals <= 1 & pvals >= 0))){
            stop("pvals must be between 0 and 1.")
        }
        if (!(class(effect.size) %in% c("integer", "numeric"))){
            stop("effect.size must be integter or numeric")
        }
        npos = length(chr)
        if ((length(pos) != npos) | (length(pvals) != npos) |
            (length(pvals) != npos)){
            stop(paste(
                "lengths of chr, pos,",
                "pvals and effect.size do not match."))
        }
        filter = (is.na(pvals)) | (is.na(effect.size)) |
            (is.infinite(pvals)) | (is.infinite(pvals))
        if (any(filter)){
            message(paste0(
                "Filtering out NAs and infinite values. \n",
                "Cytosine counts before filtering: ", npos,
                ". \nCytosine counts after filtering: ",
                sum(!filter), "."))
            chr = chr[!filter]
            pos = pos[!filter]
            pvals = pvals[!filter]
            effect.size = effect.size[!filter]
        }
        norm_stats <- qnorm(1-pvals/2)*sign(effect.size)
        stat_res <- GenomicRanges::GRangesList(lapply(unique(chr), function(o){
            GRanges(
                seqnames = chr[chr == o],
                ranges = IRanges(start = pos[chr == o], width = 1),
                stat = norm_stats[chr == o],
                pval = pvals[chr == o])
        }))
        methcpObj <- new(
            "MethCP",
            group1 = "notApplicable",
            group2 = "notApplicable",
            test = test,
            stat = stat_res)
    }

#' @title The show method
#'
#' @description
#' Print \code{MethCP} object information.
#'
#' @param object a \code{MethCP} object.
#'
#' @return No value will be returned.
#'
#' @export
setMethod("show", signature("MethCP"), function(object){
    cat(paste0(
        "MethCP object with ", length(stat(object)), " chromosomes, ",
        sum(vapply(stat(object), length, FUN.VALUE = numeric(1))),
        " methylation loci\n"))
    cat(paste0("test: ", test(object), "\n"))
    cat(paste0(
        "group1: ", paste(group1(object), collapse = " "),
        "\ngroup2: ", paste(group2(object), collapse = " ")))
    if (length(segmentation(object)) == 0){
        cat("\nhas not been segmented")
    } else {
        cat("\nhas been segmented")
    }
})