From 0cca47e4bb0590d2d8870011780e0d16eaba6389 Mon Sep 17 00:00:00 2001 From: mlweilert Date: Thu, 7 Nov 2019 09:20:49 -0600 Subject: [PATCH] Fix: Filter bad chroms when making ExoData object. --- R/AllClasses.R | 63 ++++++++++++++++++++++++++++---------------------- 1 file changed, 35 insertions(+), 28 deletions(-) diff --git a/R/AllClasses.R b/R/AllClasses.R index b41efd1..06e3b69 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -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 @@ -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 @@ -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 @@ -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 @@ -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.") } @@ -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") @@ -103,12 +103,12 @@ 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{ @@ -116,43 +116,50 @@ ExoData <- function(file = NULL, reads = NULL , height = 1 , 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 @@ -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!")