@@ -205,6 +205,13 @@ normArgSingleInteger <- function(x) {

     stop("'", name, "' should be a single, non-NA integer")

   x

 }

+normArgSingleNumber <- function(x) {

+    name <- deparse(substitute(x))

+    x <- as.numeric(x)

+    if (!isSingleNumber(x))

+        stop("'", name, "' should be a single, non-NA number")

+    x

+}

 normArgTRUEorFALSE <- function(x) {

   name <- deparse(substitute(x))

   if (!isTRUEorFALSE(x))

@@ -227,11 +234,11 @@ normArgSingleCharacterOrNULL <- function(x) {

                          maximum_nhits = 1000000L,

                          concordant_only = FALSE, unique_only = FALSE,

                          primary_only = FALSE, ignore_duplicates = FALSE,

-                         min_depth = 0L, variant_strand = 0L,

+                         min_depth = 0L, variant_strand = 0L, variant_pct = 0,

                          ignore_query_Ns = FALSE,

                          indels = FALSE,

                          blocksize = 1000L, verbosep = FALSE,

-                         include_soft_clips = 0L,

+                         min_softclip = 0L, max_softclip = 0L,

                          exon_iit = NULL, xs = FALSE, read_pos = FALSE,

                          min_base_quality = 0L, noncovered = FALSE, nm = FALSE)

 {

@@ -258,11 +265,13 @@ normArgSingleCharacterOrNULL <- function(x) {

         normArgTRUEorFALSE(ignore_duplicates),

         normArgSingleInteger(min_depth),

         normArgSingleInteger(variant_strand),

+        normArgSingleNumber(variant_pct),

         normArgTRUEorFALSE(ignore_query_Ns),

         normArgTRUEorFALSE(indels),

         normArgSingleInteger(blocksize),

         normArgTRUEorFALSE(verbosep),

-        normArgSingleInteger(include_soft_clips),

+        normArgSingleInteger(min_softclip),

+        normArgSingleInteger(max_softclip),

         normArgSingleCharacterOrNULL(exon_iit),

         normArgTRUEorFALSE(xs),

         normArgTRUEorFALSE(read_pos),

@@ -152,7 +152,8 @@ variantSummary <- function(x, read_pos_breaks = NULL,

   gr <- with(tally,

              VRanges(seqnames,

                      IRanges(pos,

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

+                             width = ifelse(nchar(alt) == 0 & !is.na(alt),

+                                            nchar(ref), 1L)),

                      ref, alt,

                      count.total, count.ref, count,

                      seqlengths = seqlengths(genome)))

@@ -173,7 +174,7 @@ checkTallyConsistency <- function(x) {

 }

 normalizeIndelAlleles <- function(x, genome) {

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

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

   if (any(is.indel)) {

     indels <- x[is.indel]

     flanks <- flank(indels, 1)

@@ -7,13 +7,13 @@

 ### Raw tally result

 ###

-setClass("TallyIIT", representation(ptr = "externalptr",

-                                    genome = "GmapGenome",

+setClass("TallyIIT", representation(genome = "GmapGenome",

                                     bam = "BamFile",

                                     xs = "logical",

                                     read_pos = "logical",

                                     nm = "logical",

-                                    codon = "logical"))

+                                    codon = "logical"),

+         contains = "IIT")

 TallyIIT <- function(ptr, genome, bam, xs, read_pos, nm, codon) {

   new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

@@ -110,7 +110,7 @@ variantSummary <- function(x, read_pos_breaks = NULL,

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

                    "n.read.pos", "n.read.pos.ref",

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

+                   "count", "count.ref", "raw.count.total", "count.total",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

                    "count.del.plus", "count.del.minus",

@@ -280,6 +280,7 @@ variantSummaryColumnDescriptions <-

   desc <- c(

     n.read.pos = "Number of unique read positions for the ALT",

     n.read.pos.ref = "Number of unique read positions for the REF",

+    raw.count.total = "Total depth, before quality filtering",  

     count.plus = "Positive strand ALT count",

     count.plus.ref = "Positive strand REF count",

     count.minus = "Negative strand ALT count",

@@ -12,11 +12,12 @@ setClass("TallyIIT", representation(ptr = "externalptr",

                                     bam = "BamFile",

                                     xs = "logical",

                                     read_pos = "logical",

-                                    nm = "logical"))

+                                    nm = "logical",

+                                    codon = "logical"))

-TallyIIT <- function(ptr, genome, bam, xs, read_pos, nm) {

+TallyIIT <- function(ptr, genome, bam, xs, read_pos, nm, codon) {

   new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

-      read_pos = read_pos, nm = nm)

+      read_pos = read_pos, nm = nm, codon = codon)

 }

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

@@ -73,7 +74,8 @@ setMethod("bam_tally", "GmapBamReader",

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

                      as(x, "BamFile"), xs=param_list$xs,

-                     read_pos=param_list$read_pos, nm=param_list$nm)

+                     read_pos=param_list$read_pos, nm=param_list$nm,

+                     codon=!is.null(param_list$exon_iit))

           })

 variantSummary <- function(x, read_pos_breaks = NULL,

@@ -111,7 +113,7 @@ variantSummary <- function(x, read_pos_breaks = NULL,

                    "count", "count.ref", "count.total",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

-                   "del.count.plus", "del.count.minus",

+                   "count.del.plus", "count.del.minus",

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

                    "read.pos.var", "read.pos.var.ref",

                    "mdfne", "mdfne.ref", "codon.strand",

@@ -119,8 +121,6 @@ variantSummary <- function(x, read_pos_breaks = NULL,

                    "count.xs.minus", "count.xs.minus.ref",

                    "count.high.nm", "count.high.nm.ref")

-  tally$codon.strand <- structure(tally$codon.strand, class="factor",

-                                  levels=strand())

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -129,7 +129,10 @@ variantSummary <- function(x, read_pos_breaks = NULL,

     tally_names <- c(tally_names, break_names)

   }

   names(tally) <- tally_names

-  tally <- Filter(Negate(is.null), tally)  

+  tally$codon.strand <- if (x@codon) {

+      structure(tally$codon.strand, class="factor", levels=levels(strand()))

+  }

+  tally <- Filter(Negate(is.null), tally)

   if (!keep_ref_rows) {

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

@@ -143,7 +146,8 @@ variantSummary <- function(x, read_pos_breaks = NULL,

   genome <- genome(x)

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

   metacols <- DataFrame(tally[meta_names])

-  mcols(metacols) <- variantSummaryColumnDescriptions(read_pos_breaks)

+  mcols(metacols) <- variantSummaryColumnDescriptions(read_pos_breaks, x@xs,

+                                                      high_nm_score, x@codon)

   gr <- with(tally,

              VRanges(seqnames,

@@ -270,32 +274,37 @@ normArgSingleCharacterOrNULL <- function(x) {

 ### Column metadata

 ###

-variantSummaryColumnDescriptions <- function(read_pos_breaks) {

+variantSummaryColumnDescriptions <-

+    function(read_pos_breaks, xs, high_nm, codon)

+{

   desc <- c(

     n.read.pos = "Number of unique read positions for the ALT",

     n.read.pos.ref = "Number of unique read positions for the REF",

-    raw.count = "Raw ALT count",

-    raw.count.ref = "Raw REF count",

-    raw.count.total = "Raw total count",

-    mean.quality = "Average ALT base quality",

-    mean.quality.ref = "Average REF base quality",

-    count.plus = "Raw positive strand ALT count",

-    count.plus.ref = "Raw positive strand REF count",

-    count.minus = "Raw negative strand ALT count",

-    count.minus.ref = "Raw negative strand REF count",

+    count.plus = "Positive strand ALT count",

+    count.plus.ref = "Positive strand REF count",

+    count.minus = "Negative strand ALT count",

+    count.minus.ref = "Negative strand REF count",

+    count.del.plus = "Plus strand deletion count",

+    count.del.minus = "Count strand deletion count",

     read.pos.mean = "Average read position for the ALT",

     read.pos.mean.ref = "Average read position for the ALT",

     read.pos.var = "Variance in read position for the ALT",

     read.pos.var.ref = "Variance in read position for the REF",

     mdfne = "Median distance from nearest end of read for the ALT",

     mdfne.ref = "Median distance from nearest end of read for the REF",

-    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)",

-    del.count.plus = "plus strand deletion count",

-    del.count.minus = "minus strand deletion count",

-    count.xs.plus = "Plus strand XS counts",

-    count.xs.plus.ref = "Plus strand reference XS counts",

-    count.xs.minus = "Minus strand XS counts",

-    count.xs.minus.ref = "Minus strand reference XS counts")

+      if (codon) {

+          c(codon.strand = "Strand of transcription for the codon")

+      },

+      if (xs) {

+          c(count.xs.plus = "Plus strand XS counts",

+            count.xs.plus.ref = "Plus strand reference XS counts",

+            count.xs.minus = "Minus strand XS counts",

+            count.xs.minus.ref = "Minus strand reference XS counts")

+      },

+      count.high.nm = paste("Count of reads with >=", high_nm,

+          "NM score for the ALT"),

+      count.high.nm.ref = paste("Count of reads with >=", high_nm,

+          "NM score for the REF"))

   if (length(read_pos_breaks) > 0L) {

     break_desc <- paste0("Raw ALT count in read position range [",

                          head(read_pos_breaks, -1), ",",

@@ -11,11 +11,12 @@ setClass("TallyIIT", representation(ptr = "externalptr",

                                     genome = "GmapGenome",

                                     bam = "BamFile",

                                     xs = "logical",

-                                    read_pos = "logical"))

+                                    read_pos = "logical",

+                                    nm = "logical"))

-TallyIIT <- function(ptr, genome, bam, xs, read_pos) {

+TallyIIT <- function(ptr, genome, bam, xs, read_pos, nm) {

   new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

-      read_pos = read_pos)

+      read_pos = read_pos, nm = nm)

 }

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

@@ -72,16 +73,22 @@ setMethod("bam_tally", "GmapBamReader",

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

                      as(x, "BamFile"), xs=param_list$xs,

-                     read_pos=param_list$read_pos)

+                     read_pos=param_list$read_pos, nm=param_list$nm)

           })

-variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

-                           keep_ref_rows = FALSE, read_length = NA_integer_)

+variantSummary <- function(x, read_pos_breaks = NULL,

+                           keep_ref_rows = FALSE, read_length = NA_integer_,

+                           high_nm_score = NA_integer_)

 {

   read_length <- as.integer(read_length)

   if (length(read_length) != 1L) {

     stop("'read_length' must be a single integer")

   }

+  high_nm_score <- as.integer(high_nm_score)

+  stopifnot(length(high_nm_score) == 1L)

+  if (!is.na(high_nm_score) && !x@nm) {

+      stop("'high_nm_score is not NA but NM scores were not tallied")

+  }

   if (length(read_pos_breaks) > 0L) {

     if (!x@read_pos) {

       stop("'read_pos_breaks' non-empty, but read positions were not tallied")

@@ -97,25 +104,23 @@ variantSummary <- 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, read_length, x@xs)

+                 NULL, read_length, x@xs, high_nm_score)

   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", "count.ref", "count.total",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

                    "del.count.plus", "del.count.minus",

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

                    "read.pos.var", "read.pos.var.ref",

-                   "mdfne", "mdfne.ref", "codon.dir",

+                   "mdfne", "mdfne.ref", "codon.strand",

                    "count.xs.plus", "count.xs.plus.ref",

-                   "count.xs.minus", "count.xs.minus.ref")

-  

+                   "count.xs.minus", "count.xs.minus.ref",

+                   "count.high.nm", "count.high.nm.ref")

+

+  tally$codon.strand <- structure(tally$codon.strand, class="factor",

+                                  levels=strand())

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -124,6 +129,7 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

     tally_names <- c(tally_names, break_names)

   }

   names(tally) <- tally_names

+  tally <- Filter(Negate(is.null), tally)  

   if (!keep_ref_rows) {

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

@@ -131,23 +137,20 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

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

   }

-  meta_names <- setdiff(tally_names,

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

-                          "high.quality.ref", "high.quality.total"))

+  meta_names <- setdiff(names(tally),

+                        c("seqnames", "pos", "ref", "alt", "count",

+                          "count.ref", "count.total"))

   genome <- genome(x)

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

   metacols <- DataFrame(tally[meta_names])

   mcols(metacols) <- variantSummaryColumnDescriptions(read_pos_breaks)

-  samecols = tally[["ref"]] != tally[["alt"]]

   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),

+                     count.total, count.ref, count,

                      seqlengths = seqlengths(genome)))

   mcols(gr) <- metacols

   checkTallyConsistency(gr)

