From fe4f7559604c9fe2b874dd64a2059fa68d275e20 Mon Sep 17 00:00:00 2001 From: Philipp Rentzsch Date: Mon, 24 Oct 2016 10:10:54 +0200 Subject: [PATCH 1/2] added remove NA to getMeth to allow returning NA for regions if any CpG's coverage is 0 --- R/BSseq_utils.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/BSseq_utils.R b/R/BSseq_utils.R index 43fe02e..6891850 100644 --- a/R/BSseq_utils.R +++ b/R/BSseq_utils.R @@ -35,7 +35,8 @@ orderBSseq <- function(BSseq, seqOrder = NULL) { getMeth <- function(BSseq, regions = NULL, type = c("smooth", "raw"), - what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95) { + what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95, + na.rm = TRUE) { p.conf <- function(p, n, alpha) { z <- abs(qnorm((1 - alpha)/2, mean = 0, sd = 1)) upper <- (p + z^2/(2*n) + z*sqrt( (p*(1-p) + z^2/(4*n)) / n)) / @@ -96,12 +97,12 @@ getMeth <- function(BSseq, regions = NULL, type = c("smooth", "raw"), if(what == "perRegion") { if(type == "smooth") { out <- lapply(mmsplit, function(xx) { - colMeans(getBSseq(BSseq, "trans")(getBSseq(BSseq, "coef")[xx,,drop = FALSE]), na.rm = TRUE) + colMeans(getBSseq(BSseq, "trans")(getBSseq(BSseq, "coef")[xx,,drop = FALSE]), na.rm = na.rm) }) } if(type == "raw") { out <- lapply(mmsplit, function(xx) { - colMeans(getBSseq(BSseq, "M")[xx,,drop = FALSE] / getBSseq(BSseq, "Cov")[xx,,drop = FALSE], na.rm = TRUE) + colMeans(getBSseq(BSseq, "M")[xx,,drop = FALSE] / getBSseq(BSseq, "Cov")[xx,,drop = FALSE], na.rm = na.rm) }) } out <- do.call(rbind, out) From 0e4c82a670e2b08dac7f0bbd47d8b46996ce5d6c Mon Sep 17 00:00:00 2001 From: Philipp Rentzsch Date: Mon, 24 Oct 2016 10:18:52 +0200 Subject: [PATCH 2/2] Added manual for na.rm parameter in getMeth --- man/getMeth.Rd | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/man/getMeth.Rd b/man/getMeth.Rd index 8769873..23c00d7 100644 --- a/man/getMeth.Rd +++ b/man/getMeth.Rd @@ -8,7 +8,8 @@ } \usage{ getMeth(BSseq, regions = NULL, type = c("smooth", "raw"), - what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95) + what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95, + na.rm = TRUE) } \arguments{ \item{BSseq}{An object of class \code{BSseq}.} @@ -20,6 +21,7 @@ getMeth(BSseq, regions = NULL, type = c("smooth", "raw"), methylation estimates (see below). This is only supported if \code{what} is equal to \code{perBase}.} \item{alpha}{alpha value for the confidence interval.} + \item{na.rm}{whether to ignore NA methylation values for region averages} } \note{ A \code{BSseq} object needs to be smoothed by the function @@ -44,7 +46,10 @@ getMeth(BSseq, regions = NULL, type = c("smooth", "raw"), corresponding to a genomic region (and each row of the matrix being a loci inside the region). If \code{what = "perRegion"} the function returns a matrix, with each row corresponding to a region and - containing the average methylation level in that region. + containing the average methylation level in that region. If + \code{na.rm = FALSE}, the average methylation level of regions is given + as NA if any CpG in the region has a coverage of 0 and therefore + methylation level of NA. } \references{ A Agresti and B Coull.