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

         methods to work with GMAP and GSNAP from within R. In addition,

         it provides methods to tally alignment results on a

         per-nucleotide basis using the bam_tally tool.

-Version: 1.7.2

+Version: 1.7.3

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

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

         GenomicFeatures, Biostrings, VariantAnnotation (>= 1.11.1),

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

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

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

@@ -243,6 +244,8 @@ normArgTRUEorFALSE <- function(x) {

 variantSummaryColumnDescriptions <- function(read_pos_breaks) {

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

@@ -53,6 +53,8 @@ variantSummary(x, read_pos_breaks = NULL, high_base_quality = 0L,

   after quality filtering (except for indels, for which there is no

   quality filtering). The following \code{elementMetadata}

   columns are also present:

+  \item{n.read.pos}{The number of unique read positions for the alt allele.}

+  \item{n.read.pos.ref}{The number of unique read positions for the ref allele.}

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

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

   \item{raw.count.ref}{The number of reads with the reference allele.}

@@ -4,7 +4,7 @@

 #include "iit.h"

 #include "bytestream.h"

-enum { SEQNAMES, POS, REF, READ, COUNT, COUNT_REF,

+enum { SEQNAMES, POS, REF, READ, N_CYCLES, N_CYCLES_REF, COUNT, COUNT_REF,

        COUNT_TOTAL, HIGH_QUALITY, HIGH_QUALITY_REF, HIGH_QUALITY_TOTAL,

        MEAN_QUALITY, MEAN_QUALITY_REF, COUNT_PLUS, COUNT_PLUS_REF, COUNT_MINUS,

        COUNT_MINUS_REF, READ_POS_MEAN, READ_POS_MEAN_REF, READ_POS_VAR,

@@ -15,6 +15,8 @@ typedef struct TallyTable {

   int *pos;

   SEXP ref_R;

   SEXP read_R;

+  int *n_cycles;

+  int *n_cycles_ref;

   int *count;

   int *count_ref;

   int *count_total;

@@ -52,6 +54,8 @@ static SEXP R_TallyTable_new(int n_rows, int n_cycle_bins) {

   SET_VECTOR_ELT(tally_R, POS, allocVector(INTSXP, n_rows));

   SET_VECTOR_ELT(tally_R, REF, allocVector(STRSXP, n_rows));

   SET_VECTOR_ELT(tally_R, READ, allocVector(STRSXP, n_rows));

+  SET_VECTOR_ELT(tally_R, N_CYCLES, allocVector(INTSXP, n_rows));

+  SET_VECTOR_ELT(tally_R, N_CYCLES_REF, allocVector(INTSXP, n_rows));

   SET_VECTOR_ELT(tally_R, COUNT, allocVector(INTSXP, n_rows));

   SET_VECTOR_ELT(tally_R, COUNT_REF, allocVector(INTSXP, n_rows));

   SET_VECTOR_ELT(tally_R, COUNT_TOTAL, allocVector(INTSXP, n_rows));

@@ -88,6 +92,8 @@ static TallyTable *TallyTable_new(SEXP tally_R) {

   tally->pos = INTEGER(VECTOR_ELT(tally_R, POS));

   tally->ref_R = VECTOR_ELT(tally_R, REF);

   tally->read_R = VECTOR_ELT(tally_R, READ);

+  tally->n_cycles = INTEGER(VECTOR_ELT(tally_R, N_CYCLES));

+  tally->n_cycles_ref = INTEGER(VECTOR_ELT(tally_R, N_CYCLES_REF));

   tally->count = INTEGER(VECTOR_ELT(tally_R, COUNT));

   tally->count_ref = INTEGER(VECTOR_ELT(tally_R, COUNT_REF));

   tally->count_total = INTEGER(VECTOR_ELT(tally_R, COUNT_TOTAL));

@@ -123,7 +129,7 @@ read_total_counts(unsigned char **bytes, int row, int *count_total) {

 static void

 read_cycle_counts(unsigned char **bytes, int row, TallyParam param,

-                  double *read_pos_mean, double *read_pos_var,

+                  int *n_cycles, double *read_pos_mean, double *read_pos_var,

                   double *mdfne, int **cycle_bins)

 {

   int n_cycle_breaks = param.n_cycle_bins + 1;

@@ -133,9 +139,9 @@ read_cycle_counts(unsigned char **bytes, int row, TallyParam param,

   if (param.mdfne_buf != NULL) {

     memset(param.mdfne_buf, 0, sizeof(int) * midpoint);

   }

-  int n_cycles = read_int(bytes);

+  n_cycles[row] = read_int(bytes);

-  for (int index = 0; index < n_cycles; index++) {

+  for (int index = 0; index < n_cycles[row]; index++) {

     int cycle = abs(read_int(bytes));

     int count = read_int(bytes);

     int bin = 0;

@@ -145,7 +151,7 @@ read_cycle_counts(unsigned char **bytes, int row, TallyParam param,

     count_sum += count;

     weighted_sum += cycle * count;

     weighted_sum_sq += cycle * cycle * count;

-    if (param.cycle_breaks != NULL) {

+    if (param.n_cycle_bins > 0) {

       while(n_cycle_breaks > bin &&

             cycle > param.cycle_breaks[bin])

         bin++;

@@ -211,7 +217,7 @@ parse_indels(unsigned char *bytes, int row,

       SET_STRING_ELT(tally->read_R, row, R_BlankString);

       SET_STRING_ELT(tally->ref_R, row, seq_R);

     }

-    read_cycle_counts(&bytes, row, param,

+    read_cycle_counts(&bytes, row, param, tally->n_cycles,

                       tally->read_pos_mean, tally->read_pos_var,

                       tally->mdfne, tally->cycle_bins);

     /* quality is per-base and so not relevant for indels */

@@ -221,6 +227,7 @@ parse_indels(unsigned char *bytes, int row,

     tally->high_quality_ref[row] = NA_INTEGER;

     tally->high_quality_total[row] = NA_INTEGER;

     /* no position from which to tabulate the cycles */

+    tally->n_cycles_ref[row] = NA_INTEGER;

     tally->read_pos_mean_ref[row] = NA_REAL;

     tally->read_pos_var_ref[row] = NA_REAL;

     tally->mdfne_ref[row] = NA_REAL;

@@ -275,13 +282,14 @@ parse_alleles(unsigned char *bytes, int row, int ref_row,

                                      tally->count_plus, tally->count_minus,

                                      tally->count);

   for (int allele = 0; allele < n_alleles; allele++, row++) {

+    tally->n_cycles[row] = 0;

     for (int b = 0; b < param.n_cycle_bins; b++) {

       tally->cycle_bins[b][row] = 0;

     }

     tally->high_quality[row] = 0;

     tally->mean_quality[row] = R_NaN;

     if (tally->count[row] > 0) {

-      read_cycle_counts(&bytes, row, param,

+      read_cycle_counts(&bytes, row, param, tally->n_cycles,

                         tally->read_pos_mean, tally->read_pos_var,

                         tally->mdfne, tally->cycle_bins);

       read_quality_counts(&bytes, row, tally->high_quality, tally->mean_quality,

@@ -297,6 +305,7 @@ parse_alleles(unsigned char *bytes, int row, int ref_row,

     high_quality_total += tally->high_quality[r];

   }

   for (int r = ref_row; r < row; r++) {

+    tally->n_cycles_ref[r] = tally->n_cycles[ref_row];

     tally->count_total[r] = tally->count_total[ref_row];

     tally->mean_quality_ref[r] = tally->mean_quality[ref_row];

     tally->high_quality_ref[r] = tally->high_quality[ref_row];

@@ -313,6 +322,7 @@ parse_alleles(unsigned char *bytes, int row, int ref_row,

   if (have_ref_row) {

     /* clear the 'alt' columns for the 'ref' row with NAs */

     SET_STRING_ELT(tally->read_R, ref_row, NA_STRING);

+    tally->n_cycles[ref_row] = NA_INTEGER;

     tally->mean_quality[ref_row] = NA_REAL;

     tally->high_quality[ref_row] = NA_REAL;

     tally->count_plus[ref_row] = NA_INTEGER;