@@ -13,8 +13,8 @@ Description: GSNAP and GMAP are a pair of tools to align short-read

 Version: 1.3.2

 Depends: R (>= 2.15.0), methods, GenomicRanges

 Imports: IRanges, Rsamtools (>= 1.7.4), rtracklayer (>= 1.17.15), GenomicRanges,

-         GenomicFeatures, Biostrings, VariantAnnotation (>= 1.7.17), tools,

-         Biobase, BSgenome

+         GenomicFeatures, Biostrings, VariantAnnotation (>= 1.7.34), tools,

+         Biobase, BSgenome, methods

 Suggests: RUnit, BSgenome.Dmelanogaster.UCSC.dm3,

 	  BSgenome.Scerevisiae.UCSC.sacCer3, VariantAnnotation,

           org.Hs.eg.db, TxDb.Hsapiens.UCSC.hg19.knownGene,

@@ -6,6 +6,7 @@ importFrom(tools, file_path_as_absolute, file_ext, file_path_sans_ext,

            list_files_with_exts)

 importFrom(Biobase, createPackage)

 import(IRanges)

+import(methods)

 import(GenomicRanges)

 importFrom(Biostrings, getSeq, readDNAStringSet, DNAStringSet)

 importFrom(GenomicRanges, genome, seqinfo)

@@ -16,15 +17,16 @@ importFrom(GenomicFeatures, transcripts, exons)

 importClassesFrom(rtracklayer, RTLFile, FastaFile)

 importFrom(rtracklayer, "referenceSequence<-", import, export, FastaFile)

 importMethodsFrom(rtracklayer, export)

-importFrom(VariantAnnotation, readVcf, ScanVcfParam, fixed)

-importClassesFrom(VariantAnnotation, "VCF")

+importFrom(VariantAnnotation, readVcf, ScanVcfParam, fixed, VRanges, ref, alt,

+           "ref<-", "alt<-")

+importClassesFrom(VariantAnnotation, "VCF", "VRanges")

 importFrom(BSgenome, getSeq, providerVersion)

 ## public API

 export(bam_tally, GmapGenome, GmapGenomeDirectory, GsnapParam,

        directory, BamTallyParam, makeGmapGenomePackage, GsnapOutput,

-       GmapSnps, GmapSnpDirectory, TP53Genome, TP53Which)