@@ -160,11 +163,8 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

 checkTallyConsistency <- function(x) {

   with(mcols(x), {

-    stopifnot(all(raw.count + raw.count.ref <= raw.count.total, na.rm=TRUE))

-    stopifnot(all(altDepth(x) <= raw.count, na.rm=TRUE))

-    stopifnot(all(refDepth(x) <= raw.count.ref, na.rm=TRUE))

-    stopifnot(all(count.plus + count.minus == raw.count, na.rm=TRUE))

-    stopifnot(all(count.plus.ref + count.minus.ref == raw.count.ref))

+    stopifnot(all(count.plus + count.minus == altDepth(x), na.rm=TRUE))

+    stopifnot(all(count.plus.ref + count.minus.ref == refDepth(x)))

   })

 }

@@ -228,7 +228,7 @@ normArgSingleCharacterOrNULL <- function(x) {

                          blocksize = 1000L, verbosep = FALSE,

                          include_soft_clips = 0L,

                          exon_iit = NULL, xs = FALSE, read_pos = FALSE,

-                         min_base_quality = 0L, noncovered = FALSE)

+                         min_base_quality = 0L, noncovered = FALSE, nm = FALSE)

 {

   if (!is(bamreader, "GmapBamReader"))

     stop("'bamreader' must be a GmapBamReader")

@@ -262,7 +262,8 @@ normArgSingleCharacterOrNULL <- function(x) {

         normArgTRUEorFALSE(xs),

         normArgTRUEorFALSE(read_pos),

         normArgSingleInteger(min_base_quality),

-        normArgTRUEorFALSE(noncovered))

+        normArgTRUEorFALSE(noncovered),

+        normArgTRUEorFALSE(nm))

 }

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

