@@ -14,16 +14,14 @@

 }

 \usage{

-  BamTallyParam(genome, which = RangesList(),

+  BamTallyParam(genome, which = GRanges(),

                 desired_read_group = NULL,

-                read_pos_breaks = NULL,

-                high_base_quality = 0L,

                 minimum_mapq = 0L,

                 concordant_only = FALSE, unique_only = FALSE,

                 primary_only = FALSE, ignore_duplicates = FALSE,

                 min_depth = 0L, variant_strand = 0L,

                 ignore_query_Ns = FALSE,

-                indels = TRUE)

+                indels = FALSE)

 }

 \arguments{

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

@@ -33,11 +31,6 @@

   }

   \item{desired_read_group}{The name of the read group to which to limit

     the tallying; if not NULL, must be a single, non-NA string.}

-  \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.}

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

     at all.}

   \item{concordant_only}{Consider only what gnsap

@@ -12,54 +12,73 @@

 }

 \usage{

-GsnapParam(genome, unique_only = FALSE,

-   max_mismatches = NULL, suboptimal_levels = 0L,

-   mode = "standard", snps = NULL,

-   npaths = if (unique_only) 1L else 100L,

-   quiet_if_excessive = unique_only, nofails = unique_only,

-   split_output = !unique_only,

-   novelsplicing = FALSE, splicing = NULL, 

-   nthreads = 1L, part = NULL, batch = "2", ...)

+GsnapParam(genome, unique_only = FALSE, molecule = c("RNA", "DNA"),

+           max_mismatches = NULL,

+           suboptimal_levels = 0L, mode = "standard",

+           snps = NULL,

+           npaths = if (unique_only) 1L else 100L,

+           quiet_if_excessive = unique_only, nofails = unique_only,

+           split_output = !unique_only,

+           novelsplicing = FALSE, splicing = NULL, 

+           nthreads = 1L, part = NULL, batch = "2",

+           terminal_threshold = if (molecule == "DNA") 1000L else 2L,

+           gmap_mode = if (molecule == "DNA") "none" else

+                       "pairsearch,terminal,improve", ...)

 }

 \arguments{

-\item{genome}{A GmapGenome object to align against}

-\item{unique_only}{Whether only alignments with a unique match should be

-  output. The default is FALSE.}

-\item{max_mismatches}{The maximum number of mismatches to allow per

-  alignment. If NULL, then the value defaults to ((readlength + 2) / 12 - 2))}

-\item{suboptimal_levels}{Report suboptimal hits beyond best hit. The

-  default is 0L.}

-\item{mode}{The alignment mode. It can be "standard", "cmet-stranded",

-  "cmet-nonstranded", "atoi-stranded", or "atoi-nonstranded". The default 

-  is "standard".}

-\item{snps}{If not NULL, then a GmapSnps object. Provided SNPs will not

-  count as mismatches.}

-\item{npaths}{The maximum number of paths to print.}

-\item{quiet_if_excessive}{If more than maximum number of paths are

-  found, then no alignment from the read will be in the output.}

-\item{nofails}{Exclude failed alignments from output}

-\item{split_output}{Basename for multiple-file output, separately for nomapping,

-  halfmapping_uniq, halfmapping_mult, unpaired_uniq, unpaired_mult, paired_uniq,

-  paired_mult, concordant_uniq, and concordant_mult results (up to 9

-  files, or 10 if --fails-as-input is selected, or 3 for single-end reads)}

-\item{novelsplicing}{Logical indicating whether to look for novel

-  splicing events. FALSE is the default.}

-\item{splicing}{If not NULL, a GmapSplices object. NULL is the default.}

-\item{nthreads}{The number of worker threads gsnap should use to align.}

-\item{part}{If not NULL, then process only the i-th out of every n

-  sequences e.g., 0/100 or 99/100 (useful for distributing jobs to a

-  computer farm). If NULL, then all sequences are processed. NULL is the default.}

-\item{batch}{This argument allows control over gsnap's memory mapping

-  and allocation. The default is mode 2.

-  Mode 0: \{offsets=allocate, positions=mmap, genome=mmap\},

-  Mode 1: \{offsets=allocate, positions=mmap & preload,genome=mmap\},

-  Mode 2: \{offsets=allocate, positions=mmap & preload,genome=mmap &

-  preload\}, 

-  Mode 3: \{offsets=allocate, positions=allocate,genome=mmap & preload\},

-  Mode 4: \{offsets=allocate, positions=allocate,genome=allocate\}}

+  \item{genome}{A GmapGenome object to align against}

+  \item{unique_only}{Whether only alignments with a unique match should be

+    output. The default is FALSE.}

+  \item{molecule}{The type of molecule sequenced; used to determine

+    appropriate parameter defaults.}

+  \item{max_mismatches}{The maximum number of mismatches to allow per

+    alignment. If NULL, then the value defaults to ((readlength + 2) / 12 - 2))}

+  \item{suboptimal_levels}{Report suboptimal hits beyond best hit. The

+    default is 0L.}

+  \item{mode}{The alignment mode. It can be "standard", "cmet-stranded",

+    "cmet-nonstranded", "atoi-stranded", or "atoi-nonstranded". The default 

+    is "standard".}

+  \item{snps}{If not NULL, then a GmapSnps object. Provided SNPs will not

+    count as mismatches.}