+       GmapSnps, GmapSnpDirectory, TP53Genome, TP53Which, summarizeVariants)

 exportClasses(GmapGenome, GmapGenomeDirectory, GmapSnpDirectory,

               GsnapOutput, GmapSnps, BamTallyParam, GsnapParam)

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

 setClass("BamTallyParam",

          representation(genome = "GmapGenome",

                         which = "RangesList",

-                        cycle_breaks = "integerORNULL",

-                        high_base_quality = "integer",

                         minimum_mapq = "integer",

                         concordant_only = "logical",

                         unique_only = "logical",

@@ -24,17 +22,18 @@ setClass("BamTallyParam",

 ### Constructor

 ###

-normArgWhich <- function(x) {

+normArgWhich <- function(x, genome) {

   if (is(x, "GenomicRanges"))

     x <- split(ranges(x), seqnames(x))

   else if (!is(x, "RangesList"))

     stop("'which' must be a GenomicRanges or RangesList")

+  si <- seqinfo(genome)

+  seqinfo(x, new2old = match(seqlevels(si), seqlevels(x))) <-

+    merge(si, seqinfo(x))

   x

 }

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

-                          cycle_breaks = NULL,

-                          high_base_quality = 0L,

                           minimum_mapq = 0L,

                           concordant_only = FALSE, unique_only = FALSE,

                           primary_only = FALSE, ignore_duplicates = FALSE,

@@ -42,12 +41,29 @@ BamTallyParam <- function(genome, which = RangesList(),

                           ignore_query_Ns = FALSE,

                           indels = FALSE)

 {

+  if (!isSingleNumber(minimum_mapq) || minimum_mapq < 0)

+    stop("minimum_mapq must be a single, non-negative, non-NA number")

+  if (!isTRUEorFALSE(concordant_only))

+    stop("concordant_only must be TRUE or FALSE")

+  if (!isTRUEorFALSE(unique_only))

+    stop("unique_only must be TRUE or FALSE")

+  if (!isTRUEorFALSE(primary_only))

+    stop("primary_only must be TRUE or FALSE")

+  if (!isTRUEorFALSE(ignore_duplicates))

+    stop("ignore_duplicates must be TRUE or FALSE")

+  if (!isSingleNumber(min_depth) || minimum_mapq < 0)

+    stop("min_depth must be a single, non-negative, non-NA number")

+  if (!variant_strand %in% c(0, 1, 2))

+    stop("variant_strand must be one of 0, 1, or 2")

+  if (!isTRUEorFALSE(ignore_query_Ns))

+    stop("ignore_query_Ns must be TRUE or FALSE")

+  if (!isTRUEorFALSE(indels))

+    stop("indels must be TRUE or FALSE")

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

   params <- mget(args, environment())

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

-  params$which <- normArgWhich(which)

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

-                      "variant_strand")

+  params$which <- normArgWhich(which, params$genome)

+  integer_params <- c("minimum_mapq", "min_depth", "variant_strand")

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

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

 }

@@ -3,6 +3,22 @@

 ### -------------------------------------------------------------------------

 ###

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

+### Raw tally result

+###

+

+setClass("TallyIIT", representation(ptr = "externalptr", genome = "GmapGenome"))

+

+TallyIIT <- function(ptr, genome) {

+  new("TallyIIT", ptr = ptr, genome = genome)

+}

+

+setMethod("genome", "TallyIIT", function(x) x@genome)

+

+setMethod("show", "TallyIIT", function(object) {

+  cat("Tally IIT object for genome '", genome(genome(object)), "'\n", sep = "")

+})

+

 ## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### High-level wrapper

 ###

@@ -45,40 +61,71 @@ setMethod("bam_tally", "GmapBamReader",

             }

             param_list$genome <- NULL

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

-            tally_names <- c("seqnames", "pos", "ref", "alt", "ncycles",

-                             "ncycles.ref", "count", "count.ref", "count.total",

-                             "high.quality", "high.quality.ref",

-                             "high.quality.total", "mean.quality",

-                             "mean.quality.ref",

-                             "count.pos", "count.pos.ref",

-                             "count.neg", "count.neg.ref",

-                             "read.pos.mean", "read.pos.mean.ref",

-                             "read.pos.var", "read.pos.var.ref")

-            cycle_breaks <- param_list$cycle_breaks

-            if (!is.null(cycle_breaks)) {

-              cycle_breaks <- as.integer(cycle_breaks)

-              break_names <- paste("cycleCount", head(cycle_breaks, -1),

-                                    tail(cycle_breaks, -1), sep = ".")

-              tally_names <- c(tally_names, break_names)

-            }

-            names(tally) <- tally_names

-            if (param_list$variant_strand > 0) {

-              variant_rows <- !is.na(tally$alt)

-              tally <- lapply(tally, `[`, variant_rows)

-            }

-            gr <- GRanges(tally$seqnames,

-                          IRanges(tally$pos,

-                                  width = rep.int(1L, length(tally$pos))),

-                          strand = Rle("+", length(tally$pos)),

-                          location = paste(tally$seqnames, tally$pos,

-                            sep = ":"),

-                          DataFrame(tail(tally, -2)),

-                          seqlengths = seqlengths(genome))

-            seqinfo(gr) <- seqinfo(genome)

-            gr

+            TallyIIT(do.call(.bam_tally_C, c(list(x), param_list)), genome)

           })

+summarizeVariants <- function(x, read_pos_breaks = NULL, high_base_quality = 0L)

+{

+  tally <- .Call(R_tally_iit_parse, x@ptr,

+                 read_pos_breaks,

+                 normArgSingleInteger(high_base_quality),

+                 NULL)

+

+  tally_names <- c("seqnames", "pos", "ref", "alt", "n.read.pos",

+                   "n.read.pos.ref", "raw.count", "raw.count.ref",

+                   "raw.count.total",

+                   "high.quality", "high.quality.ref",

+                   "high.quality.total", "mean.quality",

+                   "mean.quality.ref",

+                   "count.pos", "count.pos.ref",

+                   "count.neg", "count.neg.ref",

+                   "read.pos.mean", "read.pos.mean.ref",

+                   "read.pos.var", "read.pos.var.ref")

+  if (!is.null(read_pos_breaks)) {

+    read_pos_breaks <- as.integer(read_pos_breaks)

+    break_names <- paste("readPosCount", head(read_pos_breaks, -1),

+                         tail(read_pos_breaks, -1), sep = ".")

+    tally_names <- c(tally_names, break_names)

+  }

+  names(tally) <- tally_names

+

+  variant_rows <- !is.na(tally$alt)

+  if (!all(variant_rows))

+    tally <- lapply(tally, `[`, variant_rows)

+  

+  meta_names <- setdiff(tally_names,

+                        c("seqnames", "pos", "ref", "alt", "high.quality",

+                          "high.quality.total"))

+  genome <- genome(x)

+  indel <- nchar(tally$ref) == 0L | nchar(tally$alt) == 0L

+  gr <- with(tally,

+             VRanges(seqnames,

+                     IRanges(pos,

+                             width = ifelse(nchar(alt) == 0, nchar(ref), 1L)),

+                     ref, alt,

+                     ifelse(indel, raw.count.total, high.quality.total),

+                     ifelse(indel, raw.count.ref, high.quality.ref),

+                     ifelse(indel, raw.count, high.quality),

+                     DataFrame(tally[meta_names]),

+                     seqlengths = seqlengths(genome)))

+  seqinfo(gr) <- seqinfo(genome)

+  gr <- normalizeIndelAlleles(gr, genome)

+  gr

+}

+

+normalizeIndelAlleles <- function(x, genome) {

+  is.indel <- nchar(ref(x)) == 0L | nchar(alt(x)) == 0L

+  if (any(is.indel)) {

+    indels <- x[is.indel]

+    indels <- shift(indels, -1)

+    anchor <- getSeq(genome, indels)

+    ref(indels) <- paste0(anchor, ref(indels))

+    alt(indels) <- paste0(anchor, alt(indels))

+    x[is.indel] <- indels

+  }

+  x

+}

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Low-level wrappers

 ###

@@ -98,7 +145,7 @@ normArgTRUEorFALSE <- function(x) {

 }

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

-                         which = NULL, cycle_breaks = NULL,

+                         which = NULL, read_pos_breaks = NULL,

                          high_base_quality = 0L, alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

                          maximum_nhits = 1000000L,

@@ -123,32 +170,28 @@ normArgTRUEorFALSE <- function(x) {

     stop("'genome_dir' must be NULL or a single, non-NA string")

   if (!is.null(db) && !IRanges:::isSingleString(db))

     stop("'db' must be NULL or a single, non-NA string")

-  if (!is.null(cycle_breaks)) {

-    cycle_breaks <- as.integer(cycle_breaks)

-    if (any(is.na(cycle_breaks)))

-      stop("'cycle_breaks' should not contain missing values")

-    if (length(cycle_breaks) < 2)

-      stop("'cycle_breaks' needs at least two elements to define a bin")

-    if (is.unsorted(cycle_breaks))

-      stop("'cycle_breaks' must be sorted")

+  if (!is.null(read_pos_breaks)) {

+    read_pos_breaks <- as.integer(read_pos_breaks)

+    if (any(is.na(read_pos_breaks)))

+      stop("'read_pos_breaks' should not contain missing values")

+    if (length(read_pos_breaks) < 2)

+      stop("'read_pos_breaks' needs at least two elements to define a bin")

+    if (is.unsorted(read_pos_breaks))

+      stop("'read_pos_breaks' must be sorted")

   }

-  iit <- .Call(R_Bamtally_iit, bamreader@.extptr, genome_dir, db, which,

-               normArgSingleInteger(alloclength),

-               normArgSingleInteger(minimum_mapq),

-               normArgSingleInteger(good_unique_mapq),

-               normArgSingleInteger(maximum_nhits),

-               normArgTRUEorFALSE(concordant_only),

-               normArgTRUEorFALSE(unique_only),

-               normArgTRUEorFALSE(primary_only),

-               normArgTRUEorFALSE(ignore_duplicates),

-               normArgSingleInteger(min_depth),

-               normArgSingleInteger(variant_strand),

-               normArgTRUEorFALSE(ignore_query_Ns),

-               normArgTRUEorFALSE(indels),

-               normArgSingleInteger(blocksize),

-               normArgTRUEorFALSE(verbosep))

-  .Call(R_tally_iit_parse, iit,

-        cycle_breaks,

-        normArgSingleInteger(high_base_quality),

-        NULL)

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

+        normArgSingleInteger(alloclength),

+        normArgSingleInteger(minimum_mapq),

+        normArgSingleInteger(good_unique_mapq),

+        normArgSingleInteger(maximum_nhits),

+        normArgTRUEorFALSE(concordant_only),

+        normArgTRUEorFALSE(unique_only),

+        normArgTRUEorFALSE(primary_only),

+        normArgTRUEorFALSE(ignore_duplicates),

+        normArgSingleInteger(min_depth),

+        normArgSingleInteger(variant_strand),

+        normArgTRUEorFALSE(ignore_query_Ns),

+        normArgTRUEorFALSE(indels),

+        normArgSingleInteger(blocksize),

+        normArgTRUEorFALSE(verbosep))

 }

@@ -22,14 +22,14 @@ test_bam_tally_C <- function() {

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

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

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

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

   genome <- GmapGenome("hg19_IGIS21")

   which <- RangesList("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,

-                                    cycle_breaks = c(0L, 10L, 75L),

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

                                     indels = TRUE))

 }