@@ -82,9 +82,9 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

   if (length(read_length) != 1L) {

     stop("'read_length' must be a single integer")

   }

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

+  if (length(read_pos_breaks) > 0L) {

     if (!x@read_pos) {

-      stop("'read_pos_breaks' non-NULL, but read positions were not tallied")

+      stop("'read_pos_breaks' non-empty, but read positions were not tallied")

     }

     read_pos_breaks <- as.integer(read_pos_breaks)

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

@@ -71,7 +71,8 @@ setMethod("bam_tally", "GmapBamReader",

             param_list$genome <- NULL

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

-                     as(x, "BamFile"), param_list[c("xs", "read_pos")])

+                     as(x, "BamFile"), xs=param_list$xs,

+                     read_pos=param_list$read_pos)

           })

 variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

@@ -9,10 +9,13 @@

 setClass("TallyIIT", representation(ptr = "externalptr",

                                     genome = "GmapGenome",

-                                    bam = "BamFile"))

+                                    bam = "BamFile",

+                                    xs = "logical",

+                                    read_pos = "logical"))

-TallyIIT <- function(ptr, genome, bam) {

-  new("TallyIIT", ptr = ptr, genome = genome, bam = bam)

+TallyIIT <- function(ptr, genome, bam, xs, read_pos) {

+  new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

+      read_pos = read_pos)

 }

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

