Skip to content

Fix: Filter bad chroms when making ExoData object. #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mlweilert wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fix: Filter bad chroms when making ExoData object. #9

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 35 additions & 28 deletions R/AllClasses.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
##' @importFrom biovizBase flatGrl
##' @importFrom Rsamtools ScanBamParam scanBamFlag
##' @importFrom GenomicAlignments readGAlignments
##' @importFrom GenomicAlignments readGAlignments
##' @importFrom GenomicRanges GRanges
##' @importFrom IRanges slice
##' @importFrom data.table ":=" rbindlist .N data.table as.data.table
Expand Down Expand Up @@ -38,7 +38,7 @@ setValidity("ExoData",

##' ExoData object and constructors
##'
##' \code{ExoData} is a subclass of \code{GenomicRanges}, used to asses the
##' \code{ExoData} is a subclass of \code{GenomicRanges}, used to asses the
##' quality of ChIP-exo/nexus sample.
##'
##' @param file a character value with location of the bam file with the aligned
Expand All @@ -49,9 +49,9 @@ setValidity("ExoData",
##' of the experiment into a set of regions.
##' @param mc.cores a numeric value with the number of cores to use,
##' i.e. at most how many child processes will be run simultaneously.
##' @param nregions a numeric value indicating the number of regions sampled to
##' @param nregions a numeric value indicating the number of regions sampled to
##' estimate the quality parameter distributions. The default value is 1e3.
##' @param ntimes a numeric value indicating the number of times that regions are
##' @param ntimes a numeric value indicating the number of times that regions are
##' sampled to estimate the quality parameter distributions. The default value
##' is 1e2.
##' @param save.reads a logical value to indicate if the reads are stored in the
Expand All @@ -67,11 +67,11 @@ setValidity("ExoData",
##' @docType class
##'
##' @examples
##'
##'
##' files <- list.files(system.file("extdata",package = "ChIPexoQualExample"),
##' full.names = TRUE)
##' ExoData(files[5],mc.cores = 2L)
##'
##'
##'
##' @rdname ExoData
##' @export
Expand All @@ -80,7 +80,7 @@ ExoData <- function(file = NULL, reads = NULL , height = 1 ,
save.reads = FALSE,nregions = 1e3,ntimes = 1e2,
verbose = TRUE)
{

if(!is.null(file) & !is.null(reads)){
stop("Both 'file' and 'reads' are available, can't use both.")
}
Expand All @@ -89,12 +89,12 @@ ExoData <- function(file = NULL, reads = NULL , height = 1 ,
}

if(!is.null(file))stopifnot(is.character( file),file.exists(file))

if(!is.null(reads))stopifnot(class(reads) %in% c("GAlignments","GRanges"))

stopifnot(is.numeric(height),height >= 1)
stopifnot(is.logical(save.reads))

if(verbose){
if(is.null(file))message("Using 'reads' argument")
else message("Using 'file' argument")
Expand All @@ -103,56 +103,63 @@ ExoData <- function(file = NULL, reads = NULL , height = 1 ,
if(verbose) message("Creating ExoData object using aligned reads in ",file)
}
if(verbose) message("Keeping reads in object: ",ifelse(save.reads,"Yes","No"))

if(!is.null(file)){
if(verbose) message("Loading experiment reads")
paramFwd <- ScanBamParam(flag = scanBamFlag(isMinusStrand = FALSE))
paramRev <- ScanBamParam(flag = scanBamFlag(isMinusStrand = TRUE))
paramRev <- ScanBamParam(flag = scanBamFlag(isMinusStrand = TRUE))

fwdReads <- readGAlignments(file,param = paramFwd)
revReads <- readGAlignments(file,param = paramRev)
}else{
fwdReads <- subset(reads,as.character(strand(reads)) == "+")
revReads <- subset(reads,as.character(strand(reads)) == "-")
file <- ""
}

if(class(fwdReads) == "GAlignments")fwdReads <- as(fwdReads,"GRanges")
if(class(revReads) == "GAlignments")revReads <- as(revReads,"GRanges")


if(verbose) message("Removing chromosomes with insufficient coverage to build ExoData class")
fwdChr <- unique(seqnames(fwdReads))
revChr <- unique(seqnames(revReads))
bothChr <- intersect(fwdChr, revChr) #ensure there is coverage across both strands per chromosome
fwdReads<-fwdReads[grep(pattern = paste("^", bothChr, "$", collapse = "|", sep = ""), x = seqnames(fwdReads))]
revReads<-revReads[grep(pattern = paste("^", bothChr, "$", collapse = "|", sep = ""), x = seqnames(revReads))]

cover <- coverage(fwdReads) + coverage(revReads)

rlist <- slice(cover,lower = height,rangesOnly = TRUE)

if(any(vapply(rlist,length,0L) == 0)){
chr <- names(which(vapply(rlist,length,0L) > 0))
rlist <- rlist[chr]
}else{
chr <- names(rlist)
}

fwdReads <- split(fwdReads,as.character(seqnames(fwdReads)))
revReads <- split(revReads,as.character(seqnames(revReads)))

fwdReads <- fwdReads[chr]
revReads <- revReads[chr]

if(verbose) message("Calculating summary statistics")

if(Sys.info()[["sysname"]] == "Windows"){
snow <- SnowParam(workers = mc.cores, type = "SOCK")
stats <- bpmapply(calculateSummary,rlist,fwdReads,revReads,
BPPARAM = snow,SIMPLIFY = FALSE)
BPPARAM = snow,SIMPLIFY = FALSE)
}else{
stats <- bpmapply(calculateSummary,rlist,fwdReads,revReads,
BPPARAM = MulticoreParam(workers = mc.cores),
SIMPLIFY = FALSE)
SIMPLIFY = FALSE)
}

regions <- as(rlist,"GRanges")
mcols(regions) <- do.call(rbind,stats)
nreads <- sum(vapply(fwdReads,length,1)) + sum(vapply(revReads, length, 1))

if(save.reads){
fwdReads <- biovizBase::flatGrl(fwdReads)
mcols(fwdReads) <- NULL
Expand All @@ -163,21 +170,21 @@ ExoData <- function(file = NULL, reads = NULL , height = 1 ,
reads <- GRanges()
}
depth <- uniquePos <- width <- NULL

if(verbose) message("Calculating quality scores distribution")

DT <- formatRegions(regions)

paramList <- lapply(seq_len(ntimes),calculateParamDist,DT,nregions)

param <- rbindlist(paramList)

estimate <- term <- NULL

rm(paramList,DT)
paramDist <- list("beta1" = param[term == "uniquePos",(estimate)] ,
"beta2" = param[term == "width",(estimate)])

metadata(regions) <- list("file"=file,"nreads"=nreads,
"ntimes"=ntimes,"nregions"=nregions)
if(verbose) message("Done!")
Expand Down