#' create annotation
#'
#' This function creates a GRanges object from an input rnb.set object which then is input
#' for the IHS score calculation
#'
#' @param rnb.path path to an RnBSet object stored in disk
#' @param use.sex.chromosomes Flag indicating if the sex chromosomes are to be included into the analysis
#'
#' @return GRanges object containing the CpG sites of interest
#' @noRd
create.annotation <- function(rnb.path,use.sex.chromosomes=FALSE){
if(inherits(rnb.path,"RnBSet")){
rnb.set <- rnb.path
}else if(file.exists(rnb.path)){
rnb.set <- load.rnb.set(rnb.path)
}else{
stop("Invalid value for rnb.set")
}
anno <- annotation(rnb.set)
if(hasCovg(rnb.set)){
coverage <- covg(rnb.set)
mean_coverage <- rowMeans(coverage,na.rm=TRUE)
rm(coverage)
high_enough <- mean_coverage >= get.option('coverage.threshold')
anno <- anno[high_enough,]
}
if(!use.sex.chromosomes){
if(all(c("chrX","chrY")%in%anno$Chromosome)){
anno <- anno[anno$Chromosome!="chrX",]
anno <- anno[anno$Chromosome!="chrY",]
}else if(all(c("X","Y")%in%anno$Chromosome)){
anno <- anno[anno$Chromosome!="X",]
anno <- anno[anno$Chromosome!="Y",]
}
}
anno <- GRanges(Rle(anno$Chromosome),IRanges(start=anno$Start,end=anno$End),anno$Strand)
names <- paste(as.character(seqnames(anno)),start(ranges(anno)),end(ranges(anno)),sep="_")
unique <- unique(names)
match <- match(unique,names)
anno <- anno[match]
return(anno)
}
#' rnb.calculate.fdrp
#'
#' This function calculates the FDRP scores for a given RnBSet and a bam file containing the reads.
#'
#' @param rnb.set RnBSet object containing the CpG sites for which calculation should be conducted
#' @param bam.path path to the bam file containing the reads
#' @param log.path path to the log directory, if not existing it is created
#' @param cores cores available for the analysis
#' @param use.sex.chromosomes Flag indicating if scores are also to be computed for the sex chromosomes
#' @param ... Further parameters passed to \link{calculate.fdrp}
#'
#' @return FDRP scores for the CpG sites in the RnBSet with a higher coverage than COVERAGE.THRESHOLD
#'
#' @author Michael Scherer
#' @examples
#' \dontrun{
#' example.rnb.set <- system.file(file.path("extData","small_rnbSet.zip"),package="WSH")
#' example.bam <- system.file(file.path("extData","small_example.bam"),package="WSH")
#' fdrp <- rnb.calculate.fdrp(rnb.set=example.rnb.set,bam.path=example.bam)
#' }
#' @export
rnb.calculate.fdrp <- function(rnb.set,
bam.path,
log.path=getwd(),
cores=1,
use.sex.chromosomes = FALSE,
...){
logger.start("Computing FDRP from RnBSet object")
anno <- create.annotation(rnb.set,use.sex.chromosomes = use.sex.chromosomes)
fdrps <- calculate.fdrp(bam.path,anno,log.path=log.path,cores=cores,use.sex.chromosomes = use.sex.chromosomes,...)
logger.completed()
return(fdrps)
}
#' rnb.calculate.qfdrp
#'
#' This function calculates the qFDRP scores for a given RnBSet and a bam file containing the reads.
#'
#' @param rnb.set RnBSet object containing the CpG sites for which calculation should be conducted
#' @param bam.path path to the bam file containing the reads
#' @param log.path path to the log directory, if not existing it is created
#' @param cores cores available for the analysis
#' @param use.sex.chromosomes Flag indicating if scores are also to be computed for the sex chromosomes
#' @param ... Further parameters passed to \link{calculate.qfdrp}
#'
#' @return qFDRP scores for the CpG sites in the RnBSet with a higher coverage than COVERAGE.THRESHOLD
#'
#' @author Michael Scherer
#' @examples
#' \dontrun{
#' example.rnb.set <- system.file(file.path("extData","small_rnbSet.zip"),package="WSH")
#' example.bam <- system.file(file.path("extData","small_example.bam"),package="WSH")
#' fdrp <- rnb.calculate.qfdrp(rnb.set=example.rnb.set,bam.path=example.bam)
#' }
#' @export
rnb.calculate.qfdrp <- function(rnb.set,
bam.path,
log.path=getwd(),
cores=1,
use.sex.chromosomes = FALSE,
...){
logger.start("Computing qFDRP from RnBSet object")
anno <- create.annotation(rnb.set,use.sex.chromosomes=use.sex.chromosomes)
qfdrps <- calculate.qfdrp(bam.path,anno,log.path=log.path,cores=cores,use.sex.chromosomes=use.sex.chromosomes,...)
logger.completed()
return(qfdrps)
}
#' rnb.calculate.pdr
#'
#' This function calculates the PDR scores for a given RnBSet and a bam file containing the reads.
#'
#' @param rnb.set RnBSet object containing the CpG sites for which calculation should be conducted
#' @param bam.path path to the bam file containing the reads
#' @param log.path path to the log directory, if not existing it is created
#' @param cores cores available for the analysis
#' @param use.sex.chromosomes Flag indicating if scores are also to be computed for the sex chromosomes
#' @param ... Further parameters passed to \link{calculate.pdr}
#'
#' @return PDR scores for the CpG sites in the RnBSet with a higher coverage than COVERAGE.THRESHOLD
#'
#' @author Michael Scherer
#' @examples
#' \dontrun{
#' example.rnb.set <- system.file(file.path("extData","small_rnbSet.zip"),package="WSH")
#' example.bam <- system.file(file.path("extData","small_example.bam"),package="WSH")
#' fdrp <- rnb.calculate.pdr(rnb.set=example.rnb.set,bam.path=example.bam)
#' }
#' @export
rnb.calculate.pdr <- function(rnb.set,
bam.path,
log.path=getwd(),
cores=1,
use.sex.chromosomes = FALSE,
...){
logger.start("Computing PDR from RnBSet object")
anno <- create.annotation(rnb.set,use.sex.chromosomes=use.sex.chromosomes)
pdrs <- calculate.pdr(bam.path,anno,log.path=log.path,cores=cores,use.sex.chromosomes = use.sex.chromosomes,...)
logger.completed()
return(pdrs)
}
#' rnb.calculate.mhl
#'
#' This function calculates the MHL scores for a given RnBSet and a bam file containing the reads.
#'
#' @param rnb.set RnBSet object containing the CpG sites for which calculation should be conducted
#' @param bam.path path to the bam file containing the reads
#' @param roi.path path to a directory to which output can be added
#'
#' @return MHL scores for the CpG sites in the RnBSet with a higher coverage than COVERAGE.THRESHOLD
#'
#' @export
rnb.calculate.mhl <- function(rnb.set,bam.path,roi.path=getwd()){
logger.start("Computing MHL from RnBSet object")
anno <- create.annotation(rnb.set)
roi.df <- data.frame(Chromosome=seqnames(anno),Start=start(anno),End=end(anno))
roi.location <- file.path(roi.path,"roi.bed")
write.table(roi.df,roi.location,sep="\t",quote=F,row.names=F)
mhl <- calculate.mhl(roi.location,bam.path)
logger.completed()
return(mhl)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.