@@ -68,7 +71,7 @@ setMethod("bam_tally", "GmapBamReader",

             param_list$genome <- NULL

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

-                     as(x, "BamFile"))

+                     as(x, "BamFile"), param_list[c("xs", "read_pos")])

           })

 variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

@@ -78,10 +81,23 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

   if (length(read_length) != 1L) {

     stop("'read_length' must be a single integer")

   }

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

+    if (!x@read_pos) {

+      stop("'read_pos_breaks' non-NULL, but read positions were not tallied")

+    }

+    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")

+  }

+

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

                  read_pos_breaks,

                  normArgSingleInteger(high_base_quality),

-                 NULL, read_length)

+                 NULL, read_length, x@xs)

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

                    "n.read.pos", "n.read.pos.ref",

@@ -92,9 +108,13 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

                    "mean.quality.ref",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

+                   "del.count.plus", "del.count.minus",

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

                    "read.pos.var", "read.pos.var.ref",

-                   "mdfne", "mdfne.ref", "codon.dir")

+                   "mdfne", "mdfne.ref", "codon.dir",

+                   "count.xs.plus", "count.xs.plus.ref",

+                   "count.xs.minus", "count.xs.minus.ref")

+  

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -194,7 +214,7 @@ normArgSingleCharacterOrNULL <- function(x) {

 }

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