@@ -1,7 +1,19 @@

+library(gmapR)

+

+test_bam_tally <- function() {

+  param <- BamTallyParam(TP53Genome(), TP53Which(), indels = TRUE)

+  bam <- system.file("extdata/H1993.analyzed.bam", 

+                     package="LungCancerLines", mustWork=TRUE)

+  tallies <- bam_tally(bam, param)

+  variants <- summarizeVariants(tallies, NULL, 0L)

+}

+

 test_IITs_not_created <- function() {

   gmapGenome <- GmapGenome(genome="NoGenome", directory = getwd())

   btp <- BamTallyParam(genome=gmapGenome,

                        which = RangesList())

-  bam_file <- system.file("extdata/test_data_aln/test_data_aln.concordant_uniq.bam", package="gmapR", mustWork=TRUE)

+  bam_file <-

+    system.file("extdata/test_data_aln/test_data_aln.concordant_uniq.bam",

+                package="gmapR", mustWork=TRUE)

   checkException(bam_tally(x=bam_file, param=btp))

 }

@@ -14,7 +14,7 @@

 }

 \usage{

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

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

                 high_base_quality = 0L,

                 minimum_mapq = 0L,

                 concordant_only = FALSE, unique_only = FALSE,

@@ -29,8 +29,8 @@

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

     default, the entire genome is processed.

   }

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

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

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

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

     counts are returned.}

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

     read to be counted as high quality.}

