R/prepareForHDF5.R

In h5vc: Managing alignment tallies using a hdf5 backend

prepareForHDF5 <- function( counts, reference ){
  stopifnot(class(reference) == "DNAString")
  ret <- list()
  ret$Counts <- counts[-c(grep("[-+]", dimnames(counts)[[1]])),,,,drop=FALSE]
  ret$Coverages <- apply( ret$Counts, c(2,3,4), sum)
  ret$Deletions <- apply( counts[c("-.front","-.mid","-.back"),,,,drop=FALSE], c(2,3,4), sum )
  ret$Insertions <- apply( counts[c("+.front","+.mid","+.back"),,,,drop=FALSE], c(2,3,4), sum )
  ref <- encodeDNAString( reference ) + 1
  for( base in 1:4){
    idxs = which( ref == base )
    for( strand in 1:2){
      ret$Counts[base,, strand,idxs] <- 0
      ret$Counts[base + 4,, strand,idxs] <- 0
      ret$Counts[base + 8,, strand,idxs] <- 0 
    }
  }
  ret$Reference <- ref - 1 #the tally file encodes the reference as A=0,C=1,G=2,T=3
  return(ret)
}