R/coverage.R

In h5vc: Managing alignment tallies using a hdf5 backend

binnedCoverage <- function( data, sampledata, gccount = FALSE ){
  ret = rowSums( data$Coverages[sampledata$Column,,] )
  names(ret) = sampledata$Sample[order(sampledata$Column)]
  ret = as.data.frame(t(ret))
  tmp <- strsplit( data$h5dapplyInfo$Group, split = "/")[[1]]
  ret$Chrom <- tmp[length(tmp)]
  ret$Start = data$h5dapplyInfo$Blockstart
  ret$End = data$h5dapplyInfo$Blockend
  if( gccount ){
    stopifnot( "Reference" %in% names(data) )
    gcc <- sum( table(data$Reference)[c("1","2")] )
    ret$GCCount <- ifelse( is.na(gcc), 0, gcc )
  }
  ret
}

binnedAFs <- function( data, sampledata, normalise = TRUE, binWidth = 0.05, minCov = 10, minCount = 2 ){
  covs = apply( data$Coverages[sampledata$Column,,], c(1,3), sum)
  counts = apply( data$Counts[,sampledata$Column,,], c(2,4), sum)
  counts[counts < minCount] <- NA #filtering for support
  covs[covs < minCov] <- NA #filtering for power
  afs <- counts / covs
  rownames(afs) = sampledata$Sample[order(sampledata$Column)]
  ret <- apply(afs, 1, function(x) hist(x, breaks = seq(0,1,binWidth), plot = FALSE)$counts)
  rownames(ret) <- as.character(seq(binWidth,1,binWidth))
  if( normalise ){
    ret <- apply( ret, 2, function(x) x / sum(x))
  }
  return(ret)
}