-                         which = NULL, read_pos_breaks = NULL,

+                         which = NULL,

                          high_base_quality = 0L, desired_read_group = NULL,

                          alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

@@ -206,7 +226,8 @@ normArgSingleCharacterOrNULL <- function(x) {

                          indels = FALSE,

                          blocksize = 1000L, verbosep = FALSE,

                          include_soft_clips = 0L,

-                         exon_iit = NULL)

+                         exon_iit = NULL, xs = FALSE, read_pos = FALSE,

+                         min_base_quality = 0L, noncovered = FALSE)

 {

   if (!is(bamreader, "GmapBamReader"))

     stop("'bamreader' must be a GmapBamReader")

@@ -219,15 +240,6 @@ normArgSingleCharacterOrNULL <- function(x) {

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

   if (!is.null(desired_read_group) && !isSingleString(desired_read_group))

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

-  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")

-  }

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

         desired_read_group,

         normArgSingleInteger(alloclength),

@@ -245,7 +257,11 @@ normArgSingleCharacterOrNULL <- function(x) {

         normArgSingleInteger(blocksize),

         normArgTRUEorFALSE(verbosep),

         normArgSingleInteger(include_soft_clips),

-        normArgSingleCharacterOrNULL(exon_iit))

+        normArgSingleCharacterOrNULL(exon_iit),

+        normArgTRUEorFALSE(xs),

+        normArgTRUEorFALSE(read_pos),

+        normArgSingleInteger(min_base_quality),

+        normArgTRUEorFALSE(noncovered))

 }

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

@@ -271,7 +287,13 @@ variantSummaryColumnDescriptions <- function(read_pos_breaks) {

     read.pos.var.ref = "Variance in read position for the REF",

     mdfne = "Median distance from nearest end of read for the ALT",

     mdfne.ref = "Median distance from nearest end of read for the REF",

-    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)" )

+    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)",

+    del.count.plus = "plus strand deletion count",

+    del.count.minus = "minus strand deletion count",

+    count.xs.plus = "Plus strand XS counts",

+    count.xs.plus.ref = "Plus strand reference XS counts",

+    count.xs.minus = "Minus strand XS counts",

+    count.xs.minus.ref = "Minus strand reference XS counts")

   if (length(read_pos_breaks) > 0L) {

     break_desc <- paste0("Raw ALT count in read position range [",

                          head(read_pos_breaks, -1), ",",

@@ -9,13 +9,10 @@

 setClass("TallyIIT", representation(ptr = "externalptr",

                                     genome = "GmapGenome",

-                                    bam = "BamFile",

-                                    xs = "logical",

-                                    read_pos = "logical"))

+                                    bam = "BamFile"))

-TallyIIT <- function(ptr, genome, bam, xs, read_pos) {

-  new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

-      read_pos = read_pos)

