diff --git a/MANUAL.markdown b/MANUAL.markdown index a0f454c..6ca9ad5 100644 --- a/MANUAL.markdown +++ b/MANUAL.markdown @@ -179,10 +179,10 @@ policy, which is similar to [Maq]'s default policy. number 5 or greater, set with [`-l`]) on the high-quality (left) end of the read. The first `L` bases are called the "seed". - 2. The sum of the [Phred quality] values at all mismatched positions - may not exceed `E` (set with [`-e`]). Where qualities are - unavailable (e.g. if the reads are from a FASTA file), the Phred - quality defaults to 40. + 2. The sum of the [Phred quality] values at *all* mismatched positions + (not just in the seed) may not exceed `E` (set with [`-e`]). Where + qualities are unavailable (e.g. if the reads are from a FASTA + file), the [Phred quality] defaults to 40. The [`-n`] option is mutually exclusive with the [`-v`] option. @@ -414,12 +414,12 @@ Colorspace Alignment [Colorspace alignment]: #colorspace-alignment -As of version 0.12.0, `bowtie` aligns colorspace reads against -colorspace references when [`-C`] is specified. Colorspace is the +As of version 0.12.0, `bowtie` can align colorspace reads against a +colorspace index when [`-C`] is specified. Colorspace is the characteristic output format of Applied Biosystems' SOLiD system. In a -colorspace read, each character is a color (rather than a nucleotide) +colorspace read, each character is a color rather than a nucleotide, where a color encodes a class of dinucleotides. E.g. the color blue -encodes "one of the dinucleotides: AA, CC, GG, TT." Colorspace has the +encodes any of the dinucleotides: AA, CC, GG, TT. Colorspace has the advantage of (often) being able to distinguish sequencing errors from SNPs once the read has been aligned. See ABI's [Principles of Di-Base Sequencing] application note for details. @@ -836,7 +836,21 @@ Default: off. #### Alignment - +
+ + +
+ +[`-v`]: #bowtie-options-v + + -v + + + +Report alignments with at most `` mismatches. [`-e`] and [`-l`] +options are ignored and quality values have no effect on what +alignments are valid. [`-v`] is mutually exclusive with [`-n`]. + +
[`-n`/`--seedmms`]: #bowtie-options-n [`-n`]: #bowtie-options-n @@ -859,8 +873,9 @@ exclusive with the [`-v`] option. -Maximum permitted total of quality values at mismatched read positions. -The default is 70. Like [Maq], `bowtie` rounds quality values to the +Maximum permitted total of quality values at *all* mismatched read +positions throughout the entire alignment, not just in the "seed". The +default is 70. Like [Maq], `bowtie` rounds quality values to the nearest 10 and saturates at 30; rounding can be disabled with [`--nomaqround`]. @@ -891,18 +906,6 @@ internally rounds values to the nearest 10, with a maximum of 30. By default, `bowtie` also rounds this way. [`--nomaqround`] prevents this rounding in `bowtie`. -
- -[`-v`]: #bowtie-options-v - - -v - - - -Report alignments with at most `` mismatches. [`-e`] and [`-l`] -options are ignored and quality values have no effect on what -alignments are valid. [`-v`] is mutually exclusive with [`-n`]. -
[`-I`/`--minins`]: #bowtie-options-I @@ -1123,13 +1126,12 @@ tuning] section for details). Behaves like [`-m`] except that if a read has more than `` reportable alignments, one is reported at random. In [default -output mode], the selected alignment's 7th column is set to `` to -indicate the read has at least `` valid alignments. In +output mode], the selected alignment's 7th column is set to ``+1 to +indicate the read has at least ``+1 valid alignments. In [`-S`/`--sam`] mode, the selected alignment is given a `MAPQ` (mapping -quality) of 0. Randomly-selected alignments do not count toward the -"reads with at least one reported alignment" total reported by -`bowtie`. This option requires [`--best`] mode; if `-M` is specified -without [`--best`], [`--best`] is enabled automatically. +quality) of 0 and the `XM:I` field is set to ``+1. This option +requires [`--best`]; if specified without [`--best`], [`--best`] is enabled +automatically. [default output mode]: #default-bowtie-output @@ -1173,36 +1175,8 @@ best stratum. By default, Bowtie reports all reportable alignments regardless of whether they fall into multiple strata. When [`--strata`] is specified, [`--best`] must also be specified. -
- -[`--snpphred`]: #bowtie-options-snpphred - - --snpphred - - - -When decoding colorspace alignments, use `` as the SNP penalty. -This should be set to the user's best guess of the true ratio of SNPs -per base in the subject genome, converted to the [Phred quality] scale. -E.g., if the user expects about 1 SNP every 1,000 positions, -`--snpphred` should be set to 30 (which is also the default). To -specify the fraction directly, use [`--snpfrac`]. - -
- -[`--snpfrac`]: #bowtie-options-snpfrac - - --snpfrac - - - -When decoding colorspace alignments, use `` as the estimated ratio -of SNPs per base. For best decoding results, this should be set to the -user's best guess of the true ratio. `bowtie` internally converts the -ratio to a [Phred quality], and behaves as if that quality had been set -via the [`--snpphred`] option. Default: 0.001. - -
+
#### Output @@ -1335,47 +1309,101 @@ quality fields will be omitted. See [Default Bowtie output] for field descriptions. This option is ignored if the output mode is [`-S`/`--sam`]. - + + -[`--colseq`]: #bowtie-options-colseq +[`--fullref`]: #bowtie-options-fullref + + --fullref + + + +Print the full refernce sequence name, including whitespace, in +alignment output. By default `bowtie` prints everything up to but not +including the first whitespace. + + + +#### Colorspace + + + + + + +
+ +[`--snpphred`]: #bowtie-options-snpphred - --colseq + --snpphred + + + +When decoding colorspace alignments, use `` as the SNP penalty. +This should be set to the user's best guess of the true ratio of SNPs +per base in the subject genome, converted to the [Phred quality] scale. +E.g., if the user expects about 1 SNP every 1,000 positions, +`--snpphred` should be set to 30 (which is also the default). To +specify the fraction directly, use [`--snpfrac`]. + +
+ +[`--snpfrac`]: #bowtie-options-snpfrac + + --snpfrac + + + +When decoding colorspace alignments, use `` as the estimated ratio +of SNPs per base. For best decoding results, this should be set to the +user's best guess of the true ratio. `bowtie` internally converts the +ratio to a [Phred quality], and behaves as if that quality had been set +via the [`--snpphred`] option. Default: 0.001. + +
+ +[`--col-cseq`]: #bowtie-options-col-cseq + + --col-cseq If reads are in colorspace and the [default output mode] is active, -`--colseq` causes the reads' color sequence to appear in the +`--col-cseq` causes the reads' color sequence to appear in the read-sequence column (column 5) instead of the decoded nucleotide sequence. See the [Decoding colorspace alignments] section for details about decoding. This option is ignored in [`-S`/`--sam`] mode. -
+
-[`--colqual`]: #bowtie-options-colqual +[`--col-cqual`]: #bowtie-options-col-cqual - --colqual + --col-cqual If reads are in colorspace and the [default output mode] is active, -`--colqual` causes the reads' original (color) quality sequence to +`--col-cqual` causes the reads' original (color) quality sequence to appear in the quality column (column 6) instead of the decoded qualities. See the [Colorspace alignment] section for details about decoding. This option is ignored in [`-S`/`--sam`] mode. -
+
-[`--fullref`]: #bowtie-options-fullref +[`--col-keepends`]: #bowtie-options-col-keepends - --fullref + --col-keepends -Print the full refernce sequence name, including whitespace, in -alignment output. By default `bowtie` prints everything up to but not -including the first whitespace. +When decoding colorpsace alignments, `bowtie` trims off a nucleotide +and quality from the left and right edges of the alignment. This is +because those nucleotides are supported by only one color, in contrast +to the middle nucleotides which are supported by two. Specify +`--col-keepends` to keep the extreme-end nucleotides and qualities. -
+ + #### SAM @@ -1391,12 +1419,13 @@ including the first whitespace. Print alignments in [SAM] format. See the [SAM output] section of the -manual for details. To suppress all SAM headers, use [`--sam-nohead`]. -To suppress just the `@SQ` headers (e.g. if the alignment is against a -very large number of reference sequences), use [`--sam-nosq`]. -`bowtie` does not write BAM files directly, but SAM output can be -converted to BAM on the fly by piping `bowtie`'s output to -`samtools view`. [`-S`/`--sam`] is not compatible with [`--refout`]. +manual for details. To suppress all SAM headers, use [`--sam-nohead`] +in addition to `-S/--sam`. To suppress just the `@SQ` headers (e.g. if +the alignment is against a very large number of reference sequences), +use [`--sam-nosq`] in addition to `-S/--sam`. `bowtie` does not write +BAM files directly, but SAM output can be converted to BAM on the fly +by piping `bowtie`'s output to `samtools view`. [`-S`/`--sam`] is not +compatible with [`--refout`]. [SAM output]: #sam-bowtie-output @@ -1421,6 +1450,8 @@ See the [SAM Spec][SAM] for details about the `MAPQ` field Default: 255. Suppress header lines (starting with `@`) when output is [`-S`/`--sam`]. +This must be specified *in addition to* [`-S`/`--sam`]. `--sam-nohead` +is ignored unless [`-S`/`--sam`] is also specified. @@ -1430,7 +1461,9 @@ Suppress header lines (starting with `@`) when output is [`-S`/`--sam`]. -Suppress `@SQ` header lines when output is [`-S`/`--sam`]. +Suppress `@SQ` header lines when output is [`-S`/`--sam`]. This must be +specified *in addition to* [`-S`/`--sam`]. `--sam-nosq` is ignored +unless [`-S`/`--sam`] is also specified. @@ -1445,7 +1478,8 @@ field on the `@RG` header line. Specify `--sam-RG` multiple times to set multiple fields. See the [SAM Spec][SAM] for details about what fields are legal. Note that, if any `@RG` fields are set using this option, the `ID` and `SM` fields must both be among them to make the `@RG` line -legal according to the [SAM Spec][SAM]. +legal according to the [SAM Spec][SAM]. `--sam-RG` is ignored unless +[`-S`/`--sam`] is also specified. @@ -1590,7 +1624,7 @@ Default `bowtie` output If the read was in colorspace, then the sequence shown in this column is the sequence of *decoded nucleotides*, not the original colors. See the [Colorspace alignment] section for details about - decoding. To display colors instead, use the [`--colseq`] option. + decoding. To display colors instead, use the [`--col-cseq`] option. 6. ASCII-encoded read qualities (reversed if orientation is `-`). The encoded quality values are on the Phred scale and the encoding is @@ -1599,9 +1633,12 @@ Default `bowtie` output If the read was in colorspace, then the qualities shown in this column are the *decoded qualities*, not the original qualities. See the [Colorspace alignment] section for details about decoding. - To display colors instead, use the [`--colqual`] option. + To display colors instead, use the [`--col-cqual`] option. -7. Number of other instances where the same read aligns against the +7. If [`-M`] was specified and the [`-M`] ceiling was exceeded for this + read, this column contains the + +Number of other instances where the same read aligns against the same reference characters as were aligned against in this alignment. This is *not* the number of other places the read aligns with the same number of mismatches. The number in this column is generally @@ -1785,9 +1822,15 @@ right, the fields are: For a read with no reported alignments, `` is 0 if the read had - no alignments, or 1 if the read had alignments that were suppressed - by the [`-m`] option. - + no alignments. If [`-m`] was specified and the read's alignments + were supressed because the [`-m`] ceiling was exceeded, `` equals + the [`-m`] ceiling + 1, to indicate that there were at least that + many valid alignments (but they were suppressed). In [`-M`] mode, if + the alignment was randomly selected because the [`-M`] ceiling was + exceeded, `` equals the [`-M`] ceiling + 1, to indicate that there + were at least that many valid alignments (but only one was + reported). + [SAM format specification]: http://samtools.sf.net/SAM1.pdf diff --git a/aligner.h b/aligner.h index fabc936..ba41827 100644 --- a/aligner.h +++ b/aligner.h @@ -20,6 +20,7 @@ #include "ref_aligner.h" #include "reference.h" #include "aligner_metrics.h" +#include "search_globals.h" /** * State machine for carrying out an alignment, which usually consists @@ -332,7 +333,6 @@ class UnpairedAlignerV2 : public Aligner { HitSinkPerThread* sinkPt, vector >& os, const BitPairReference* refs, - int snpPhred, bool rangeMode, bool verbose, bool quiet, @@ -342,7 +342,6 @@ class UnpairedAlignerV2 : public Aligner { AlignerMetrics *metrics = NULL) : Aligner(true, rangeMode), refs_(refs), - snpPhred_(snpPhred), doneFirst_(true), firstIsFw_(true), chase_(false), @@ -422,7 +421,8 @@ class UnpairedAlignerV2 : public Aligner { (ebwtFw? &bufa_->qualRev : &bufa_->qual), &bufa_->name, bufa_->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, ra.ebwt->rmap(), ebwtFw, @@ -518,8 +518,6 @@ class UnpairedAlignerV2 : public Aligner { // Reference sequences (needed for colorspace decoding) const BitPairReference* refs_; - int snpPhred_; - // Progress state bool doneFirst_; bool firstIsFw_; @@ -579,7 +577,6 @@ class PairedBWAlignerV1 : public Aligner { uint32_t mixedThresh, uint32_t mixedAttemptLim, const BitPairReference* refs, - int snpPhred, bool rangeMode, bool verbose, bool quiet, @@ -587,7 +584,7 @@ class PairedBWAlignerV1 : public Aligner { ChunkPool *pool, int *btCnt) : Aligner(true, rangeMode), - refs_(refs), snpPhred_(snpPhred), + refs_(refs), patsrc_(NULL), qlen1_(0), qlen2_(0), doneFw_(true), doneFwFirst_(true), chase1Fw_(false), chase1Rc_(false), @@ -835,7 +832,8 @@ class PairedBWAlignerV1 : public Aligner { (ebwtFwL? &bufL->qualRev : &bufL->qual), &bufL->name, bufL->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, rmap, ebwtFwL, @@ -865,7 +863,8 @@ class PairedBWAlignerV1 : public Aligner { (ebwtFwR? &bufR->qualRev : &bufR->qual), &bufR->name, bufR->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, rmap, ebwtFwR, @@ -1249,8 +1248,6 @@ class PairedBWAlignerV1 : public Aligner { const BitPairReference* refs_; - int snpPhred_; - PatternSourcePerThread *patsrc_; uint32_t qlen1_; uint32_t qlen2_; @@ -1431,7 +1428,6 @@ class PairedBWAlignerV2 : public Aligner { uint32_t maxInsert, uint32_t mixedAttemptLim, const BitPairReference* refs, - int snpPhred, bool rangeMode, bool verbose, bool quiet, @@ -1440,7 +1436,6 @@ class PairedBWAlignerV2 : public Aligner { int *btCnt) : Aligner(true, rangeMode), refs_(refs), - snpPhred_(snpPhred), patsrc_(NULL), qlen1_(0), qlen2_(0), chase_(false), @@ -1670,7 +1665,8 @@ class PairedBWAlignerV2 : public Aligner { (ebwtFwL? &bufL->qualRev : &bufL->qual), &bufL->name, bufL->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, rmap, ebwtFwL, @@ -1700,7 +1696,8 @@ class PairedBWAlignerV2 : public Aligner { (ebwtFwR? &bufR->qualRev : &bufR->qual), &bufR->name, bufR->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, rmap, ebwtFwR, @@ -1742,7 +1739,8 @@ class PairedBWAlignerV2 : public Aligner { (ebwtFw? &buf->qualRev : &buf->qual), &buf->name, buf->color, - snpPhred_, + colorExEnds, + snpPhred, refs_, r.ebwt->rmap(), ebwtFw, @@ -1907,8 +1905,6 @@ class PairedBWAlignerV2 : public Aligner { const BitPairReference* refs_; - int snpPhred_; - PatternSourcePerThread *patsrc_; uint32_t qlen1_, qlen2_; bool chase_; diff --git a/aligner_0mm.h b/aligner_0mm.h index 312c2d3..603f62e 100644 --- a/aligner_0mm.h +++ b/aligner_0mm.h @@ -25,7 +25,6 @@ class UnpairedExactAlignerV1Factory : public AlignerFactory { UnpairedExactAlignerV1Factory( Ebwt >& ebwtFw, Ebwt >* ebwtBw, - int snpPhred, bool doFw, bool doRc, HitSink& sink, @@ -45,7 +44,6 @@ class UnpairedExactAlignerV1Factory : public AlignerFactory { uint32_t seed) : ebwtFw_(ebwtFw), ebwtBw_(ebwtBw), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), sink_(sink), sinkPtFactory_(sinkPtFactory), @@ -112,14 +110,13 @@ class UnpairedExactAlignerV1Factory : public AlignerFactory { return new UnpairedAlignerV2( params, dr, rchase, - sink_, sinkPtFactory_, sinkPt, os_, refs_, snpPhred_, + sink_, sinkPtFactory_, sinkPt, os_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL, NULL); } private: Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; - int snpPhred_; bool doFw_; bool doRc_; HitSink& sink_; @@ -151,7 +148,6 @@ class PairedExactAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw, Ebwt >* ebwtBw, bool color, - int snpPhred, bool doFw, bool doRc, bool v1, @@ -181,7 +177,6 @@ class PairedExactAlignerV1Factory : public AlignerFactory { uint32_t seed) : ebwtFw_(ebwtFw), color_(color), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), v1_(v1), @@ -328,7 +323,7 @@ class PairedExactAlignerV1Factory : public AlignerFactory { refAligner, rchase, sink_, sinkPtFactory_, sinkPt, mate1fw_, mate2fw_, peInner_, peOuter_, dontReconcile_, symCeil_, mixedThresh_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, verbose_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); return al; } else { @@ -344,7 +339,7 @@ class PairedExactAlignerV1Factory : public AlignerFactory { rchase, sink_, sinkPtFactory_, sinkPt, sinkPtSe1, sinkPtSe2, mate1fw_, mate2fw_, peInner_, peOuter_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); delete drVec; return al; @@ -355,7 +350,6 @@ class PairedExactAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; bool color_; - int snpPhred_; bool doFw_; bool doRc_; bool v1_; diff --git a/aligner_1mm.h b/aligner_1mm.h index 236806f..530464f 100644 --- a/aligner_1mm.h +++ b/aligner_1mm.h @@ -25,7 +25,6 @@ class Unpaired1mmAlignerV1Factory : public AlignerFactory { Unpaired1mmAlignerV1Factory( Ebwt >& ebwtFw, Ebwt >* ebwtBw, - int snpPhred, bool doFw, bool doRc, HitSink& sink, @@ -45,7 +44,6 @@ class Unpaired1mmAlignerV1Factory : public AlignerFactory { uint32_t seed) : ebwtFw_(ebwtFw), ebwtBw_(ebwtBw), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), sink_(sink), @@ -149,14 +147,13 @@ class Unpaired1mmAlignerV1Factory : public AlignerFactory { // Set up the aligner return new UnpairedAlignerV2( params, dr, rchase, - sink_, sinkPtFactory_, sinkPt, os_, refs_, snpPhred_, + sink_, sinkPtFactory_, sinkPt, os_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL, NULL); } private: Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; - const int snpPhred_; bool doFw_; bool doRc_; HitSink& sink_; @@ -188,7 +185,6 @@ class Paired1mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw, Ebwt >* ebwtBw, bool color, - int snpPhred, bool doFw, bool doRc, bool v1, @@ -219,7 +215,6 @@ class Paired1mmAlignerV1Factory : public AlignerFactory { ebwtFw_(ebwtFw), ebwtBw_(ebwtBw), color_(color), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), v1_(v1), @@ -440,7 +435,7 @@ class Paired1mmAlignerV1Factory : public AlignerFactory { refAligner, rchase, sink_, sinkPtFactory_, sinkPt, mate1fw_, mate2fw_, peInner_, peOuter_, dontReconcile_, symCeil_, mixedThresh_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, verbose_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); delete dr1FwVec; delete dr1RcVec; @@ -456,7 +451,7 @@ class Paired1mmAlignerV1Factory : public AlignerFactory { sinkPt, sinkPtSe1, sinkPtSe2, mate1fw_, mate2fw_, peInner_, peOuter_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); delete dr1FwVec; return al; @@ -467,7 +462,6 @@ class Paired1mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; bool color_; - const int snpPhred_; bool doFw_; bool doRc_; bool v1_; diff --git a/aligner_23mm.h b/aligner_23mm.h index 15fe14c..e79b396 100644 --- a/aligner_23mm.h +++ b/aligner_23mm.h @@ -26,7 +26,6 @@ class Unpaired23mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw, Ebwt >* ebwtBw, bool two, - int snpPhred, bool doFw, bool doRc, HitSink& sink, @@ -47,7 +46,6 @@ class Unpaired23mmAlignerV1Factory : public AlignerFactory { ebwtFw_(ebwtFw), ebwtBw_(ebwtBw), two_(two), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), sink_(sink), sinkPtFactory_(sinkPtFactory), @@ -218,7 +216,7 @@ class Unpaired23mmAlignerV1Factory : public AlignerFactory { return new UnpairedAlignerV2( params, dr, rchase, - sink_, sinkPtFactory_, sinkPt, os_, refs_, snpPhred_, + sink_, sinkPtFactory_, sinkPt, os_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL, NULL); } @@ -226,7 +224,6 @@ class Unpaired23mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; bool two_; - const int snpPhred_; bool doFw_; bool doRc_; HitSink& sink_; @@ -259,7 +256,6 @@ class Paired23mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw, Ebwt >* ebwtBw, bool color, - int snpPhred, bool doFw, bool doRc, bool v1, @@ -291,7 +287,6 @@ class Paired23mmAlignerV1Factory : public AlignerFactory { ebwtFw_(ebwtFw), ebwtBw_(ebwtBw), color_(color), - snpPhred_(snpPhred), doFw_(doFw), doRc_(doRc), v1_(v1), @@ -636,7 +631,7 @@ class Paired23mmAlignerV1Factory : public AlignerFactory { refAligner, rchase, sink_, sinkPtFactory_, sinkPt, mate1fw_, mate2fw_, peInner_, peOuter_, dontReconcile_, symCeil_, mixedThresh_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, verbose_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); delete dr1FwVec; delete dr1RcVec; @@ -652,7 +647,7 @@ class Paired23mmAlignerV1Factory : public AlignerFactory { sinkPt, sinkPtSe1, sinkPtSe2, mate1fw_, mate2fw_, peInner_, peOuter_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, INT_MAX, pool_, NULL); delete dr1FwVec; return al; @@ -663,7 +658,6 @@ class Paired23mmAlignerV1Factory : public AlignerFactory { Ebwt >& ebwtFw_; Ebwt >* ebwtBw_; bool color_; - int snpPhred_; bool doFw_; bool doRc_; bool v1_; diff --git a/aligner_seed_mm.h b/aligner_seed_mm.h index 737a2a3..af11d06 100644 --- a/aligner_seed_mm.h +++ b/aligner_seed_mm.h @@ -31,7 +31,6 @@ class UnpairedSeedAlignerFactory : public AlignerFactory { uint32_t seedMms, uint32_t seedLen, int qualCutoff, - int snpPhred, int maxBts, HitSink& sink, const HitSinkPerThreadFactory& sinkPtFactory, @@ -55,7 +54,6 @@ class UnpairedSeedAlignerFactory : public AlignerFactory { seedMms_(seedMms), seedLen_(seedLen), qualCutoff_(qualCutoff), - snpPhred_(snpPhred), maxBts_(maxBts), sink_(sink), sinkPtFactory_(sinkPtFactory), @@ -537,7 +535,7 @@ class UnpairedSeedAlignerFactory : public AlignerFactory { return new UnpairedAlignerV2( params, dr, rchase, - sink_, sinkPtFactory_, sinkPt, os_, refs_, snpPhred_, + sink_, sinkPtFactory_, sinkPt, os_, refs_, rangeMode_, verbose_, quiet_, maxBts_, pool_, btCnt, metrics_); } @@ -550,7 +548,6 @@ class UnpairedSeedAlignerFactory : public AlignerFactory { const uint32_t seedMms_; const uint32_t seedLen_; const int qualCutoff_; - const int snpPhred_; const int maxBts_; HitSink& sink_; const HitSinkPerThreadFactory& sinkPtFactory_; @@ -587,7 +584,6 @@ class PairedSeedAlignerFactory : public AlignerFactory { uint32_t seedMms, uint32_t seedLen, int qualCutoff, - int snpPhred, int maxBts, HitSink& sink, const HitSinkPerThreadFactory& sinkPtFactory, @@ -622,7 +618,6 @@ class PairedSeedAlignerFactory : public AlignerFactory { seedMms_(seedMms), seedLen_(seedLen), qualCutoff_(qualCutoff), - snpPhred_(snpPhred), maxBts_(maxBts), sink_(sink), sinkPtFactory_(sinkPtFactory), @@ -1335,7 +1330,7 @@ class PairedSeedAlignerFactory : public AlignerFactory { refAligner, rchase, sink_, sinkPtFactory_, sinkPt, mate1fw_, mate2fw_, peInner_, peOuter_, dontReconcile_, symCeil_, mixedThresh_, mixedAttemptLim_, refs_, - snpPhred_, rangeMode_, verbose_, quiet_, maxBts_, pool_, + rangeMode_, verbose_, quiet_, maxBts_, pool_, btCnt); delete dr1FwVec; delete dr1RcVec; @@ -1349,7 +1344,7 @@ class PairedSeedAlignerFactory : public AlignerFactory { new TCostAwareRangeSrcDr(strandFix_, dr1FwVec, verbose_, quiet_, true), refAligner, rchase, sink_, sinkPtFactory_, sinkPt, sinkPtSe1, sinkPtSe2, mate1fw_, mate2fw_, peInner_, peOuter_, - mixedAttemptLim_, refs_, snpPhred_, rangeMode_, verbose_, + mixedAttemptLim_, refs_, rangeMode_, verbose_, quiet_, maxBts_, pool_, btCnt); delete dr1FwVec; return al; @@ -1366,7 +1361,6 @@ class PairedSeedAlignerFactory : public AlignerFactory { const uint32_t seedMms_; const uint32_t seedLen_; const int qualCutoff_; - const int snpPhred_; const int maxBts_; HitSink& sink_; const HitSinkPerThreadFactory& sinkPtFactory_; diff --git a/ebwt.h b/ebwt.h index 7a4500b..9545fee 100644 --- a/ebwt.h +++ b/ebwt.h @@ -990,8 +990,8 @@ class Ebwt { // Searching and reporting void joinedToTextOff(uint32_t qlen, uint32_t off, uint32_t& tidx, uint32_t& textoff, uint32_t& tlen) const; - inline bool report(const String& query, String* quals, String* name, bool color, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, const std::vector& refcs, size_t numMms, uint32_t off, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, const EbwtSearchParams& params) const; - inline bool reportChaseOne(const String& query, String* quals, String* name, bool color, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, const std::vector& refcs, size_t numMms, uint32_t i, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, const EbwtSearchParams& params, SideLocus *l = NULL) const; + inline bool report(const String& query, String* quals, String* name, bool color, bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, const std::vector& refcs, size_t numMms, uint32_t off, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, const EbwtSearchParams& params) const; + inline bool reportChaseOne(const String& query, String* quals, String* name, bool color, bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, const std::vector& refcs, size_t numMms, uint32_t i, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, uint16_t cost, const EbwtSearchParams& params, SideLocus *l = NULL) const; inline bool reportReconstruct(const String& query, String* quals, String* name, String& lbuf, String& rbuf, const uint32_t *mmui32, const char* refcs, size_t numMms, uint32_t i, uint32_t top, uint32_t bot, uint32_t qlen, int stratum, const EbwtSearchParams& params, SideLocus *l = NULL) const; inline int rowL(const SideLocus& l) const; inline uint32_t countUpTo(const SideLocus& l, int c) const; @@ -1143,6 +1143,7 @@ class EbwtSearchParams { String* quals, // read quality values String* name, // read name bool color, // true -> read is colorspace + bool colExEnds, // true -> exclude nucleotides at extreme ends after decoding int snpPhred, // penalty for a SNP const BitPairReference* ref, // reference (= NULL if not necessary) const ReferenceMap* rmap, // map to another reference coordinate system @@ -1210,7 +1211,8 @@ class EbwtSearchParams { } assert(ref != NULL); char read[1024]; - char rf[1024]; + uint32_t rfbuf[(1024+16)/4]; + ASSERT_ONLY(char rfbuf2[1024]); char qual[1024]; char ns[1024]; char cmm[1024]; @@ -1222,17 +1224,19 @@ class EbwtSearchParams { size_t readf = seqan::length(hit.patSeq); size_t refi = 0; size_t reff = readf + 1; - // The "Phred-scaled probability of a mutation" a la Maq/BWA bool maqRound = false; for(size_t i = 0; i < qlen + 1; i++) { if(i < qlen) { read[i] = (int)hit.patSeq[i]; qual[i] = mmPenalty(maqRound, phredCharToPhredQual(hit.quals[i])); } - int rc = ref->getBase(h.first, h.second + i); - assert_geq(rc, 0); - assert_lt(rc, 4); - rf[i] = (1 << rc); + ASSERT_ONLY(rfbuf2[i] = ref->getBase(h.first, h.second + i)); + } + int offset = ref->getStretch(rfbuf, h.first, h.second, qlen + 1); + char *rf = (char*)rfbuf + offset; + for(size_t i = 0; i < qlen + 1; i++) { + assert_eq(rf[i], rfbuf2[i]); + rf[i] = (1 << rf[i]); } decodeHit( read, // ASCII colors, '0', '1', '2', '3', '.' @@ -1248,36 +1252,45 @@ class EbwtSearchParams { nmm, // where nucleotide mismatches are in the string cmms, // number of color mismatches nmms);// number of nucleotide mismatches - seqan::resize(hit.patSeq, qlen+1); - seqan::resize(hit.quals, qlen+1); - hit.refcs.resize(qlen+1, 0); - for(size_t i = 0; i < qlen + 1; i++) { + size_t nqlen = qlen + (colExEnds ? -1 : 1); + seqan::resize(hit.patSeq, nqlen); + seqan::resize(hit.quals, nqlen); + hit.refcs.resize(nqlen); + size_t lo = colExEnds ? 1 : 0; + size_t hi = colExEnds ? qlen : qlen+1; + size_t destpos = 0; + for(size_t i = lo; i < hi; i++, destpos++) { // Set sequence character assert_leq(ns[i], 4); assert_geq(ns[i], 0); - hit.patSeq[i] = (Dna5)(int)ns[i]; + hit.patSeq[destpos] = (Dna5)(int)ns[i]; // Set initial quality - hit.quals[i] = '!'; + hit.quals[destpos] = '!'; // Color mismatches penalize quality if(i > 0) { - if(cmm[i-1] == 'M') hit.quals[i] += qual[i-1]; - else hit.quals[i] -= qual[i-1]; + if(cmm[i-1] == 'M') hit.quals[destpos] += qual[i-1]; + else hit.quals[destpos] -= qual[i-1]; } if(i < qlen) { - if(cmm[i] == 'M') hit.quals[i] += qual[i]; - else hit.quals[i] -= qual[i]; + if(cmm[i] == 'M') hit.quals[destpos] += qual[i]; + else hit.quals[destpos] -= qual[i]; + } + if(hit.quals[destpos] < '!') { + hit.quals[destpos] = '!'; } - if(hit.quals[i] < '!') hit.quals[i] = '!'; if(nmm[i] != 'M') { - uint32_t off = i; - if(!_fw) off = (qlen+1) - i - 1; - assert_lt(off, qlen+1); + uint32_t off = i - (colExEnds? 1:0); + if(!_fw) off = nqlen - off - 1; + assert_lt(off, nqlen); hit.mms.set(off); hit.refcs[off] = "ACGT"[ref->getBase(h.first, h.second+i)]; } } - qlen++; - mlen++; + if(colExEnds) { + qlen--; mlen--; + } else { + qlen++; mlen++; + } } else { // Turn the mmui32 and refcs arrays into the mm FixedBitset and // the refc vector @@ -2354,6 +2367,7 @@ inline bool Ebwt::report(const String& query, String* quals, String* name, bool color, + bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, @@ -2381,6 +2395,7 @@ inline bool Ebwt::report(const String& query, quals, // read quality values name, // read name color, // true -> read is colorspace + colExEnds, // true -> exclude nucleotides on ends snpPhred, // phred probability of SNP ref, // reference sequence rmap_, // map to another reference coordinate system @@ -2416,6 +2431,7 @@ inline bool Ebwt::reportChaseOne(const String& query, String* quals, String* name, bool color, + bool colExEnds, int snpPhred, const BitPairReference* ref, const std::vector& mmui32, @@ -2473,9 +2489,9 @@ inline bool Ebwt::reportChaseOne(const String& query, assert_eq(rcoff, off); } #endif - return report(query, quals, name, color, snpPhred, ref, mmui32, - refcs, numMms, off, top, bot, qlen, stratum, cost, - params); + return report(query, quals, name, color, colExEnds, snpPhred, ref, + mmui32, refcs, numMms, off, top, bot, qlen, stratum, + cost, params); } /** diff --git a/ebwt_search.cpp b/ebwt_search.cpp index a81c6e7..5e7384d 100644 --- a/ebwt_search.cpp +++ b/ebwt_search.cpp @@ -128,9 +128,10 @@ static bool fuzzy; static bool fullRef; static bool samNoHead; // don't print any header lines in SAM output static bool samNoSQ; // don't print @SQ header lines -static bool color; +static bool color; // true -> inputs are colorspace +bool colorExEnds; // true -> nucleotides on either end of decoded cspace alignment should be excluded static string rgs; // SAM outputs for @RG header line -static int snpPhred; // probability of SNP, for scoring colorspace alignments +int snpPhred; // probability of SNP, for scoring colorspace alignments static Bitset suppressOuts(64); // output fields to suppress static bool sampleMax; // whether to report a random alignment when maxed-out via -m/-M static int defaultMapq; // default mapping quality to print in SAM mode @@ -231,6 +232,7 @@ static void resetOptions() { samNoHead = false; // don't print any header lines in SAM output samNoSQ = false; // don't print @SQ header lines color = false; // don't align in colorspace by default + colorExEnds = true; // true -> nucleotides on either end of decoded cspace alignment should be excluded rgs = ""; // SAM outputs for @RG header line snpPhred = 30; // probability of SNP, for scoring colorspace alignments suppressOuts.clear(); // output fields to suppress @@ -320,7 +322,8 @@ enum { ARG_DEFAULT_MAPQ, ARG_COLOR_SEQ, ARG_COLOR_QUAL, - ARG_COST + ARG_COST, + ARG_COLOR_KEEP_ENDS }; static struct option long_options[] = { @@ -421,8 +424,9 @@ static struct option long_options[] = { {(char*)"snpfrac", required_argument, 0, ARG_SNPFRAC}, {(char*)"suppress", required_argument, 0, ARG_SUPPRESS_FIELDS}, {(char*)"mapq", required_argument, 0, ARG_DEFAULT_MAPQ}, - {(char*)"colseq", no_argument, 0, ARG_COLOR_SEQ}, - {(char*)"colqual", no_argument, 0, ARG_COLOR_QUAL}, + {(char*)"col-cseq", no_argument, 0, ARG_COLOR_SEQ}, + {(char*)"col-cqual", no_argument, 0, ARG_COLOR_QUAL}, + {(char*)"col-keepends", no_argument, 0, ARG_COLOR_KEEP_ENDS}, {(char*)"cost", no_argument, 0, ARG_COST}, {(char*)0, 0, 0, 0} // terminator }; @@ -459,11 +463,12 @@ static void printUsage(ostream& out) { << " --solexa1.3-quals input quals are from GA Pipeline ver. >= 1.3" << endl << " --integer-quals qualities are given as space-separated integers (not ASCII)" << endl << "Alignment:" << endl - << " -n/--seedmms max mismatches in seed (can be 0-3, default: -n 2)" << endl - << " -e/--maqerr max sum of mismatch quals (rounds like maq; default: 70)" << endl - << " -l/--seedlen seed length (default: 28)" << endl - << " --nomaqround disable Maq-like quality rounding (to nearest 10 <= 30)" << endl << " -v report end-to-end hits w/ <=v mismatches; ignore qualities" << endl + << " or" << endl + << " -n/--seedmms max mismatches in seed (can be 0-3, default: -n 2)" << endl + << " -e/--maqerr max sum of mismatch quals across alignment for -n (def: 70)" << endl + << " -l/--seedlen seed length for -n (default: 28)" << endl + << " --nomaqround disable Maq-like quality rounding for -n (nearest 10 <= 30)" << endl << " -I/--minins minimum insert size for paired-end alignment (default: 0)" << endl << " -X/--maxins maximum insert size for paired-end alignment (default: 250)" << endl << " --fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (default: --fr)" << endl @@ -481,8 +486,6 @@ static void printUsage(ostream& out) { << " --best hits guaranteed best stratum; ties broken by quality" << endl << " --strata hits in sub-optimal strata aren't reported (requires --best)" << endl //<< " --strandfix attempt to fix strand biases (def: on w/ --best, off w/o)" << endl - << " --snpphred penalty for a SNP when decoding colorspace (default: 30)" << endl - << " --snpfrac approx. fraction of SNP bases (e.g. 0.001); sets --snpphred" << endl << "Output:" << endl << " -t/--time print wall-clock time taken by search phases" << endl << " -B/--offbase leftmost ref offset = in bowtie output (default: 0)" << endl @@ -493,9 +496,14 @@ static void printUsage(ostream& out) { << " --un write unaligned reads/pairs to file(s) " << endl << " --max write reads/pairs over -m limit to file(s) " << endl << " --suppress suppresses given columns (comma-delim'ed) in default output" << endl - << " --colseq print aligned colorspace seqs as colors, not decoded bases" << endl - << " --colqual print original colorspace quals, not decoded quals" << endl << " --fullref write entire ref name (default: only up to 1st space)" << endl + << "Colorspace:" << endl + << " --snpphred Phred penalty for SNP when decoding colorspace (def: 30)" << endl + << " or" << endl + << " --snpfrac approx. fraction of SNP bases (e.g. 0.001); sets --snpphred" << endl + << " --col-cseq print aligned colorspace seqs as colors, not decoded bases" << endl + << " --col-cqual print original colorspace quals, not decoded quals" << endl + << " --col-keepends keep nucleotides at extreme ends of decoded alignment" << endl << "SAM:" << endl << " -S/--sam write hits in SAM format" << endl << " --mapq default mapping quality (MAPQ) to print for SAM alignments" << endl @@ -665,6 +673,7 @@ static void parseOptions(int argc, const char **argv) { case ARG_PHRED64: phred64Quals = true; break; case ARG_PHRED33: solexaQuals = false; phred64Quals = false; break; case ARG_NOMAQROUND: noMaqRound = true; break; + case ARG_COLOR_KEEP_ENDS: colorExEnds = false; break; case ARG_SNPPHRED: snpPhred = parseInt(0, "--snpphred must be at least 0"); break; case ARG_SNPFRAC: { double p = parse(optarg); @@ -1148,7 +1157,6 @@ static void *exactSearchWorker(void *vp) { GreedyDFSRangeSource bt( &ebwt, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh 0xffffffff, // max backtracks (no max) 0, // reportPartials (don't) @@ -1188,7 +1196,6 @@ static void *exactSearchWorkerStateful(void *vp) { UnpairedExactAlignerV1Factory alSEfact( ebwt, NULL, - snpPhred, !nofw, !norc, _sink, @@ -1209,7 +1216,6 @@ static void *exactSearchWorkerStateful(void *vp) { PairedExactAlignerV1Factory alPEfact( ebwt, NULL, - snpPhred, color, !nofw, !norc, @@ -1375,7 +1381,6 @@ static void *mismatchSearchWorkerFullStateful(void *vp) { Unpaired1mmAlignerV1Factory alSEfact( ebwtFw, &ebwtBw, - snpPhred, !nofw, !norc, _sink, @@ -1397,7 +1402,6 @@ static void *mismatchSearchWorkerFullStateful(void *vp) { ebwtFw, &ebwtBw, color, - snpPhred, !nofw, !norc, useV1, @@ -1467,7 +1471,6 @@ static void* mismatchSearchWorkerFull(void *vp){ GreedyDFSRangeSource bt( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh 0xffffffff, // max backtracks (no max) 0, // reportPartials (don't) @@ -1684,7 +1687,6 @@ static void *twoOrThreeMismatchSearchWorkerStateful(void *vp) { ebwtFw, &ebwtBw, two, - snpPhred, !nofw, !norc, _sink, @@ -1706,7 +1708,6 @@ static void *twoOrThreeMismatchSearchWorkerStateful(void *vp) { ebwtFw, &ebwtBw, color, - snpPhred, !nofw, !norc, useV1, @@ -1763,7 +1764,6 @@ static void* twoOrThreeMismatchSearchWorkerFull(void *vp) { GreedyDFSRangeSource btr1( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh // Do not impose maximums in 2/3-mismatch mode 0xffffffff, // max backtracks (no limit) @@ -1778,7 +1778,6 @@ static void* twoOrThreeMismatchSearchWorkerFull(void *vp) { GreedyDFSRangeSource bt2( &ebwtBw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh // Do not impose maximums in 2/3-mismatch mode 0xffffffff, // max backtracks (no limit) @@ -1793,7 +1792,6 @@ static void* twoOrThreeMismatchSearchWorkerFull(void *vp) { GreedyDFSRangeSource bt3( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh (none) // Do not impose maximums in 2/3-mismatch mode 0xffffffff, // max backtracks (no limit) @@ -1808,7 +1806,6 @@ static void* twoOrThreeMismatchSearchWorkerFull(void *vp) { GreedyDFSRangeSource bthh3( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) 0xffffffff, // qualThresh // Do not impose maximums in 2/3-mismatch mode 0xffffffff, // max backtracks (no limit) @@ -1954,7 +1951,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btf1( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -1968,7 +1964,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource bt1( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -1984,7 +1979,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btf2( &ebwtBw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (no) @@ -2000,7 +1994,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btr2( &ebwtBw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh (none) maxBtsBetter, // max backtracks seedMms, // report partials (up to seedMms mms) @@ -2016,7 +2009,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btf3( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh (none) maxBtsBetter, // max backtracks seedMms, // reportPartials (do) @@ -2033,7 +2025,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btr3( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -2049,7 +2040,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btr23( &ebwtFw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -2067,7 +2057,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btf4( &ebwtBw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -2083,7 +2072,6 @@ static void* seededQualSearchWorkerFull(void *vp) { GreedyDFSRangeSource btf24( &ebwtBw, params, refs, // reference sequence (for colorspace) - snpPhred, // phred probability of SNP (for colorspace) qualCutoff, // qualThresh maxBtsBetter, // max backtracks 0, // reportPartials (don't) @@ -2150,7 +2138,6 @@ static void* seededQualSearchWorkerFullStateful(void *vp) { seedMms, seedLen, qualCutoff, - snpPhred, maxBts, _sink, *sinkFact, @@ -2178,7 +2165,6 @@ static void* seededQualSearchWorkerFullStateful(void *vp) { seedMms, seedLen, qualCutoff, - snpPhred, maxBts, _sink, *sinkFact, diff --git a/ebwt_search_backtrack.h b/ebwt_search_backtrack.h index 81a3376..16ccb5c 100644 --- a/ebwt_search_backtrack.h +++ b/ebwt_search_backtrack.h @@ -9,6 +9,7 @@ #include "range.h" #include "range_source.h" #include "aligner_metrics.h" +#include "search_globals.h" /** * Class that coordinates quality- and quantity-aware backtracking over @@ -27,7 +28,6 @@ class GreedyDFSRangeSource { const Ebwt* ebwt, const EbwtSearchParams& params, const BitPairReference* refs, - int snpPhred, uint32_t qualThresh, /// max acceptable q-distance const int maxBts, /// maximum # backtracks allowed uint32_t reportPartials = 0, @@ -41,7 +41,6 @@ class GreedyDFSRangeSource { bool halfAndHalf = false, // hacky way of supporting separate revisitable regions bool maqPenalty = true) : _refs(refs), - _snpPhred(snpPhred), _qry(NULL), _qlen(0), _qual(NULL), @@ -1540,7 +1539,7 @@ class GreedyDFSRangeSource { // their indices into the query string; not in terms // of their offset from the 3' or 5' end. if(_ebwt->reportChaseOne((*_qry), _qual, _name, - _color, _snpPhred, _refs, _mms, + _color, colorExEnds, snpPhred, _refs, _mms, _refcs, stackDepth, ri, top, bot, _qlen, stratum, cost, _params)) { @@ -1691,7 +1690,6 @@ class GreedyDFSRangeSource { } const BitPairReference* _refs; // reference sequences (or NULL if not colorspace) - int _snpPhred; // phred probability of SNP String* _qry; // query (read) sequence size_t _qlen; // length of _qry String* _qual; // quality values for _qry diff --git a/search_globals.h b/search_globals.h new file mode 100644 index 0000000..26598a5 --- /dev/null +++ b/search_globals.h @@ -0,0 +1,14 @@ +/* + * search_globals.h + * + * Created on: Dec 5, 2009 + * Author: Ben Langmead + */ + +#ifndef SEARCH_GLOBALS_H_ +#define SEARCH_GLOBALS_H_ + +extern bool colorExEnds; +extern int snpPhred; + +#endif /* SEARCH_GLOBALS_H_ */