@@ -36,9 +36,9 @@

   \item{ref}{The reference base at that position.}

   \item{alt}{The base for the alternate allele, \code{NA} for the

     reference allele row.}

-  \item{ncycles}{The number of unique cycles at which the alternate allele was

+  \item{n.read.pos}{The number of unique cycles at which the alternate allele was

     observed, \code{NA} for the reference allele row.}

-  \item{ncycles.ref}{The number of unique cycles at which the reference

+  \item{n.read.pos.ref}{The number of unique cycles at which the reference

     allele was observed.}

   \item{count}{The number of reads with the alternate allele,

     \code{NA} for the reference allele row.}

@@ -64,7 +64,7 @@

     allele.}

   An additional column is present for each bin formed by

-  the \code{cycle_breaks} parameter, with the read count for that bin.

+  the \code{read_pos_breaks} parameter, with the read count for that bin.

 }

 \author{Michael Lawrence}

@@ -277,7 +277,7 @@ aligned to the TP53 region of the human genome. We'll use this file to

 demonstrate the use of \software{bam\_tally} through the

 \software{gmapR} package. The resulting data structure will contain

 the needed information such as number of alternative alleles, quality

-scores, and nucleotide cycle for each allele.

+scores, and read position for each allele.

 <<run_bamtally, eval=FALSE>>=

 library(gmapR)

@@ -288,7 +288,7 @@ bam_file <- system.file("extdata/H1993.analyzed.bam",

 breaks <- c(0L, 15L, 60L, 75L)

 bqual <- 56L

 mapq <- 13L

-param <- BamTallyParam(genome, cycle_breaks = breaks,

+param <- BamTallyParam(genome, read_pos_breaks = breaks,

                        high_base_quality = bqual,

                        minimum_mapq = mapq,

                        concordant_only = FALSE, unique_only = FALSE,