+TallyIIT <- function(ptr, genome, bam) {

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

 }

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

@@ -71,7 +68,7 @@ setMethod("bam_tally", "GmapBamReader",

             param_list$genome <- NULL

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

-                     as(x, "BamFile"), param_list[c("xs", "read_pos")])

+                     as(x, "BamFile"))

           })

 variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

@@ -81,23 +78,10 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

   if (length(read_length) != 1L) {

     stop("'read_length' must be a single integer")

   }

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

-    if (!x@read_pos) {

-      stop("'read_pos_breaks' non-NULL, but read positions were not tallied")

-    }

-    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")

-  }

-

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

                  read_pos_breaks,

                  normArgSingleInteger(high_base_quality),

-                 NULL, read_length, x@xs)

+                 NULL, read_length)

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

                    "n.read.pos", "n.read.pos.ref",

@@ -108,13 +92,9 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

                    "mean.quality.ref",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

-                   "del.count.plus", "del.count.minus",

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

                    "read.pos.var", "read.pos.var.ref",

-                   "mdfne", "mdfne.ref", "codon.dir",

-                   "count.xs.plus", "count.xs.plus.ref",

-                   "count.xs.minus", "count.xs.minus.ref")

-  

+                   "mdfne", "mdfne.ref", "codon.dir")

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -214,7 +194,7 @@ normArgSingleCharacterOrNULL <- function(x) {

 }

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

-                         which = NULL,

+                         which = NULL, read_pos_breaks = NULL,

                          high_base_quality = 0L, desired_read_group = NULL,

                          alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

@@ -226,8 +206,7 @@ normArgSingleCharacterOrNULL <- function(x) {

                          indels = FALSE,

                          blocksize = 1000L, verbosep = FALSE,

                          include_soft_clips = 0L,

-                         exon_iit = NULL, xs = FALSE, read_pos = FALSE,

-                         min_base_quality = 0L, noncovered = FALSE)

+                         exon_iit = NULL)

 {

   if (!is(bamreader, "GmapBamReader"))

     stop("'bamreader' must be a GmapBamReader")

@@ -240,6 +219,15 @@ normArgSingleCharacterOrNULL <- function(x) {

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

   if (!is.null(desired_read_group) && !isSingleString(desired_read_group))

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

+  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")

+  }

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

         desired_read_group,

         normArgSingleInteger(alloclength),

@@ -257,11 +245,7 @@ normArgSingleCharacterOrNULL <- function(x) {

         normArgSingleInteger(blocksize),

         normArgTRUEorFALSE(verbosep),

         normArgSingleInteger(include_soft_clips),

-        normArgSingleCharacterOrNULL(exon_iit),

-        normArgTRUEorFALSE(xs),

-        normArgTRUEorFALSE(read_pos),

-        normArgSingleInteger(min_base_quality),

-        normArgTRUEorFALSE(noncovered))

+        normArgSingleCharacterOrNULL(exon_iit))

 }

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

@@ -287,13 +271,7 @@ variantSummaryColumnDescriptions <- function(read_pos_breaks) {

     read.pos.var.ref = "Variance in read position for the REF",

     mdfne = "Median distance from nearest end of read for the ALT",

     mdfne.ref = "Median distance from nearest end of read for the REF",

-    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)",

-    del.count.plus = "plus strand deletion count",

-    del.count.minus = "minus strand deletion count",

-    count.xs.plus = "Plus strand XS counts",

-    count.xs.plus.ref = "Plus strand reference XS counts",

-    count.xs.minus = "Minus strand XS counts",

-    count.xs.minus.ref = "Minus strand reference XS counts")

+    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)" )

   if (length(read_pos_breaks) > 0L) {

     break_desc <- paste0("Raw ALT count in read position range [",

                          head(read_pos_breaks, -1), ",",

@@ -9,10 +9,13 @@

 setClass("TallyIIT", representation(ptr = "externalptr",

                                     genome = "GmapGenome",

-                                    bam = "BamFile"))

+                                    bam = "BamFile",

+                                    xs = "logical",

+                                    read_pos = "logical"))