+  \item{npaths}{The maximum number of paths to print.}

+  \item{quiet_if_excessive}{If more than maximum number of paths are

+    found, then no alignment from the read will be in the output.}

+  \item{nofails}{Exclude failed alignments from output}

+  \item{split_output}{Basename for multiple-file output, separately for

+    nomapping, halfmapping_uniq, halfmapping_mult, unpaired_uniq,

+    unpaired_mult, paired_uniq, paired_mult, concordant_uniq, and

+    concordant_mult results (up to 9 files, or 10 if --fails-as-input is

+    selected, or 3 for single-end reads)}

+  \item{novelsplicing}{Logical indicating whether to look for novel

+    splicing events. FALSE is the default.}

+  \item{splicing}{If not NULL, a GmapSplices object. NULL is the default.}

+  \item{nthreads}{The number of worker threads gsnap should use to align.}

+  \item{part}{If not NULL, then process only the i-th out of every n

+    sequences e.g., 0/100 or 99/100 (useful for distributing jobs to a

+    computer farm). If NULL, then all sequences are processed. NULL is

+    the default.}

+  \item{batch}{This argument allows control over gsnap's memory mapping

+    and allocation. The default is mode 2.

+    Mode 0: \{offsets=allocate, positions=mmap, genome=mmap\},

+    Mode 1: \{offsets=allocate, positions=mmap & preload,genome=mmap\},

+    Mode 2: \{offsets=allocate, positions=mmap & preload,genome=mmap &

+    preload\}, 

+    Mode 3: \{offsets=allocate, positions=allocate,genome=mmap & preload\},

+    Mode 4: \{offsets=allocate, positions=allocate,genome=allocate\}}

   \item{...}{Additional parameters for gsnap. See gsnap's full

-  documentation for those available.}

+    documentation for those available.}

+  \item{terminal_threshold}{If this number of mismatches is exceeded,

+    GSNAP will attempt to align from one of the sequence, eventually

+    giving up and discarding the rest of the sequence. This is called a

+    \dQuote{terminal alignment}. By setting this to a high value, we

+    have effectively disabled it for DNA, since terminal alignments were

+    motivated by splicing alignment problems and other special cases.

+  }

+  \item{gmap_mode}{Specifies the GMAP pipeline executed when GSNAP

+    delegates to GMAP (a Smith-Waterman aligner) in difficult cases. We

+    have disabled this for DNA, since such difficult cases are only

+    anticipated in the context of splicing or complex rearrangements.}

 }

 \seealso{

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

 \alias{bam_tally,BamFile-method}

 \alias{bam_tally,character-method}

 \alias{bam_tally,GmapBamReader-method}

+\alias{genome,TallyIIT-method}

 \alias{bam_tally}

+\alias{variantSummary}

 \title{Per-position Alignment Summaries}

@@ -17,41 +19,46 @@

 \usage{

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

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

+variantSummary(x, read_pos_breaks = NULL, high_base_quality = 0L,

+               keep_ref_rows = FALSE)

 }

 \arguments{

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

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

     parameters for the tally operation. }

+  \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.}

+  \item{keep_ref_rows}{Whether to keep the rows describing only the

+    reference calls, i.e., where ref and alt are the same. These are

+    useful when one needs the reference counts even when there are no

+    alts at that position.}

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

 }

 \value{

-  A \code{\link[GenomicRanges]{GRanges}}, with a range for each position

-  that passed the filters, and with the following \code{elementMetadata}

-  columns:

-  \item{location}{A string representation of the location, of the form

-    \dQuote{chr:pos}. This makes it easy, e.g., to check for the

-    presence of a variant in another result object.}

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

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

-    reference allele row.}

+  The \code{bam_tally} function returns an opaque pointer to a C-level

+  data structure with the class \dQuote{TallyIIT}. Currently, the only

+  operation applicable to this object is \code{variantSummary}.

+  

+  The \code{variantSummary} function returns

+  a \code{\link[VariantAnnotation]{VRanges}}, with a range for each position

+  that passed the filters. The depth columns correspond to the counts

+  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 cycles at which the alternate allele was

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

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

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

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

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

-  \item{count.total}{The total number of reads at that position,

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

+  \item{raw.count.total}{The total number of reads at that position,

     including reference and all alternates.}

-  \item{high.quality}{The number of reads for the alternate allele that were

-    above \code{high_quality_cutoff}, \code{NA} for the reference allele

-    row.}

-  \item{high.quality.ref}{The number of reads for the reference allele that were

-    above \code{high_quality_cutoff}.}

-  \item{mean.quality}{The mean mapping quality for the alternate allele,

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

   \item{mean.quality.ref}{The mean mapping quality for the reference

     allele.}

   \item{count.pos}{The number of positive strand reads for the alternate

@@ -62,9 +69,14 @@

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

   \item{count.neg.ref}{The number of negative strand reads for the reference

     allele.}

-

+  \item{read.pos.var}{Variance in the read positions for the alt allele.}

+  \item{read.pos.var.ref}{Variance in the read positions for the ref allele.}

+  

   An additional column is present for each bin formed by

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

 }

+\seealso{\code{tallyVariants} in the VariantTools package provides a

+  high-level wrapper for this functionality.}

+

 \author{Michael Lawrence}