@@ -6,9 +6,10 @@

 ###

 setClass("BamTallyParam",

-         representation(which = "RangesList",

+         representation(genome = "GmapGenome",

+                        which = "RangesList",

                         cycle_breaks = "integerORNULL",

-                        high_quality_cutoff = "integer",

+                        high_base_quality = "integer",

                         minimum_mapq = "integer",

                         concordant_only = "logical",

                         unique_only = "logical",

@@ -22,8 +23,9 @@ setClass("BamTallyParam",

 ### Constructor

 ###

-BamTallyParam <- function(which = RangesList(), cycle_breaks = NULL,

-                          high_quality_cutoff = 0L,

+BamTallyParam <- function(genome, which = RangesList(),

+                          cycle_breaks = NULL,

+                          high_base_quality = 0L,

                           minimum_mapq = 0L,

                           concordant_only = FALSE, unique_only = FALSE,

                           primary_only = FALSE,

@@ -33,8 +35,9 @@ BamTallyParam <- function(which = RangesList(), cycle_breaks = NULL,

 {

   args <- names(formals(sys.function()))

   params <- mget(args, environment())

+  params$genome <- as(genome, "GmapGenome")

   params$which <- as(which, "RangesList")

-  integer_params <- c("high_quality_cutoff", "minimum_mapq", "min_depth",

+  integer_params <- c("high_base_quality", "minimum_mapq", "min_depth",

                       "variant_strand")

   params[integer_params] <- lapply(params[integer_params], as.integer)

   do.call(new, c("BamTallyParam", params))  

@@ -49,3 +52,17 @@ setAs("BamTallyParam", "list", function(from) {

 })

 setMethod("as.list", "BamTallyParam", function(x) as(x, "list"))

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Utilities

+###

+

+flankingCycleBreaks <- function(read_length, width = 10L) {

+  if (is.na(read_length))

+    return(NULL)

+  if (read_length < 1)

+    stop("'read_length' must be >= 1 or NA")

+  if (width < 0)

+    stop("'width' must be non-negative")

+  c(0L, width, read_length - width, read_length)

+}

@@ -8,30 +8,31 @@

 ###

 setGeneric("bam_tally",

-           function(x, genome, param = BamTallyParam(), ...)

+           function(x, param, ...)

            standardGeneric("bam_tally"),

            signature = "x")

 setMethod("bam_tally", "BamFile",

-          function(x, genome, param, ...)

+          function(x, param, ...)

           {

             x <- GmapBamReader(x)

             callGeneric()

           })

 setMethod("bam_tally", "character",

-          function(x, genome, param, ...)

+          function(x, param, ...)

           {

             x <- BamFile(x)

             callGeneric()

           })

 setMethod("bam_tally", "GmapBamReader",

-          function(x, genome, param, ...)

+          function(x, param, ...)

           {

             param_list <- as.list(param)

             args <- list(...)

             param_list[names(args)] <- args

+            genome <- param_list$genome

             param_list$db <- genome(genome)

             param_list$genome_dir <- path(directory(genome))

             tally <- do.call(.bam_tally_C, c(list(x), param_list))

@@ -85,7 +86,7 @@ normArgTRUEorFALSE <- function(x) {

 .bam_tally_C <- function(bamreader, genome_dir = NULL, db = NULL,

                          which = NULL, cycle_breaks = NULL,

-                         high_quality_cutoff = 0L, alloclength = 200000L,

+                         high_base_quality = 0L, alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

                          maximum_nhits = 1000000L,

                          concordant_only = FALSE, unique_only = FALSE,

@@ -120,7 +121,7 @@ normArgTRUEorFALSE <- function(x) {

   }

   .Call(R_Bamtally_iit, bamreader@.extptr, genome_dir, db, which,

         cycle_breaks,

-        normArgSingleInteger(high_quality_cutoff),

+        normArgSingleInteger(high_base_quality),

         normArgSingleInteger(alloclength),

         normArgSingleInteger(minimum_mapq),

         normArgSingleInteger(good_unique_mapq),

@@ -3,12 +3,10 @@ test_bam_tally <- function() {

                    "test_data_aln.concordant_uniq.bam", package = "gmapR")

   bf <- BamFile(f)

-  bam_tally(bf, genome)

+  bam_tally(bf, BamTallyParam(genome))

 }

 test_bam_tally_C <- function() {

-### FIXME: switch over to using Jeremiah's BAM here, once we get the

-### human gmap index in a package.

   library(BSgenome.Dmelanogaster.UCSC.dm3)

   genome <- GmapGenome(Dmelanogaster, create = TRUE)

@@ -18,11 +16,11 @@ test_bam_tally_C <- function() {

   bf <- Rsamtools::BamFile(bam)

   which <- RangesList(chr4 = IRanges(1e6, 2e6))

-  gr <- bam_tally(bf, genome)

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

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

+  gr <- bam_tally(bf, BamTallyParam(genome))

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

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

   empty <- RangesList(chr2L = IRanges(1e6, 2e6))

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

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

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

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

 }

@@ -14,8 +14,8 @@

 }

 \usage{

-  BamTallyParam(which = RangesList(), cycle_breaks = NULL,

-                high_quality_cutoff = 0L,

+  BamTallyParam(genome, which = RangesList(), cycle_breaks = NULL,

+                high_base_quality = 0L,

                 minimum_mapq = 0L,

                 concordant_only = FALSE, unique_only = FALSE,

                 primary_only = FALSE,

@@ -24,6 +24,7 @@

                 indels = FALSE)

 }

 \arguments{

+  \item{genome}{A \code{GmapGenome} object, or something coercible to one.}

   \item{which}{A \code{RangesList} or something coercible to

     one that limits the tally to that range or set of ranges. By

     default, the entire genome is processed.

@@ -31,7 +32,7 @@

   \item{cycle_breaks}{The breaks, like those passed to \code{\link{cut}}

     for aggregating the per-cycle counts. If \code{NULL}, no per-cycle

     counts are returned.}

-  \item{high_quality_cutoff}{The minimum mapping quality for a

+  \item{high_base_quality}{The minimum mapping quality for a

     read to be counted as high quality.}

   \item{minimum_mapq}{Minimum mapping quality for a read to be counted

     at all.}

@@ -6,8 +6,8 @@

 \alias{bam_tally,GmapBamReader-method}

 \alias{bam_tally}

-\title{Create a Summarization of the Read Sequences and Qualities for Each

-  Genomic Position}

+\title{Per-position Alignment Summaries}

+

 \description{

   Given a set of alignments, for each position in the genome output

   counts for the reference allele and all alternate alleles. Often used

@@ -15,15 +15,13 @@

 }

 \usage{

-\S4method{bam_tally}{BamFile}(x, genome, param = BamTallyParam(), ...)

-\S4method{bam_tally}{character}(x, genome, param = BamTallyParam(), ...)

+\S4method{bam_tally}{BamFile}(x, param, ...)

+\S4method{bam_tally}{character}(x, param, ...)

 }

 \arguments{

   \item{x}{a \code{BamFile} object or string path to a BAM file to read}

-  \item{genome}{the \code{GmapGenome} object corresponding to the

-    alignments in the BAM file}

-  \item{param}{The \code{\linkS4class{BamTallyParam}} object with extra

+  \item{param}{The \code{\linkS4class{BamTallyParam}} object with

     parameters for the tally operation. }

   \item{...}{Arguments that override settings in \code{param}.}

 }