-TallyIIT <- function(ptr, genome, bam) {

-  new("TallyIIT", ptr = ptr, genome = genome, bam = bam)

+TallyIIT <- function(ptr, genome, bam, xs, read_pos) {

+  new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs,

+      read_pos = read_pos)

 }

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

@@ -68,7 +71,7 @@ setMethod("bam_tally", "GmapBamReader",

             param_list$genome <- NULL

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

-                     as(x, "BamFile"))

+                     as(x, "BamFile"), param_list[c("xs", "read_pos")])

           })

 variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

@@ -78,10 +81,23 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

   if (length(read_length) != 1L) {

     stop("'read_length' must be a single integer")

   }

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

+    if (!x@read_pos) {

+      stop("'read_pos_breaks' non-NULL, but read positions were not tallied")

+    }

+    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")

+  }

+

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

                  read_pos_breaks,

                  normArgSingleInteger(high_base_quality),

-                 NULL, read_length)

+                 NULL, read_length, x@xs)

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

                    "n.read.pos", "n.read.pos.ref",

@@ -92,9 +108,13 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

                    "mean.quality.ref",

                    "count.plus", "count.plus.ref",

                    "count.minus", "count.minus.ref",

+                   "del.count.plus", "del.count.minus",

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

                    "read.pos.var", "read.pos.var.ref",

-                   "mdfne", "mdfne.ref", "codon.dir")

+                   "mdfne", "mdfne.ref", "codon.dir",

+                   "count.xs.plus", "count.xs.plus.ref",

+                   "count.xs.minus", "count.xs.minus.ref")

+  

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -194,7 +214,7 @@ normArgSingleCharacterOrNULL <- function(x) {

 }

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

-                         which = NULL, read_pos_breaks = NULL,

+                         which = NULL,

                          high_base_quality = 0L, desired_read_group = NULL,

                          alloclength = 200000L,

                          minimum_mapq = 0L, good_unique_mapq = 35L,

@@ -206,7 +226,8 @@ normArgSingleCharacterOrNULL <- function(x) {

                          indels = FALSE,

                          blocksize = 1000L, verbosep = FALSE,

                          include_soft_clips = 0L,

-                         exon_iit = NULL)

+                         exon_iit = NULL, xs = FALSE, read_pos = FALSE,

+                         min_base_quality = 0L, noncovered = FALSE)

 {

   if (!is(bamreader, "GmapBamReader"))

     stop("'bamreader' must be a GmapBamReader")

@@ -219,15 +240,6 @@ normArgSingleCharacterOrNULL <- function(x) {

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

   if (!is.null(desired_read_group) && !isSingleString(desired_read_group))

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

-  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")

-  }

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

         desired_read_group,

         normArgSingleInteger(alloclength),

@@ -245,7 +257,11 @@ normArgSingleCharacterOrNULL <- function(x) {

         normArgSingleInteger(blocksize),

         normArgTRUEorFALSE(verbosep),

         normArgSingleInteger(include_soft_clips),

-        normArgSingleCharacterOrNULL(exon_iit))

+        normArgSingleCharacterOrNULL(exon_iit),

+        normArgTRUEorFALSE(xs),

+        normArgTRUEorFALSE(read_pos),

+        normArgSingleInteger(min_base_quality),

+        normArgTRUEorFALSE(noncovered))

 }

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

@@ -271,7 +287,13 @@ variantSummaryColumnDescriptions <- function(read_pos_breaks) {

     read.pos.var.ref = "Variance in read position for the REF",

     mdfne = "Median distance from nearest end of read for the ALT",

     mdfne.ref = "Median distance from nearest end of read for the REF",

-    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)" )

+    codon.dir = "Direction of transcription for the codon. 0=positive, 1=negative, 2=NA (not a codon)",

+    del.count.plus = "plus strand deletion count",

+    del.count.minus = "minus strand deletion count",

+    count.xs.plus = "Plus strand XS counts",

+    count.xs.plus.ref = "Plus strand reference XS counts",

