R/edmr.R

In Shicheng-Guo/edmr: Empirical differentially methylated regions calculation

#' Empirical differentially methylated regions
#' 
#' Comprehensive DMR analysis based on bimodal normal distribution model and 
#' weighted cost function for regional methylation analysis optimization. 
#' It captures the regional methylation modification by taking the spatial 
#' distribution of CpGs into account for the enrichment DNA methylation 
#' sequencing data so as to optimize the definition of the empirical regions. 
#' Combined with the dependent adjustment for regional p-value combination.
#' 
#' @param myDiff a \code{myDiff} object from \code{methylKit} package. Required.
#' @param step a numeric variable for calculating auto-correlation, default: 100.
#' @param dist distance cutoff to call a gap for DMR, default: "none", which will be automatically determined by the bimodal normal distribution, default: 100.
#' @param DMC.qvalue qvalue cutoff for DMC definition, default: 0.01
#' @param DMC.methdiff methylation difference cutoff for DMC definition, default: 25.
#' @param num.DMCs cutoff of the number DMCs in each region to call DMR, default: 1.
#' @param num.CpGs cutoff of the number of CpGs, default: 3.
#' @param DMR.methdiff cutoff of the DMR mean methylation difference, default=20.
#' @param plot plot the bimodal normal distribution fitting or not, default=FAlSE.
#' @param main the title of the plot, if plot=TRUE. Default=FALSE.
#' @param mode the mode of call DMRs. 1: using all CpGs together. 2: use unidirectional CpGs to call DMRs. default: 1.
#' @param ACF p-value combination test with (TRUE, default) or without (FALSE) dependency adjustment. 
#' @param fuzzypval p-value cutoff for raw regions definition.
#' @return \code{GRanges}
#' @export
#' @examples
#' library(methylKit)
#' library(GenomicRanges)
#' library(mixtools)
#' library(data.table)
#' data(edmr)
#' mydmr=edmr(chr22.myDiff, mode=1, ACF=TRUE)
#' mysigdmr=filter.dmr(mydmr)
#' mydmr2=edmr(chr22.myDiff, mode=2, ACF=TRUE)
#' mydmr3=edmr(chr22.myDiff, mode=1, ACF=FALSE)
#' mydmr4=edmr(chr22.myDiff, mode=2, ACF=FALSE)

edmr=function(myDiff, step=100, dist="none", DMC.qvalue=0.01, DMC.methdiff=25, num.DMCs=1, num.CpGs=3, DMR.methdiff=20, plot=FALSE, main="", mode=1, ACF=TRUE, fuzzypval=1){
  if(mode==1){
    myDMR=eDMR.sub(myDiff, step=step, dist=dist, DMC.qvalue=DMC.qvalue, DMC.methdiff=DMC.methdiff, num.DMCs=num.DMCs, num.CpGs=num.CpGs, 
                   DMR.methdiff=DMR.methdiff, granges=TRUE, plot=plot, main=main, direction="both", ACF=ACF)
  } else if (mode==2){
    hyper.myDMR=eDMR.sub(myDiff, step=step, dist=dist, DMC.qvalue=DMC.qvalue, DMC.methdiff=DMC.methdiff, num.DMCs=num.DMCs, num.CpGs=num.CpGs,
                         DMR.methdiff=DMR.methdiff, granges=TRUE, plot=plot, main=main, direction="hyper", ACF=ACF)
    hypo.myDMR=eDMR.sub(myDiff, step=step, dist=dist, DMC.qvalue=DMC.qvalue, DMC.methdiff=DMC.methdiff, num.DMCs=num.DMCs, num.CpGs=num.CpGs,
                        DMR.methdiff=DMR.methdiff, granges=TRUE, plot=plot, main=main, direction="hypo", ACF=ACF)
    myDMR=c(hyper.myDMR, hypo.myDMR)
  }
  else {
    stop ("mode = 1 or 2")
  }
}