R/utils.R
In kaiqiong/SOMNiBUS: Smooth modeling of bisulfite sequencing
Defines functions
splitDataByGene
splitDataByChromatin
splitDataByBed
splitDataByGRanges
splitDataByDensity
splitDataByRegion
hessianComp
Documented in
splitDataByBed
splitDataByChromatin
splitDataByDensity
splitDataByGene
splitDataByGRanges
splitDataByRegion
#' @title hessianComp computes the Hessian matrix
#'
#' @description computes the Hessian matrix for an EM estimate
#' @param w_ij the diagnal values of the weight matrix
#' @param new.par estimate of alpha
#' @param new.lambda estimate of lambda
#' @param X a vector of read depths
#' @param Y a vector of methylated counts
#' @param my.design.matrix design matrix from the final fit
#' @param gam_smoothMat the smooth matrix from the final gam fit
#' @param Z covariate matrix
#' @param pred.pi predicted methylation probability from the final fit
#' @param p0 the probability of observing a methylated read when
#' the underlying true status is unmethylated. \code{p0} is the rate of false
#' methylation calls, i.e. false positive rate.
#' @param p1 the probability of observing a methylated read when
#' the underlying true status is methylated. \code{1-p1} is the rate of false
#' non-methylation calls, i.e. false negative rate.
#' @param disp_est estimated dispersion parameter
#' @param RanEff whether a subject-level Random effect is added or not
#' @param N number of unique samples in the provided dataset
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return a Hessian matrix for all alphas (including the alphas for RE)
#' @references Zhao, Kaiqiong, et al. 'A novel statistical method for
#' modeling covariate effects in bisulfite sequencing derived measures
#' of DNA methylation.'Biometrics (2020).
#' @author Kaiqiong Zhao
#' @importFrom Matrix bdiag
#'
#' @noRd
hessianComp <- function(w_ij, new.par, new.lambda, X, Y, my.design.matrix,
gam_smoothMat, Z, pred.pi, p0, p1, disp_est,
RanEff, N, verbose = TRUE) {
t0 <- Sys.time()
if(verbose) Message("Compute the Hessian matrix for an EM estimate")
## Q1: the second partial derivative w.r.t alpha^2 Q2: the second
## derivative w.r.t alpha & alpha_star
res <- outer(seq_len(length(new.par)),
seq_len(length(new.par)),
Vectorize(function(l,m) {
sum(-X * w_ij * my.design.matrix[, m] * my.design.matrix[, l])
}))
smoth.mat <- lapply(as.list(seq_len(ncol(Z) + 1)), function(i) {
gam_smoothMat[[i]]$S[[1]] * new.lambda[i]
}) ## extract the penalty matrix
## assume the lambda for the constant of the intercept is 0 -- no
## penalization
smoth.mat[[length(smoth.mat) + 1]] <- 0
if (RanEff) {
## !!!! Otherwise, we get very wide CI
smoth.mat[[length(smoth.mat) + 1]] <- diag(N) * new.lambda[ncol(Z) + 2]
span.penal.matrix <- as.matrix(Matrix::bdiag(smoth.mat[c(length(smoth.mat) - 1,
(seq_len((length(smoth.mat) - 2))),
length(smoth.mat))]))
} else {
span.penal.matrix <- as.matrix(Matrix::bdiag(smoth.mat[c(length(smoth.mat),
(seq_len((length(smoth.mat) - 1))))]))
}
Q1_with_lambda <- res - span.penal.matrix/disp_est
Q1_no_lambda <- res
Q2 <- outer(seq_len(length(new.par)), seq_len(length(new.par)),
Vectorize(function(l, m) {
term1 <- Y * p1 * p0/(p1 * pred.pi + p0 * (1 - pred.pi))^2 + (X - Y) *
(1 - p1) * (1 - p0)/((1 - p1) * pred.pi + (1 - p0) * (1 - pred.pi))^2
sum(term1 * w_ij * my.design.matrix[, m] * my.design.matrix[, l])
}))
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(Q1_with_lambda + Q2)
}
#' @title Split methylation data into regions based on the spacing of CpGs
#' @description This function splits the methylation data into regions
#' based on the spacing of CpGs.
#'
#' @param dat a data frame with rows as individual CpGs appearing in all the
#' samples. The first 4 columns should contain the information of `Meth_Counts`
#' (methylated counts), `Total_Counts` (read depths), `Position` (Genomic
#' position for the CpG site) and `ID` (sample ID). The covariate information,
#' such as disease status or cell type composition, are listed in column 5 and
#' onwards.
#' @param gap positive integer defining the gap width
#' beyond which we consider that two regions are independent.
#' Odd and decimal values will be rounded to the next even numbers
#' (e.g. 8.2 and 8.7 become gaps of 8 and 10 respectively).
#' The default value is \code{1e+6} (1Mb).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByRegion( dat=RAdat.f, gap = 1e6, min.cpgs = 5,
#' verbose = FALSE)
#' @importFrom IRanges IRanges reduce gaps start end
#' @export
splitDataByRegion <- function(dat, gap = 1e6, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste("Partitioning the regions based on",
"the spacing of CpGs")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- 1e+6
} else {
if(!is.numeric(gap) || gap < 1){
gap <- 1e+6
w_msg <- paste("'gap' should be a positive integer.",
"We will use the default value 1e+6 (1Mb).")
Warning(w_msg)
} else {
ori_gap <- gap
gap <- round(gap)
if(gap %% 2 != 0) {
gap <- gap + 1
}
}
}
# convert dataframe in IRanges
dat_gr <- IRanges::IRanges(start = dat$Position, end = dat$Position)
# reduce first orders the ranges in x from left to right,
# then merges the overlapping or adjacent ones.
dat_gr <- IRanges::reduce(dat_gr)
# identify gap between regions
g = IRanges::gaps(dat_gr)
# extend each region with half of the maximal distance d
# and attempt reduction to identify overlapping regions
hemi_d <- ceiling(gap/2)
regions_gr <- IRanges::reduce(dat_gr + hemi_d) - hemi_d
# save the identified regions
listOfRegions <- list()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else{
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- sapply(X = listOfRegions,
FUN = function(region)
paste("region",
min(region$Position),
max(region$Position), sep = "_")
)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
# return the list of independent regions
return(listOfRegions)
}
#' @title Split methylation data into regions based on the density of CpGs
#' @description This function splits the methylation data into regions
#' based on the density of CpGs.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status or cell type composition,
#' are listed in column 5 and onwards.
#' @param window.size this positive integer defines the size of the
#' sliding window in bp. Decimal values will be rounded to the nearest integer.
#' The value should be greater than 10. The default value is \code{100} (100 bp)
#' @param by positive integer defines by how many base pairs the
#' window moves at each increment. Decimal values will be rounded to the
#' nearest integer. The default value is \code{1} (1 bp).
#' @param min.density positive integer defines the minimum density
#' threshold for each window. Decimal values will be rounded to the
#' nearest integer. The default value is \code{5} (5 CpGs/\code{window.size}).
#' @param gap positive integer defining the gap width
#' beyond which we consider that two regions are independent.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{10} (10bp).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information. The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByDensity(dat = RAdat.f, window.size = 100, by = 1,
#' min.density = 5, gap = 10, min.cpgs = 50, verbose = FALSE)
#' @importFrom IRanges IRanges reduce start end shift
#'
#' @export
splitDataByDensity <- function(dat, window.size = 100, by = 1, min.density = 5,
gap = 10, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"the density of CpGs")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(window.size)){
window.size <- 100
} else {
if(!is.numeric(window.size) || window.size < 10) {
window.size <- 100
w_msg <- paste("'window.size' should be a positive integer",
"greater than 10. We will use the default value 100.")
Warning(w_msg)
} else {
window.size <- round(window.size)
}
}
if(is.null(by)){
by <- 1
} else {
if(!is.numeric(by) || by < 1) {
by <- 1
w_msg <- paste("'by' should be a positive integer.",
"We will use the default value 1.")
Warning(w_msg)
} else {
by <- round(by)
}
}
if(is.null(min.density)){
min.density <- 5
} else {
if(!is.numeric(min.density) || min.density < 1) {
min.density <- 5
w_msg <- paste("'min.density' should be a positive integer.",
"We will use the default value 5.")
Warning(w_msg)
} else {
min.density <- round(min.density)
}
}
if(is.null(gap)){
gap <- 10
} else {
if(!is.numeric(gap) || gap < 1) {
gap <- 10
w_msg <- paste("'gap' should be a positive integer.",
"We will use the default value 10.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
if(gap >= window.size){
gap <- round(window.size/2)
w_msg <- paste0("'gap' should be lesser than the 'window.size' value. ",
"We will use half of window size (",gap,").")
Warning(w_msg)
}
# convert dataframe in IRanges
dat_gr <- unique(IRanges::IRanges(start = dat$Position, end = dat$Position))
# save the end index for the while loop
endOfRegion <- max(IRanges::end(dat_gr))
# create the starting sliding window
sliding_window <- IRanges::IRanges(start = min(IRanges::start(dat_gr)),
end = (min(IRanges::start(dat_gr)) +
window.size - 1))
cpgs <- NULL
# calculate the density in each sliding window
while(IRanges::start(sliding_window) <= endOfRegion){
# overlap region to window
hits <- IRanges::findOverlaps(query = dat_gr, subject = sliding_window)
# save density
if(is.null(x = cpgs)){
cpgs <- data.frame(start = IRanges::start(sliding_window),
end = IRanges::end(sliding_window),
density = length(hits))
} else {
cpgs <- rbind(cpgs,
data.frame(start = IRanges::start(sliding_window),
end = IRanges::end(sliding_window),
density = length(hits)))
}
# shift window
sliding_window <- IRanges::shift(x = sliding_window, shift = by)
}
# filter on density
cpgs <- cpgs[cpgs$density >= min.density,]
# transform the windows in IRanges to create regions
regions_gr <- IRanges::IRanges(start = cpgs$start, end = cpgs$end)
# each window is resized to identify independent regions
regions_gr <- IRanges::resize(x = regions_gr, width = (gap + 1))
regions_gr <- IRanges::reduce(regions_gr)
# save the identified regions
listOfRegions <- list()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else{
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- sapply(X = listOfRegions,
FUN = function(region)
paste("region",
min(region$Position),
max(region$Position), sep = "_")
)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genomic annotations
#' @description This function splits the methylation data into regions
#' based on the genomic annotations provided under the form of a GenomicRanges
#' object.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param annots GenomicRanges object containing the annotations
#' @param gap integer defining the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1} (meaning strict overlaping).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' annot <- GenomicRanges::GRanges(seqnames = "chr1", IRanges::IRanges(
#' start = c(102711720,102711844,102712006,102712503,102712702),
#' end = c(102711757,102711909,102712195,102712637,102712712)
#' ))
#' results <- splitDataByGRanges(dat = RAdat.f,
#' chr = rep(x = "chr1", times = nrow(RAdat.f)),
#' annots = annot, gap = -1, min.cpgs = 5)
#' @importFrom IRanges IRanges reduce gaps start end shift
#' @importFrom GenomicRanges GRanges findOverlaps
#' @importFrom rtracklayer import
#' @importFrom GenomeInfoDb seqlevelsStyle seqinfo
#' @importFrom S4Vectors subjectHits queryHits
#'
#' @export
splitDataByGRanges <- function(dat, chr, annots, gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
if(!inherits(annots, "GenomicRanges"))
Error("'annots' should inherit the GRanges class.")
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
# convert dataframe in GRanges
dat$Chr <- chr
dat_gr <- GenomicRanges::GRanges(seqnames = dat$Chr,
ranges = IRanges::IRanges(start = dat$Position,
end = dat$Position))
dat_gr <- unique(dat_gr)
# error handling
dat_style <- tryCatch({
# check if style has a compatible entry for the species supported by Seqname
GenomeInfoDb::seqlevelsStyle(dat_gr)
}, error=function(e) e)
if(inherits(dat_style, "error")){
stop_msg <- paste0("Incompatible chromosome name. ",
"Please see genomeStyles() for supported styles.")
Error(stop_msg)
}
# error handling
hits <- tryCatch({
# adapt the style (chromosome nomenclature) of the annotation
# to the style of the input data if needed
if(length(intersect(x = GenomeInfoDb::seqlevelsStyle(annots),
y = dat_style)) == 0)
suppressWarnings(GenomeInfoDb::seqlevelsStyle(annots) <- dat_style[1])
# find overlaps
GenomicRanges::findOverlaps(query = dat_gr, subject = annots, maxgap = gap)
}, error=function(e) e)
if(inherits(hits, "error")){
stop_msg <- paste0("Incompatible formats between ",
"the input data and the annotations.")
Error(stop_msg)
}
cpgs <- NULL
if(length(hits) > 0){
for(anno in unique(S4Vectors::subjectHits(hits))){
hit <- hits[S4Vectors::subjectHits(hits) == anno]
# save cpgs
if(is.null(x = cpgs)){
cpgs <- data.frame(chr = GenomeInfoDb::seqlevelsInUse(dat_gr[min(S4Vectors::queryHits(hit))]),
start = min(IRanges::start(dat_gr[S4Vectors::queryHits(hit)])),
end = max(IRanges::end(dat_gr[S4Vectors::queryHits(hit)])),
annotation = as.character(annots[anno]))
} else {
cpgs <- rbind(cpgs,
data.frame(chr = GenomeInfoDb::seqlevelsInUse(dat_gr[min(S4Vectors::queryHits(hit))]),
start = min(IRanges::start(dat_gr[S4Vectors::queryHits(hit)])),
end = max(IRanges::end(dat_gr[S4Vectors::queryHits(hit)])),
annotation = as.character(annots[anno])))
}
}
}
# transform the windows in GRanges to create regions
regions_gr <- GenomicRanges::GRanges(seqnames = cpgs$chr,
IRanges::IRanges(start = cpgs$start,
end = cpgs$end))
# save the identified regions
listOfRegions <- list()
listOfRegions_names <- c()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[(dat$Chr == GenomeInfoDb::seqlevelsInUse(region) &
dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region)),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",as.character(region),
"] containing ",length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",as.character(region),
"] containing ",length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else {
chrom <- gsub(x = unique(region_dat$Chr),
pattern = "chr",
replacement = "")
listOfRegions_names <- c(listOfRegions_names,
paste(paste0("chr",chrom),
min(region_dat$Position),
max(region_dat$Position), sep = "_"))
# suppress the temporary 'Chr' column
region_dat$Chr <- NULL
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- listOfRegions_names
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genomic annotations
#' @description This function splits the methylation data into regions
#' based on the genomic annotation provided under the form of a 1-based BED file
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param bed character, path to the 1-based BED file containing the annotations
#' @param gap integer defining the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1} (meaning strict overlapping).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom IRanges IRanges reduce gaps start end shift
#' @importFrom GenomicRanges GRanges findOverlaps
#' @importFrom rtracklayer import
#' @importFrom GenomeInfoDb seqlevelsStyle seqinfo
#' @importFrom S4Vectors subjectHits queryHits
#'
#' @export
splitDataByBed <- function(dat, chr, bed, gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"genomic annotations")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
if(is.null(bed) || !is.character(bed) || !file.exists(bed))
Error("'bed' should a valid file path.")
# error handling in BED import
annots <- tryCatch(
rtracklayer::import(bed),
error=function(e) e
)
if(inherits(annots, "error"))
Error("'bed' should a path to a valid BED file.")
# the BED file is converted automatically by rtracklayer in 0-based,
# we have to return back to 1-based
IRanges::start(annots) <- IRanges::start(annots) - 1
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
listOfRegions <- splitDataByGRanges(dat = dat, chr = chr, annots = annots,
gap = gap, min.cpgs = min.cpgs,
max.cpgs = max.cpgs, verbose = verbose)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the chromatin states
#' @description This function splits the methylation data into regions
#' based on the chromatin states predicted by ChromHMM software
#' (Ernst and Kellis (2012)).
#' The annotations come from the Bioconductor package `annnotatr`.
#' Chromatin states determined by chromHMM are
#' available in hg19 for nine cell lines (Gm12878, H1hesc, Hepg2, Hmec, Hsmm,
#' Huvec, K562, Nhek, and Nhlf).
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID`(sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param cell.line character defining the cell line of interest. Nine cell
#' lines are available:
#' - \code{"gm12878"}: Lymphoblastoid cells GM12878,
#' - \code{"h1hesc"}: Embryonic cells H1 hESC,
#' - \code{"hepg2"}: Liver carcinoma HepG2,
#' - \code{"hmec"}, Mammary epithelial cells HMEC,
#' - \code{"hsmm"}, Skeletal muscle myoblasts HSMM,
#' - \code{"huvec"}: Umbilical vein endothelial HUVEC,
#' - \code{"k562"}: Myelogenous leukemia K562,
#' - \code{"nhek"}: Keratinocytes NHEK,
#' - \code{"nhlf"}: Normal human lung fibroblasts NHLF.
#' @param states character vector defining the chromatin states of interest
#' among the following available options:
#' - \code{"ActivePromoter"}: Active Promoter
#' - \code{"WeakPromoter"}: Weak Promoter
#' - \code{"PoisedPromoter"}: Poised Promoter
#' - \code{"StrongEnhancer"}: Strong Enhancer
#' - \code{"WeakEnhancer"}: Weak/poised Enhancer
#' - \code{"Insulator"}: Insulator
#' - \code{"TxnTransition"}: Transcriptional Transition
#' - \code{"TxnElongation"}: Transcriptional Elongation
#' - \code{"WeakTxn"}: Weak Transcribed
#' - \code{"Repressed"}: Polycomb-Repressed
#' - \code{"Heterochrom"}: Heterochromatin; low signal
#' - \code{"RepetitiveCNV"}: Repetitive/Copy Number Variation
#' Use \code{state="all"} to select all the states simultaneously.
#' @param gap this integer defines the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1}.
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom yaml yaml.load_file
#' @importFrom BiocManager install
#' @importFrom annotatr build_annotations
#' @importFrom utils installed.packages
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByChromatin(dat = RAdat.f,
#' cell.line = "huvec", chr = rep(x = "chr4", times = nrow(RAdat.f)),
#' states = "Insulator", verbose = FALSE)
#'
#' @export
splitDataByChromatin <- function(dat, chr, cell.line, states,
gap = -1, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"chromatin states")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
# load annotation conf
annots_conf <- yaml.load_file(input = system.file("annotations.yml",
package = 'SOMNiBUS'))
# set organism and build
organism <- "human"
build <- "hg19"
# check the cell line
chromatin_annots <- annots_conf[[organism]][[build]][["chromatin"]]
supported_types <- names(chromatin_annots)
if(is.null(cell.line) || !is.character(cell.line)) {
stop_msg <- paste0("'cell.line' argument is mandatory. Choose one among ",
"these supported cell lines: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
} else {
cell.line <- tolower(cell.line)
}
if(!cell.line %in% supported_types){
stop_msg <- paste0("'",cell.line,"' cell line is not supported. ",
"Choose one among these supported cell lines: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
}
# check the chromHMM states
supported_states <- names(chromatin_annots[[cell.line]])
if(is.null(states) || !is.character(states)) {
stop_msg <- paste0("'states' argument is mandatory. Choose among ",
"these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Error(stop_msg)
}
if(sum(states %in% supported_states) == 0){
stop_msg <- paste0("None of the requested states are supported. ",
"Choose among these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Error(stop_msg)
}
unsupported_states <- unique(states[!states %in% supported_states])
if(length(unsupported_states) > 0){
w_msg <- paste0("The following unsupported states will be ignored: '",
paste(unsupported_states, collapse = "', '"),"'. ",
"Choose among these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
# retrieve annotations
states <- unique(states[!states %in% unsupported_states])
listOfRegions <- lapply(X = states, FUN = function(state){
annot <- chromatin_annots[[cell.line]][[state]]
annot <- suppressMessages(annotatr::build_annotations(genome = build,
annotations = annot))
# suppress duplicates
annot <- unique(annot)
# findoverlap
annots <- splitDataByGRanges(dat = dat, chr = chr,
annots = annot, gap = gap,
min.cpgs = min.cpgs, max.cpgs = max.cpgs,
verbose = verbose)
# overwrite the name of each independent region
if(length(annots) > 0)
names(annots) <- paste(names(annots), cell.line, state, sep = "_")
return(annots)
})
# create final listOfRegions
listOfRegions <- unlist(listOfRegions,recursive=FALSE)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genes annotations
#' @description This function splits the methylation data into regions
#' based on the genes. The annotations are coming from the Bioconductor
#' package `annnotatr`.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID`(sample ID).
#' The covariate information, such as disease status or cell type composition,
#' are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param organism character defining the organism of interest
#' Only Homo sapiens (\code{"human"}) is available.
#' Additional packages are required for Mus musculus (\code{"mouse"}),
#' Rattus norvegicus (\code{"rat"}) and Drosophila melanogaster (\code{"fly"}).
#' The matching is case-insensitive. The default value is \code{"human"}.
#' @param build character defining the version of the genome build on which the
#' methylation data have been mapped. By default, the build is set to
#' \code{"hg38"}, however the build \code{"hg19"} is also available for
#' Homo sapiens:
#' Once the additional packages are installed, the following organisms and
#' builds are available:
#' - \code{"mm9"} and \code{"mm10"} for Mus musculus;
#' - \code{"rn4"}, \code{"rn5"} and \code{"rn6"} for Rattus norvegicus;
#' - \code{"dm3"} and \code{"dm6"} for Drosophila melanogaster;
#' @param types character vector defining the type of genic annotations
#' to use among the following options:
#' - \code{"upstream"} for the annotations included 1-5Kb upstream of the TSS;
#' - \code{"promoter"} for the annotations included < 1Kb upstream of the TSS;
#' - \code{"threeprime"} for the annotations included in 3' UTR;
#' - \code{"fiveprime"} for the annotations included in the 5' UTR;
#' - \code{"exon"} for the annotations included in the exons;
#' - \code{"intron"} for the annotations included in the introns;
#' - \code{"all"} for all the annotations aforementioned.
#' The default value is \code{"promoter"}.
#' @param gap this integer defines the maximum gap allowed between two regions
#' to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1}.
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom yaml yaml.load_file
#' @importFrom BiocManager install
#' @importFrom annotatr build_annotations
#' @importFrom utils installed.packages
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' # Add a column containing the chromosome information
#' RAdat$Chr <- "chr4"
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByGene(dat = RAdat.f,
#' chr = rep(x = "chr1", times = nrow(RAdat.f)), verbose = FALSE)
#'
#' @export
splitDataByGene <- function(dat, chr, organism = "human",
build = "hg38", types = "promoter", gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"genic annotations")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
# load annotation conf
annots_conf <- yaml.load_file(input = system.file("annotations.yml",
package = 'SOMNiBUS'))
# check organism
if(is.null(organism)){
organism <- "human"
build <- "hg38"
} else{
if(!is.character(organism)) {
organism <- "human"
build <- "hg38"
w_msg <- paste("'organism' should be a non-null character.",
"We will use the default value human",
"and the build will be set at hg38.")
Warning(w_msg)
} else {
organism <- tolower(organism)
}
}
if(!organism %in% names(annots_conf)){
stop_msg <- paste0("'",organism,"' organism is not supported. ",
"Choose one among these supported types: '",
paste(names(annots_conf), collapse = "', '"),"'.")
Error(stop_msg)
}
# check build
if(is.null(build)){
build <- annots_conf[[organism]]$default
} else {
if(!is.character(build)) {
build <- annots_conf[[organism]]$default
warning_msg <- paste0("'build' should be a non-null character.",
"We will use the default",
"value for ",organism, ": ", build,".")
Warning(warning_msg)
} else {
build <- tolower(build)
}
}
supported_builds <- names(annots_conf[[organism]])
supported_builds <- supported_builds[!supported_builds %in% "default"]
if(!build %in% supported_builds){
stop_msg <- paste0("'",build,"' build is not supported. Choose one among these ",
"supported builds: '",
paste(supported_builds, collapse = "', '"),"'.")
Error(stop_msg)
}
# check the annotation types
if(is.null(types)){
types <- "promoter"
} else {
if(!is.character(types)) {
types <- "promoter"
w_msg <- paste("'type' should be a non-null character.",
"We will use the default value promoter.")
Warning(w_msg)
} else {
types <- tolower(types)
}
}
supported_types <- names(annots_conf[[organism]][[build]][["genes"]])
supported_types <- supported_types[!supported_types %in% "dependencies"]
if(sum(types %in% supported_types) == 0){
stop_msg <- paste0("None of the requested types are supported. ",
"Choose among these supported types: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
}
unsupported_types <- unique(types[!types %in% supported_types])
if(length(unsupported_types) > 0){
w_msg <- paste0("The following unsupported types will be ignored: '",
paste(unsupported_types, collapse = "', '"),"'. ",
"Choose among these supported types: '",
paste(supported_types, collapse = "', '"),"'.")
Warning(w_msg)
}
# retrieve annotations
types <- unique(types[!types %in% unsupported_types])
# check if the annotation package is installed, if not install stop
annotation_packages <- annots_conf[[organism]][[build]]$genes$dependencies
missing_packages <- annotation_packages[! annotation_packages %in%
installed.packages()[,"Package"]]
if (length(missing_packages) > 0){
stop_msg <- paste0("Please install the following package(s) in order to ",
"perform the requested partiniong: '",
paste(missing_packages, collapse = "', '"),"'.")
Error(stop_msg)
}
listOfRegions <- lapply(X = types, FUN = function(type){
# retrieve annotations
annot <- annots_conf[[organism]][[build]]$genes[[type]]
annot <- suppressMessages(annotatr::build_annotations(genome = build,
annotations = annot))
# suppress duplicates
annot <- unique(annot)
# findoverlap
annots <- splitDataByGRanges(dat = dat, chr = chr, annots = annot,
gap = gap, min.cpgs = min.cpgs,
max.cpgs = max.cpgs, verbose = verbose)
# overwrite the name of each independent region
if(length(annots) > 0)
names(annots) <- paste(names(annots), type, sep = "_")
return(annots)
})
# create final listOfRegions
listOfRegions <- unlist(listOfRegions,recursive=FALSE)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
kaiqiong/SOMNiBUS documentation built on Feb. 24, 2023, 5:38 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions splitDataByGene splitDataByChromatin splitDataByBed splitDataByGRanges splitDataByDensity splitDataByRegion hessianComp
Documented in splitDataByBed splitDataByChromatin splitDataByDensity splitDataByGene splitDataByGRanges splitDataByRegion
#' @title hessianComp computes the Hessian matrix
#'
#' @description computes the Hessian matrix for an EM estimate
#' @param w_ij the diagnal values of the weight matrix
#' @param new.par estimate of alpha
#' @param new.lambda estimate of lambda
#' @param X a vector of read depths
#' @param Y a vector of methylated counts
#' @param my.design.matrix design matrix from the final fit
#' @param gam_smoothMat the smooth matrix from the final gam fit
#' @param Z covariate matrix
#' @param pred.pi predicted methylation probability from the final fit
#' @param p0 the probability of observing a methylated read when
#' the underlying true status is unmethylated. \code{p0} is the rate of false
#' methylation calls, i.e. false positive rate.
#' @param p1 the probability of observing a methylated read when
#' the underlying true status is methylated. \code{1-p1} is the rate of false
#' non-methylation calls, i.e. false negative rate.
#' @param disp_est estimated dispersion parameter
#' @param RanEff whether a subject-level Random effect is added or not
#' @param N number of unique samples in the provided dataset
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return a Hessian matrix for all alphas (including the alphas for RE)
#' @references Zhao, Kaiqiong, et al. 'A novel statistical method for
#' modeling covariate effects in bisulfite sequencing derived measures
#' of DNA methylation.'Biometrics (2020).
#' @author Kaiqiong Zhao
#' @importFrom Matrix bdiag
#'
#' @noRd
hessianComp <- function(w_ij, new.par, new.lambda, X, Y, my.design.matrix,
gam_smoothMat, Z, pred.pi, p0, p1, disp_est,
RanEff, N, verbose = TRUE) {
t0 <- Sys.time()
if(verbose) Message("Compute the Hessian matrix for an EM estimate")
## Q1: the second partial derivative w.r.t alpha^2 Q2: the second
## derivative w.r.t alpha & alpha_star
res <- outer(seq_len(length(new.par)),
seq_len(length(new.par)),
Vectorize(function(l,m) {
sum(-X * w_ij * my.design.matrix[, m] * my.design.matrix[, l])
}))
smoth.mat <- lapply(as.list(seq_len(ncol(Z) + 1)), function(i) {
gam_smoothMat[[i]]$S[[1]] * new.lambda[i]
}) ## extract the penalty matrix
## assume the lambda for the constant of the intercept is 0 -- no
## penalization
smoth.mat[[length(smoth.mat) + 1]] <- 0
if (RanEff) {
## !!!! Otherwise, we get very wide CI
smoth.mat[[length(smoth.mat) + 1]] <- diag(N) * new.lambda[ncol(Z) + 2]
span.penal.matrix <- as.matrix(Matrix::bdiag(smoth.mat[c(length(smoth.mat) - 1,
(seq_len((length(smoth.mat) - 2))),
length(smoth.mat))]))
} else {
span.penal.matrix <- as.matrix(Matrix::bdiag(smoth.mat[c(length(smoth.mat),
(seq_len((length(smoth.mat) - 1))))]))
}
Q1_with_lambda <- res - span.penal.matrix/disp_est
Q1_no_lambda <- res
Q2 <- outer(seq_len(length(new.par)), seq_len(length(new.par)),
Vectorize(function(l, m) {
term1 <- Y * p1 * p0/(p1 * pred.pi + p0 * (1 - pred.pi))^2 + (X - Y) *
(1 - p1) * (1 - p0)/((1 - p1) * pred.pi + (1 - p0) * (1 - pred.pi))^2
sum(term1 * w_ij * my.design.matrix[, m] * my.design.matrix[, l])
}))
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(Q1_with_lambda + Q2)
}
#' @title Split methylation data into regions based on the spacing of CpGs
#' @description This function splits the methylation data into regions
#' based on the spacing of CpGs.
#'
#' @param dat a data frame with rows as individual CpGs appearing in all the
#' samples. The first 4 columns should contain the information of `Meth_Counts`
#' (methylated counts), `Total_Counts` (read depths), `Position` (Genomic
#' position for the CpG site) and `ID` (sample ID). The covariate information,
#' such as disease status or cell type composition, are listed in column 5 and
#' onwards.
#' @param gap positive integer defining the gap width
#' beyond which we consider that two regions are independent.
#' Odd and decimal values will be rounded to the next even numbers
#' (e.g. 8.2 and 8.7 become gaps of 8 and 10 respectively).
#' The default value is \code{1e+6} (1Mb).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByRegion( dat=RAdat.f, gap = 1e6, min.cpgs = 5,
#' verbose = FALSE)
#' @importFrom IRanges IRanges reduce gaps start end
#' @export
splitDataByRegion <- function(dat, gap = 1e6, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste("Partitioning the regions based on",
"the spacing of CpGs")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- 1e+6
} else {
if(!is.numeric(gap) || gap < 1){
gap <- 1e+6
w_msg <- paste("'gap' should be a positive integer.",
"We will use the default value 1e+6 (1Mb).")
Warning(w_msg)
} else {
ori_gap <- gap
gap <- round(gap)
if(gap %% 2 != 0) {
gap <- gap + 1
}
}
}
# convert dataframe in IRanges
dat_gr <- IRanges::IRanges(start = dat$Position, end = dat$Position)
# reduce first orders the ranges in x from left to right,
# then merges the overlapping or adjacent ones.
dat_gr <- IRanges::reduce(dat_gr)
# identify gap between regions
g = IRanges::gaps(dat_gr)
# extend each region with half of the maximal distance d
# and attempt reduction to identify overlapping regions
hemi_d <- ceiling(gap/2)
regions_gr <- IRanges::reduce(dat_gr + hemi_d) - hemi_d
# save the identified regions
listOfRegions <- list()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else{
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- sapply(X = listOfRegions,
FUN = function(region)
paste("region",
min(region$Position),
max(region$Position), sep = "_")
)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
# return the list of independent regions
return(listOfRegions)
}
#' @title Split methylation data into regions based on the density of CpGs
#' @description This function splits the methylation data into regions
#' based on the density of CpGs.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status or cell type composition,
#' are listed in column 5 and onwards.
#' @param window.size this positive integer defines the size of the
#' sliding window in bp. Decimal values will be rounded to the nearest integer.
#' The value should be greater than 10. The default value is \code{100} (100 bp)
#' @param by positive integer defines by how many base pairs the
#' window moves at each increment. Decimal values will be rounded to the
#' nearest integer. The default value is \code{1} (1 bp).
#' @param min.density positive integer defines the minimum density
#' threshold for each window. Decimal values will be rounded to the
#' nearest integer. The default value is \code{5} (5 CpGs/\code{window.size}).
#' @param gap positive integer defining the gap width
#' beyond which we consider that two regions are independent.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{10} (10bp).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information. The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByDensity(dat = RAdat.f, window.size = 100, by = 1,
#' min.density = 5, gap = 10, min.cpgs = 50, verbose = FALSE)
#' @importFrom IRanges IRanges reduce start end shift
#'
#' @export
splitDataByDensity <- function(dat, window.size = 100, by = 1, min.density = 5,
gap = 10, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"the density of CpGs")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(window.size)){
window.size <- 100
} else {
if(!is.numeric(window.size) || window.size < 10) {
window.size <- 100
w_msg <- paste("'window.size' should be a positive integer",
"greater than 10. We will use the default value 100.")
Warning(w_msg)
} else {
window.size <- round(window.size)
}
}
if(is.null(by)){
by <- 1
} else {
if(!is.numeric(by) || by < 1) {
by <- 1
w_msg <- paste("'by' should be a positive integer.",
"We will use the default value 1.")
Warning(w_msg)
} else {
by <- round(by)
}
}
if(is.null(min.density)){
min.density <- 5
} else {
if(!is.numeric(min.density) || min.density < 1) {
min.density <- 5
w_msg <- paste("'min.density' should be a positive integer.",
"We will use the default value 5.")
Warning(w_msg)
} else {
min.density <- round(min.density)
}
}
if(is.null(gap)){
gap <- 10
} else {
if(!is.numeric(gap) || gap < 1) {
gap <- 10
w_msg <- paste("'gap' should be a positive integer.",
"We will use the default value 10.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
if(gap >= window.size){
gap <- round(window.size/2)
w_msg <- paste0("'gap' should be lesser than the 'window.size' value. ",
"We will use half of window size (",gap,").")
Warning(w_msg)
}
# convert dataframe in IRanges
dat_gr <- unique(IRanges::IRanges(start = dat$Position, end = dat$Position))
# save the end index for the while loop
endOfRegion <- max(IRanges::end(dat_gr))
# create the starting sliding window
sliding_window <- IRanges::IRanges(start = min(IRanges::start(dat_gr)),
end = (min(IRanges::start(dat_gr)) +
window.size - 1))
cpgs <- NULL
# calculate the density in each sliding window
while(IRanges::start(sliding_window) <= endOfRegion){
# overlap region to window
hits <- IRanges::findOverlaps(query = dat_gr, subject = sliding_window)
# save density
if(is.null(x = cpgs)){
cpgs <- data.frame(start = IRanges::start(sliding_window),
end = IRanges::end(sliding_window),
density = length(hits))
} else {
cpgs <- rbind(cpgs,
data.frame(start = IRanges::start(sliding_window),
end = IRanges::end(sliding_window),
density = length(hits)))
}
# shift window
sliding_window <- IRanges::shift(x = sliding_window, shift = by)
}
# filter on density
cpgs <- cpgs[cpgs$density >= min.density,]
# transform the windows in IRanges to create regions
regions_gr <- IRanges::IRanges(start = cpgs$start, end = cpgs$end)
# each window is resized to identify independent regions
regions_gr <- IRanges::resize(x = regions_gr, width = (gap + 1))
regions_gr <- IRanges::reduce(regions_gr)
# save the identified regions
listOfRegions <- list()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",IRanges::start(region),"-",
IRanges::end(region),"] containing ",
length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else{
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- sapply(X = listOfRegions,
FUN = function(region)
paste("region",
min(region$Position),
max(region$Position), sep = "_")
)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genomic annotations
#' @description This function splits the methylation data into regions
#' based on the genomic annotations provided under the form of a GenomicRanges
#' object.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param annots GenomicRanges object containing the annotations
#' @param gap integer defining the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1} (meaning strict overlaping).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' annot <- GenomicRanges::GRanges(seqnames = "chr1", IRanges::IRanges(
#' start = c(102711720,102711844,102712006,102712503,102712702),
#' end = c(102711757,102711909,102712195,102712637,102712712)
#' ))
#' results <- splitDataByGRanges(dat = RAdat.f,
#' chr = rep(x = "chr1", times = nrow(RAdat.f)),
#' annots = annot, gap = -1, min.cpgs = 5)
#' @importFrom IRanges IRanges reduce gaps start end shift
#' @importFrom GenomicRanges GRanges findOverlaps
#' @importFrom rtracklayer import
#' @importFrom GenomeInfoDb seqlevelsStyle seqinfo
#' @importFrom S4Vectors subjectHits queryHits
#'
#' @export
splitDataByGRanges <- function(dat, chr, annots, gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
if(!inherits(annots, "GenomicRanges"))
Error("'annots' should inherit the GRanges class.")
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
# convert dataframe in GRanges
dat$Chr <- chr
dat_gr <- GenomicRanges::GRanges(seqnames = dat$Chr,
ranges = IRanges::IRanges(start = dat$Position,
end = dat$Position))
dat_gr <- unique(dat_gr)
# error handling
dat_style <- tryCatch({
# check if style has a compatible entry for the species supported by Seqname
GenomeInfoDb::seqlevelsStyle(dat_gr)
}, error=function(e) e)
if(inherits(dat_style, "error")){
stop_msg <- paste0("Incompatible chromosome name. ",
"Please see genomeStyles() for supported styles.")
Error(stop_msg)
}
# error handling
hits <- tryCatch({
# adapt the style (chromosome nomenclature) of the annotation
# to the style of the input data if needed
if(length(intersect(x = GenomeInfoDb::seqlevelsStyle(annots),
y = dat_style)) == 0)
suppressWarnings(GenomeInfoDb::seqlevelsStyle(annots) <- dat_style[1])
# find overlaps
GenomicRanges::findOverlaps(query = dat_gr, subject = annots, maxgap = gap)
}, error=function(e) e)
if(inherits(hits, "error")){
stop_msg <- paste0("Incompatible formats between ",
"the input data and the annotations.")
Error(stop_msg)
}
cpgs <- NULL
if(length(hits) > 0){
for(anno in unique(S4Vectors::subjectHits(hits))){
hit <- hits[S4Vectors::subjectHits(hits) == anno]
# save cpgs
if(is.null(x = cpgs)){
cpgs <- data.frame(chr = GenomeInfoDb::seqlevelsInUse(dat_gr[min(S4Vectors::queryHits(hit))]),
start = min(IRanges::start(dat_gr[S4Vectors::queryHits(hit)])),
end = max(IRanges::end(dat_gr[S4Vectors::queryHits(hit)])),
annotation = as.character(annots[anno]))
} else {
cpgs <- rbind(cpgs,
data.frame(chr = GenomeInfoDb::seqlevelsInUse(dat_gr[min(S4Vectors::queryHits(hit))]),
start = min(IRanges::start(dat_gr[S4Vectors::queryHits(hit)])),
end = max(IRanges::end(dat_gr[S4Vectors::queryHits(hit)])),
annotation = as.character(annots[anno])))
}
}
}
# transform the windows in GRanges to create regions
regions_gr <- GenomicRanges::GRanges(seqnames = cpgs$chr,
IRanges::IRanges(start = cpgs$start,
end = cpgs$end))
# save the identified regions
listOfRegions <- list()
listOfRegions_names <- c()
for(i in seq_along(regions_gr)){
region = regions_gr[i]
region_dat <- dat[(dat$Chr == GenomeInfoDb::seqlevelsInUse(region) &
dat$Position >= IRanges::start(region) &
dat$Position <= IRanges::end(region)),]
if(length(unique(region_dat$Position)) > max.cpgs){
msg <- paste0("Dropped region [",as.character(region),
"] containing ",length(unique(region_dat$Position)),
" measurement points: above maximal ",
"region size accepted for analysis (",max.cpgs,").")
if(verbose) Message(msg)
} else {
if(length(unique(region_dat$Position)) < min.cpgs) {
msg <- paste0("Dropped region [",as.character(region),
"] containing ",length(unique(region_dat$Position)),
" measurement points: below minimal ",
"region size accepted for analysis (",min.cpgs,").")
if(verbose) Message(msg)
} else {
chrom <- gsub(x = unique(region_dat$Chr),
pattern = "chr",
replacement = "")
listOfRegions_names <- c(listOfRegions_names,
paste(paste0("chr",chrom),
min(region_dat$Position),
max(region_dat$Position), sep = "_"))
# suppress the temporary 'Chr' column
region_dat$Chr <- NULL
listOfRegions[[length(listOfRegions) + 1]] <- region_dat
}
}
}
# name each independent region
names(listOfRegions) <- listOfRegions_names
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genomic annotations
#' @description This function splits the methylation data into regions
#' based on the genomic annotation provided under the form of a 1-based BED file
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID` (sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param bed character, path to the 1-based BED file containing the annotations
#' @param gap integer defining the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1} (meaning strict overlapping).
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom IRanges IRanges reduce gaps start end shift
#' @importFrom GenomicRanges GRanges findOverlaps
#' @importFrom rtracklayer import
#' @importFrom GenomeInfoDb seqlevelsStyle seqinfo
#' @importFrom S4Vectors subjectHits queryHits
#'
#' @export
splitDataByBed <- function(dat, chr, bed, gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"genomic annotations")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
if(is.null(bed) || !is.character(bed) || !file.exists(bed))
Error("'bed' should a valid file path.")
# error handling in BED import
annots <- tryCatch(
rtracklayer::import(bed),
error=function(e) e
)
if(inherits(annots, "error"))
Error("'bed' should a path to a valid BED file.")
# the BED file is converted automatically by rtracklayer in 0-based,
# we have to return back to 1-based
IRanges::start(annots) <- IRanges::start(annots) - 1
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
listOfRegions <- splitDataByGRanges(dat = dat, chr = chr, annots = annots,
gap = gap, min.cpgs = min.cpgs,
max.cpgs = max.cpgs, verbose = verbose)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the chromatin states
#' @description This function splits the methylation data into regions
#' based on the chromatin states predicted by ChromHMM software
#' (Ernst and Kellis (2012)).
#' The annotations come from the Bioconductor package `annnotatr`.
#' Chromatin states determined by chromHMM are
#' available in hg19 for nine cell lines (Gm12878, H1hesc, Hepg2, Hmec, Hsmm,
#' Huvec, K562, Nhek, and Nhlf).
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID`(sample ID).
#' The covariate information, such as disease status
#' or cell type composition, are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param cell.line character defining the cell line of interest. Nine cell
#' lines are available:
#' - \code{"gm12878"}: Lymphoblastoid cells GM12878,
#' - \code{"h1hesc"}: Embryonic cells H1 hESC,
#' - \code{"hepg2"}: Liver carcinoma HepG2,
#' - \code{"hmec"}, Mammary epithelial cells HMEC,
#' - \code{"hsmm"}, Skeletal muscle myoblasts HSMM,
#' - \code{"huvec"}: Umbilical vein endothelial HUVEC,
#' - \code{"k562"}: Myelogenous leukemia K562,
#' - \code{"nhek"}: Keratinocytes NHEK,
#' - \code{"nhlf"}: Normal human lung fibroblasts NHLF.
#' @param states character vector defining the chromatin states of interest
#' among the following available options:
#' - \code{"ActivePromoter"}: Active Promoter
#' - \code{"WeakPromoter"}: Weak Promoter
#' - \code{"PoisedPromoter"}: Poised Promoter
#' - \code{"StrongEnhancer"}: Strong Enhancer
#' - \code{"WeakEnhancer"}: Weak/poised Enhancer
#' - \code{"Insulator"}: Insulator
#' - \code{"TxnTransition"}: Transcriptional Transition
#' - \code{"TxnElongation"}: Transcriptional Elongation
#' - \code{"WeakTxn"}: Weak Transcribed
#' - \code{"Repressed"}: Polycomb-Repressed
#' - \code{"Heterochrom"}: Heterochromatin; low signal
#' - \code{"RepetitiveCNV"}: Repetitive/Copy Number Variation
#' Use \code{state="all"} to select all the states simultaneously.
#' @param gap this integer defines the maximum gap that is allowed between
#' two regions to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1}.
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom yaml yaml.load_file
#' @importFrom BiocManager install
#' @importFrom annotatr build_annotations
#' @importFrom utils installed.packages
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByChromatin(dat = RAdat.f,
#' cell.line = "huvec", chr = rep(x = "chr4", times = nrow(RAdat.f)),
#' states = "Insulator", verbose = FALSE)
#'
#' @export
splitDataByChromatin <- function(dat, chr, cell.line, states,
gap = -1, min.cpgs = 50, max.cpgs = 2000,
verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"chromatin states")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
# load annotation conf
annots_conf <- yaml.load_file(input = system.file("annotations.yml",
package = 'SOMNiBUS'))
# set organism and build
organism <- "human"
build <- "hg19"
# check the cell line
chromatin_annots <- annots_conf[[organism]][[build]][["chromatin"]]
supported_types <- names(chromatin_annots)
if(is.null(cell.line) || !is.character(cell.line)) {
stop_msg <- paste0("'cell.line' argument is mandatory. Choose one among ",
"these supported cell lines: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
} else {
cell.line <- tolower(cell.line)
}
if(!cell.line %in% supported_types){
stop_msg <- paste0("'",cell.line,"' cell line is not supported. ",
"Choose one among these supported cell lines: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
}
# check the chromHMM states
supported_states <- names(chromatin_annots[[cell.line]])
if(is.null(states) || !is.character(states)) {
stop_msg <- paste0("'states' argument is mandatory. Choose among ",
"these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Error(stop_msg)
}
if(sum(states %in% supported_states) == 0){
stop_msg <- paste0("None of the requested states are supported. ",
"Choose among these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Error(stop_msg)
}
unsupported_states <- unique(states[!states %in% supported_states])
if(length(unsupported_states) > 0){
w_msg <- paste0("The following unsupported states will be ignored: '",
paste(unsupported_states, collapse = "', '"),"'. ",
"Choose among these supported states: '",
paste(supported_states, collapse = "', '"),"'.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
# retrieve annotations
states <- unique(states[!states %in% unsupported_states])
listOfRegions <- lapply(X = states, FUN = function(state){
annot <- chromatin_annots[[cell.line]][[state]]
annot <- suppressMessages(annotatr::build_annotations(genome = build,
annotations = annot))
# suppress duplicates
annot <- unique(annot)
# findoverlap
annots <- splitDataByGRanges(dat = dat, chr = chr,
annots = annot, gap = gap,
min.cpgs = min.cpgs, max.cpgs = max.cpgs,
verbose = verbose)
# overwrite the name of each independent region
if(length(annots) > 0)
names(annots) <- paste(names(annots), cell.line, state, sep = "_")
return(annots)
})
# create final listOfRegions
listOfRegions <- unlist(listOfRegions,recursive=FALSE)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
#' @title Split methylation data into regions based on the genes annotations
#' @description This function splits the methylation data into regions
#' based on the genes. The annotations are coming from the Bioconductor
#' package `annnotatr`.
#'
#' @param dat a data frame with rows as individual CpGs appearing
#' in all the samples. The first 4 columns should contain the information of
#' `Meth_Counts` (methylated counts), `Total_Counts` (read depths),
#' `Position` (Genomic position for the CpG site) and `ID`(sample ID).
#' The covariate information, such as disease status or cell type composition,
#' are listed in column 5 and onwards.
#' @param chr character vector containing the chromosome information. Its length
#' should be equal to the number of rows in \code{dat}.
#' @param organism character defining the organism of interest
#' Only Homo sapiens (\code{"human"}) is available.
#' Additional packages are required for Mus musculus (\code{"mouse"}),
#' Rattus norvegicus (\code{"rat"}) and Drosophila melanogaster (\code{"fly"}).
#' The matching is case-insensitive. The default value is \code{"human"}.
#' @param build character defining the version of the genome build on which the
#' methylation data have been mapped. By default, the build is set to
#' \code{"hg38"}, however the build \code{"hg19"} is also available for
#' Homo sapiens:
#' Once the additional packages are installed, the following organisms and
#' builds are available:
#' - \code{"mm9"} and \code{"mm10"} for Mus musculus;
#' - \code{"rn4"}, \code{"rn5"} and \code{"rn6"} for Rattus norvegicus;
#' - \code{"dm3"} and \code{"dm6"} for Drosophila melanogaster;
#' @param types character vector defining the type of genic annotations
#' to use among the following options:
#' - \code{"upstream"} for the annotations included 1-5Kb upstream of the TSS;
#' - \code{"promoter"} for the annotations included < 1Kb upstream of the TSS;
#' - \code{"threeprime"} for the annotations included in 3' UTR;
#' - \code{"fiveprime"} for the annotations included in the 5' UTR;
#' - \code{"exon"} for the annotations included in the exons;
#' - \code{"intron"} for the annotations included in the introns;
#' - \code{"all"} for all the annotations aforementioned.
#' The default value is \code{"promoter"}.
#' @param gap this integer defines the maximum gap allowed between two regions
#' to be considered as overlapping.
#' According to the \code{GenomicRanges::findOverlaps} function,
#' the gap between 2 ranges is the number of positions that separate them.
#' The gap between 2 adjacent ranges is 0. By convention when one range has
#' its start or end strictly inside the other (i.e. non-disjoint ranges),
#' the gap is considered to be -1.
#' Decimal values will be rounded to the nearest integer.
#' The default value is \code{-1}.
#' @param min.cpgs positive integer defining the minimum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 50.
#' @param max.cpgs positive integer defining the maximum number of
#' CpGs within a region for the algorithm to perform optimally.
#' The default value is 2000.
#' @param verbose logical indicates if the algorithm should provide progress
#' report information.
#' The default value is TRUE.
#' @return A named \code{list} of \code{data.frame} containing the data of each
#' independent region.
#' @author Audrey Lemaçon
#'
#' @importFrom yaml yaml.load_file
#' @importFrom BiocManager install
#' @importFrom annotatr build_annotations
#' @importFrom utils installed.packages
#' @examples
#' #------------------------------------------------------------#
#' data(RAdat)
#' # Add a column containing the chromosome information
#' RAdat$Chr <- "chr4"
#' RAdat.f <- na.omit(RAdat[RAdat$Total_Counts != 0, ])
#' results <- splitDataByGene(dat = RAdat.f,
#' chr = rep(x = "chr1", times = nrow(RAdat.f)), verbose = FALSE)
#'
#' @export
splitDataByGene <- function(dat, chr, organism = "human",
build = "hg38", types = "promoter", gap = -1,
min.cpgs = 50, max.cpgs = 2000, verbose = TRUE){
t0 <- Sys.time()
msg <- paste0("Partitioning the regions based on ",
"genic annotations")
if(verbose) Message(msg)
# test the arguments
if(is.null(min.cpgs)){
min.cpgs <- 50
} else {
if(!is.numeric(min.cpgs) || min.cpgs < 1) {
min.cpgs <- 50
w_msg <- paste("'min.cpgs' should be a positive integer.",
"We will use the default value 50.")
Warning(w_msg)
} else {
min.cpgs <- round(min.cpgs)
}
}
if(is.null(max.cpgs)){
max.cpgs <- 2000
} else {
if(!is.numeric(max.cpgs) || max.cpgs < 1) {
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a positive integer.",
"We will use the default value 2000.")
Warning(w_msg)
} else {
max.cpgs <- round(max.cpgs)
}
}
if(max.cpgs <= min.cpgs){
min.cpgs <- 50
max.cpgs <- 2000
w_msg <- paste("'max.cpgs' should be a higher than 'min.cpgs'.",
"We will use the default values:",
"min.cpgs=50 and max.cpgs=2000.")
Warning(w_msg)
}
if(is.null(gap)){
gap <- -1
} else {
if(!is.numeric(gap) || gap < -1) {
gap <- -1
w_msg <- paste("'gap' should be an integer >= -1.",
"We will use the default value -1.")
Warning(w_msg)
} else {
gap <- round(gap)
}
}
if(is.null(chr) || length(chr) != nrow(dat)) {
e_msg <- paste("'chr' should be a character vector which length equals the",
"number of rows in 'dat'.")
Error(e_msg)
}
# load annotation conf
annots_conf <- yaml.load_file(input = system.file("annotations.yml",
package = 'SOMNiBUS'))
# check organism
if(is.null(organism)){
organism <- "human"
build <- "hg38"
} else{
if(!is.character(organism)) {
organism <- "human"
build <- "hg38"
w_msg <- paste("'organism' should be a non-null character.",
"We will use the default value human",
"and the build will be set at hg38.")
Warning(w_msg)
} else {
organism <- tolower(organism)
}
}
if(!organism %in% names(annots_conf)){
stop_msg <- paste0("'",organism,"' organism is not supported. ",
"Choose one among these supported types: '",
paste(names(annots_conf), collapse = "', '"),"'.")
Error(stop_msg)
}
# check build
if(is.null(build)){
build <- annots_conf[[organism]]$default
} else {
if(!is.character(build)) {
build <- annots_conf[[organism]]$default
warning_msg <- paste0("'build' should be a non-null character.",
"We will use the default",
"value for ",organism, ": ", build,".")
Warning(warning_msg)
} else {
build <- tolower(build)
}
}
supported_builds <- names(annots_conf[[organism]])
supported_builds <- supported_builds[!supported_builds %in% "default"]
if(!build %in% supported_builds){
stop_msg <- paste0("'",build,"' build is not supported. Choose one among these ",
"supported builds: '",
paste(supported_builds, collapse = "', '"),"'.")
Error(stop_msg)
}
# check the annotation types
if(is.null(types)){
types <- "promoter"
} else {
if(!is.character(types)) {
types <- "promoter"
w_msg <- paste("'type' should be a non-null character.",
"We will use the default value promoter.")
Warning(w_msg)
} else {
types <- tolower(types)
}
}
supported_types <- names(annots_conf[[organism]][[build]][["genes"]])
supported_types <- supported_types[!supported_types %in% "dependencies"]
if(sum(types %in% supported_types) == 0){
stop_msg <- paste0("None of the requested types are supported. ",
"Choose among these supported types: '",
paste(supported_types, collapse = "', '"),"'.")
Error(stop_msg)
}
unsupported_types <- unique(types[!types %in% supported_types])
if(length(unsupported_types) > 0){
w_msg <- paste0("The following unsupported types will be ignored: '",
paste(unsupported_types, collapse = "', '"),"'. ",
"Choose among these supported types: '",
paste(supported_types, collapse = "', '"),"'.")
Warning(w_msg)
}
# retrieve annotations
types <- unique(types[!types %in% unsupported_types])
# check if the annotation package is installed, if not install stop
annotation_packages <- annots_conf[[organism]][[build]]$genes$dependencies
missing_packages <- annotation_packages[! annotation_packages %in%
installed.packages()[,"Package"]]
if (length(missing_packages) > 0){
stop_msg <- paste0("Please install the following package(s) in order to ",
"perform the requested partiniong: '",
paste(missing_packages, collapse = "', '"),"'.")
Error(stop_msg)
}
listOfRegions <- lapply(X = types, FUN = function(type){
# retrieve annotations
annot <- annots_conf[[organism]][[build]]$genes[[type]]
annot <- suppressMessages(annotatr::build_annotations(genome = build,
annotations = annot))
# suppress duplicates
annot <- unique(annot)
# findoverlap
annots <- splitDataByGRanges(dat = dat, chr = chr, annots = annot,
gap = gap, min.cpgs = min.cpgs,
max.cpgs = max.cpgs, verbose = verbose)
# overwrite the name of each independent region
if(length(annots) > 0)
names(annots) <- paste(names(annots), type, sep = "_")
return(annots)
})
# create final listOfRegions
listOfRegions <- unlist(listOfRegions,recursive=FALSE)
msg <- paste("Process completed in",
format(Sys.time() - t0, digits = 2))
if(verbose) Message(msg, step = "Finished")
return(listOfRegions)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.