+    count.xs.minus = "Minus strand XS counts",

+    count.xs.minus.ref = "Minus strand reference XS counts")

   if (length(read_pos_breaks) > 0L) {

     break_desc <- paste0("Raw ALT count in read position range [",

                          head(read_pos_breaks, -1), ",",

@@ -186,9 +186,9 @@ normArgTRUEorFALSE <- function(x) {

   x

 }

-normArgSingleCharacter <- function(x) {

+normArgSingleCharacterOrNULL <- function(x) {

   name <- deparse(substitute(x))

-  if (!is(x, "character") || length(x) != 1)

+  if (!is.null(x) && (!is(x, "character") || length(x) != 1))

     stop("'", name, "' should be a single character value")

   x

 }

@@ -206,7 +206,7 @@ normArgSingleCharacter <- function(x) {

                          indels = FALSE,

                          blocksize = 1000L, verbosep = FALSE,

                          include_soft_clips = 0L,

-                         cds_iit)

+                         exon_iit = NULL)

 {

   if (!is(bamreader, "GmapBamReader"))

     stop("'bamreader' must be a GmapBamReader")

@@ -245,7 +245,7 @@ normArgSingleCharacter <- function(x) {

         normArgSingleInteger(blocksize),

         normArgTRUEorFALSE(verbosep),

         normArgSingleInteger(include_soft_clips),

-        normArgSingleCharacter(cds_iit))

+        normArgSingleCharacterOrNULL(exon_iit))

 }

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

@@ -9,11 +9,10 @@

 setClass("TallyIIT", representation(ptr = "externalptr",

                                     genome = "GmapGenome",

-                                    bam = "BamFile",

-                                    xs = "logical"))

+                                    bam = "BamFile"))

-TallyIIT <- function(ptr, genome, bam, xs) {

-  new("TallyIIT", ptr = ptr, genome = genome, bam = bam, xs = xs)

+TallyIIT <- function(ptr, genome, bam) {

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

 }

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

@@ -67,8 +66,9 @@ setMethod("bam_tally", "GmapBamReader",

             }

             param_list$genome <- NULL

+            

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

-                     as(x, "BamFile"), param_list$count_xs)

+                     as(x, "BamFile"))

           })

 variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

@@ -81,12 +81,12 @@ variantSummary <- 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, read_length, x@xs)

+                 NULL, read_length)

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

                    "n.read.pos", "n.read.pos.ref",

                    "raw.count", "raw.count.ref",

-                    "raw.count.total",

+                   "raw.count.total",

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

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

                    "mean.quality.ref",

@@ -94,9 +94,7 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

                    "count.minus", "count.minus.ref",

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

                    "read.pos.var", "read.pos.var.ref",

-                   "mdfne", "mdfne.ref",

-                   if (x@xs) c("count.xs.plus", "count.xs.plus.ref",

-                               "count.xs.minus", "count.xs.minus.ref"))

+                   "mdfne", "mdfne.ref", "codon.dir")

   break_names <- character()

   if (length(read_pos_breaks) > 0L) {

     read_pos_breaks <- as.integer(read_pos_breaks)

@@ -118,8 +116,9 @@ variantSummary <- function(x, read_pos_breaks = NULL, high_base_quality = 0L,

   genome <- genome(x)

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

   metacols <- DataFrame(tally[meta_names])

-  mcols(metacols) <-

-    variantSummaryColumnDescriptions(read_pos_breaks)[meta_names,,drop=FALSE]

+  mcols(metacols) <- variantSummaryColumnDescriptions(read_pos_breaks)

+

+  samecols = tally[["ref"]] != tally[["alt"]]

   gr <- with(tally,

              VRanges(seqnames,

                      IRanges(pos,

@@ -187,6 +186,13 @@ normArgTRUEorFALSE <- function(x) {

   x

 }

+normArgSingleCharacter <- function(x) {

+  name <- deparse(substitute(x))

+  if (!is(x, "character") || length(x) != 1)

+    stop("'", name, "' should be a single character value")

+  x

+}

+