test_bam_tally <- function() {
   f <- system.file("extdata", "test_data_aln",
                    "test_data_aln.concordant_uniq.bam", package = "gmapR")
   bf <- BamFile(f)
   

   bam_tally(bf, BamTallyParam(genome))

 }
 
 test_bam_tally_C <- function() {
   library(BSgenome.Dmelanogaster.UCSC.dm3)
 
   genome <- GmapGenome(Dmelanogaster, create = TRUE)
 
   library(pasillaBamSubset)
   bam <- untreated3_chr4()
   
   bf <- Rsamtools::BamFile(bam)

   which <- GRanges("chr4", IRanges(1e6, 1.3e6))

   gr <- bam_tally(bf, BamTallyParam(genome, indels=TRUE, which = which))

   gr <- bam_tally(bf, BamTallyParam(genome, variant_strand = 1L))
   gr <- bam_tally(bf, BamTallyParam(genome, which = which, variant_strand = 1L))

   empty <- GRanges("chr2L", IRanges(1e6, 2e6))

   gr <- bam_tally(bf, BamTallyParam(genome, which = empty))

   gr <- bam_tally(bf, BamTallyParam(genome, read_pos_breaks = c(1, 15, 30, 40)))

 
   genome <- GmapGenome("hg19_IGIS21")

   which <- GRanges("1", IRanges(1e6, 2e6))

   bam <- "~/share/data/R1047_LIB6635_SAM636095_L1_NXG2449.analyzed.bam"
   bf <- Rsamtools::BamFile(bam)

   gr <- bam_tally(bf, BamTallyParam(genome, which = which, variant_strand = 1L,

                                     read_pos_breaks = c(0L, 10L, 75L),

                                     indels = TRUE))

 }