Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions R/BSseq_utils.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)) /
Expand Down Expand Up @@ -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)
Expand Down
9 changes: 7 additions & 2 deletions man/getMeth.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -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}.}
Expand All @@ -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
Expand All @@ -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.
Expand Down