R/poolCoverage.R
7d29a0be
 #' Pool coverage from multiple samples
346ab895
 #'
7d29a0be
 #' Averages the coverage of a list of samples.
346ab895
 #'
 #'
d05fa238
 #' @param all.data List of normals, read with \code{\link{readCoverageFile}}.
7d29a0be
 #' @param remove.chrs Remove these chromosomes from the pool.
 #' @param w \code{numeric(length(all.data))} vector of weights. If \code{NULL},
 #' weight all samples equally.
 #' @return A \code{data.frame} with the averaged coverage over all normals.
 #' @author Markus Riester
d05fa238
 #' @seealso \code{\link{readCoverageFile}}
7d29a0be
 #' @examples
346ab895
 #'
4a7f4816
 #' normal.coverage.file <- system.file("extdata", "example_normal.txt.gz",
346ab895
 #'     package = "PureCN")
4a7f4816
 #' normal2.coverage.file <- system.file("extdata", "example_normal2.txt.gz",
346ab895
 #'     package = "PureCN")
7d29a0be
 #' normal.coverage.files <- c(normal.coverage.file, normal2.coverage.file)
ec6f11d9
 #' pool <- poolCoverage(lapply(normal.coverage.files, readCoverageFile),
346ab895
 #'      remove.chrs = c("chrX", "chrY"))
 #'
7d29a0be
 #' @export poolCoverage
e769c0fc
 poolCoverage <- function(all.data, remove.chrs=c(), w = NULL) {
44326645
     pool <- all.data[[1]]
 
ea58a451
     if (length(all.data) == 1) {
         return(.removeChr(pool, remove.chrs))
     }
8a972c87
     if (is.null(w)) {
346ab895
         w <- rep(1, length(all.data))
8a972c87
     } else if (length(w) != length(all.data)) {
         .stopUserError("all.data and w have different lengths.")
     }
 
     pool$coverage <- 0
     pool$counts <- 0
ea58a451
 
8a972c87
     for (i in seq_along(all.data)) {
ea58a451
         pool$coverage <- pool$coverage + (w[i] * all.data[[i]]$coverage)
8a972c87
         pool$counts <- pool$counts + (w[i] * all.data[[i]]$counts)
ea58a451
     }
8a972c87
     pool <- .addAverageCoverage(pool)
ea58a451
     return(.removeChr(pool, remove.chrs))
7d29a0be
 }
ea58a451
 
346ab895
 .removeChr <- function(pool, remove.chrs = c()) {
44326645
     idx <- seqnames(pool) %in% remove.chrs
     if (sum(idx)) {
         pool[idx]$coverage <- NA
         pool[idx]$average.coverage <- NA
     }
ea58a451
     pool
 }