This is an automated email from the git hooks/post-receive script. tille pushed a commit to branch master in repository last-align.
commit e10f0c8552f21df3f5824c50e0cb9e159aa789d5 Author: Andreas Tille <[email protected]> Date: Wed May 18 18:06:58 2016 +0200 Imported Upstream version 737 --- ChangeLog.txt | 43 +++++++++- doc/last-tuning.html | 13 ++- doc/last-tuning.txt | 14 +++- doc/lastal.html | 151 ++++++++++++++++++++-------------- doc/lastal.txt | 137 ++++++++++++++++++------------- src/Alignment.cc | 137 +++++++++++++++++++++---------- src/Alignment.hh | 21 +++-- src/GreedyXdropAligner.cc | 203 ++++++++++++++++++++++++++++++++++++++++++++++ src/GreedyXdropAligner.hh | 69 ++++++++++++++++ src/LastEvaluer.cc | 45 +++++++--- src/LastEvaluer.hh | 3 +- src/LastalArguments.cc | 103 ++++++++++++++--------- src/LastalArguments.hh | 6 +- src/lastal.cc | 68 +++++++++++----- src/makefile | 25 +++--- src/version.hh | 2 +- 16 files changed, 789 insertions(+), 251 deletions(-) diff --git a/ChangeLog.txt b/ChangeLog.txt index 401615e..42d258b 100644 --- a/ChangeLog.txt +++ b/ChangeLog.txt @@ -1,10 +1,51 @@ +2016-03-31 Martin C. Frith <Martin C. Frith> + + * doc/lastal.txt, src/LastalArguments.cc: + Tried to document lastal options more clearly. + [9ecd0219342e] [tip] + +2016-03-30 Martin C. Frith <Martin C. Frith> + + * test/last-test.out, test/last-test.sh: + Added a test for minimum-difference alignment. + [3eebddccaf4e] + + * doc/last-tuning.txt, doc/lastal.txt, src/Alignment.cc, + src/Alignment.hh, src/GreedyXdropAligner.cc, + src/GreedyXdropAligner.hh, src/LastalArguments.cc, + src/LastalArguments.hh, src/lastal.cc, src/makefile, test/last- + test.out: + Added lastal option -M for fast-but-crude minimum-difference + alignment! + [a1fc3309828f] + + * doc/lastal.txt, src/LastEvaluer.cc, src/LastEvaluer.hh, + src/LastalArguments.cc, src/LastalArguments.hh, src/lastal.cc, test + /last-test.out, test/last-test.sh: + Added DNA-versus-protein alignment without frameshifts (-F0). + [cf14b8c39bf9] + +2016-03-29 Martin C. Frith <Martin C. Frith> + + * doc/lastal.txt, src/Alignment.cc, src/Alignment.hh, + src/LastalArguments.cc, src/lastal.cc: + ! Made lastal's lowercase/repeat masking faster & slightly + different. + [a282cc2f6937] + +2016-03-28 Martin C. Frith <Martin C. Frith> + + * doc/lastal.txt, src/Alignment.cc, src/Alignment.hh, src/lastal.cc: + Added warning about overlap alignment (-T1) to doc; refactoring. + [87f9c581e48a] + 2016-03-11 Martin C. Frith <Martin C. Frith> * doc/last-tuning.txt, doc/lastal.txt, doc/lastdb.txt, src/LastalArguments.cc, src/LastdbArguments.cc, src/lastal.cc, test /last-test.out, test/last-test.sh: Added minimizer doc & test, tweaked lastal verbose messages. - [323c152700e9] [tip] + [323c152700e9] 2016-03-10 Martin C. Frith <Martin C. Frith> diff --git a/doc/last-tuning.html b/doc/last-tuning.html index 047f87b..fc16735 100644 --- a/doc/last-tuning.html +++ b/doc/last-tuning.html @@ -415,11 +415,18 @@ halving it.</p> example, -C2 makes it discard alignments (before gapped extension) whose query coordinates lie in those of 2 or more stronger alignments.</p> </div> +<div class="section" id="id2"> +<h3>lastal -M</h3> +<p>This option requests "minimum-difference" alignment, which is <strong>faster +but cruder</strong> than standard gapped alignment. This treats all matches +the same, and minimizes the number of differences (mismatches plus +gaps).</p> +</div> <div class="section" id="lastal-j1"> <h3>lastal -j1</h3> -<p>This option requests <strong>gapless</strong> alignment, which is <strong>faster</strong>. (You -could get the same effect by using very high gap costs, but -j1 is -faster because it skips the gapping phase entirely.)</p> +<p>This option requests <strong>gapless</strong> alignment, which is even <strong>faster</strong>. +(You could get the same effect by using very high gap costs, but -j1 +is faster because it skips the gapping phase entirely.)</p> </div> <div class="section" id="lastal-f"> <h3>lastal -f</h3> diff --git a/doc/last-tuning.txt b/doc/last-tuning.txt index f3fc430..90af428 100644 --- a/doc/last-tuning.txt +++ b/doc/last-tuning.txt @@ -117,12 +117,20 @@ This option (gapless alignment culling) can make lastal **faster** but example, -C2 makes it discard alignments (before gapped extension) whose query coordinates lie in those of 2 or more stronger alignments. +lastal -M +--------- + +This option requests "minimum-difference" alignment, which is **faster +but cruder** than standard gapped alignment. This treats all matches +the same, and minimizes the number of differences (mismatches plus +gaps). + lastal -j1 ---------- -This option requests **gapless** alignment, which is **faster**. (You -could get the same effect by using very high gap costs, but -j1 is -faster because it skips the gapping phase entirely.) +This option requests **gapless** alignment, which is even **faster**. +(You could get the same effect by using very high gap costs, but -j1 +is faster because it skips the gapping phase entirely.) lastal -f --------- diff --git a/doc/lastal.html b/doc/lastal.html index 52bd570..319ba5d 100644 --- a/doc/lastal.html +++ b/doc/lastal.html @@ -336,17 +336,12 @@ zcat seqs.fasta.gz | lastal humanDb > myalns.maf <div class="section" id="steps-in-lastal"> <h2>Steps in lastal</h2> <ol class="arabic"> -<li><p class="first">Find similar regions</p> -</li> -</ol> -<blockquote> -<ol class="loweralpha"> <li><p class="first">Find initial matches. For each possible start position in the -query: find the shortest match with length ≥ l that <em>either</em> -occurs ≤ m times in the reference, <em>or</em> has length L.</p> +query: find the shortest match with length ≥ l that <em>either</em> occurs +≤ m times in the reference, <em>or</em> has length L.</p> </li> -<li><p class="first">Extend a gapless alignment from each initial match, and keep -those with score ≥ d.</p> +<li><p class="first">Extend a gapless alignment from each initial match, and keep those +with score ≥ d.</p> </li> <li><p class="first">Define cores: find the longest run of identical matches in each gapless alignment.</p> @@ -354,18 +349,6 @@ gapless alignment.</p> <li><p class="first">Extend a gapped alignment from either side of each core, and keep those with score ≥ e.</p> </li> -</ol> -</blockquote> -<ol class="arabic" start="2"> -<li><p class="first">Align them</p> -</li> -</ol> -<blockquote> -<ol class="loweralpha" start="5"> -<li><p class="first">"Final" alignments: extend gapped alignments again, but with -different masking and/or maximum score drop parameters. This -step is performed only if u=2 or if z differs from x.</p> -</li> <li><p class="first">Non-redundantize the alignments: if several alignments share an endpoint (same coordinates in both sequences), remove all but one highest-scoring one.</p> @@ -376,7 +359,6 @@ highest-scoring one.</p> gamma-centroid or LAMA (OFF by default).</p> </li> </ol> -</blockquote> </div> <div class="section" id="options"> <h2>Options</h2> @@ -493,7 +475,9 @@ standard affine gaps.</td></tr> <kbd><span class="option">-F <var>COST</var></span></kbd></td> <td><p class="first">Align DNA queries to protein reference sequences, using the specified frameshift cost. A value of 15 seems to be -reasonable. The output looks like this:</p> +reasonable. (As a special case, -F0 means DNA-versus-protein +alignment without frameshifts, which is faster.) The output +looks like this:</p> <pre class="literal-block"> a score=108 s prot 2 40 + 649 FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD @@ -519,7 +503,10 @@ regions in alignments) and speed (the smaller the faster).</td></tr> <td>Maximum score drop for gapless alignments.</td></tr> <tr><td class="option-group"> <kbd><span class="option">-z <var>DROP</var></span></kbd></td> -<td>Maximum score drop for final gapped alignments.</td></tr> +<td>Maximum score drop for final gapped alignments. Setting z +different from x causes lastal to extend gapless alignments +twice: first with a maximum score drop of x, and then with a +(presumably higher) maximum score drop of z.</td></tr> <tr><td class="option-group"> <kbd><span class="option">-d <var>SCORE</var></span></kbd></td> <td>Minimum score for gapless alignments.</td></tr> @@ -553,6 +540,44 @@ only affects the default value of -e, so if you specify -e then </table> </blockquote> </div> +<div class="section" id="initial-match-options"> +<h3>Initial-match options</h3> +<blockquote> +<table class="docutils option-list" frame="void" rules="none"> +<col class="option" /> +<col class="description" /> +<tbody valign="top"> +<tr><td class="option-group"> +<kbd><span class="option">-m <var>MULTIPLICITY</var></span></kbd></td> +<td><p class="first">Maximum multiplicity for initial matches. Each initial match is +lengthened until it occurs at most this many times in the +reference.</p> +<p class="last">If the reference was split into volumes by lastdb, then lastal +uses one volume at a time. The maximum multiplicity then +applies to each volume, not the whole reference. This is why +voluming changes the results.</p> +</td></tr> +<tr><td class="option-group"> +<kbd><span class="option">-l <var>LENGTH</var></span></kbd></td> +<td>Minimum length for initial matches. Length means the number of +letters spanned by the match.</td></tr> +<tr><td class="option-group"> +<kbd><span class="option">-L <var>LENGTH</var></span></kbd></td> +<td>Maximum length for initial matches.</td></tr> +<tr><td class="option-group"> +<kbd><span class="option">-k <var>STEP</var></span></kbd></td> +<td>Look for initial matches starting only at every STEP-th position +in the query. This makes lastal faster but less sensitive.</td></tr> +<tr><td class="option-group"> +<kbd><span class="option">-W <var>SIZE</var></span></kbd></td> +<td>Look for initial matches starting only at query positions that +are "minimum" in any window of SIZE consecutive positions (see +<a class="reference external" href="lastdb.html">lastdb.html</a>). By default, this parameter takes the same +value as was used for lastdb -W.</td></tr> +</tbody> +</table> +</blockquote> +</div> <div class="section" id="miscellaneous-options"> <h3>Miscellaneous options</h3> <blockquote> @@ -574,31 +599,38 @@ aligned to either strand of the query. 1 means that the matrix applies to either strand of the reference aligned to the forward strand of the query.</td></tr> <tr><td class="option-group"> +<kbd><span class="option">-M</span></kbd></td> +<td><p class="first">Find minimum-difference alignments, which is faster but cruder. +This treats all matches the same, and minimizes the number of +differences (mismatches plus gaps).</p> +<ul class="last"> +<li><p class="first">Any substitution score matrix will be ignored. The +substitution scores are set by the match score (r) and the +mismatch cost (q).</p> +</li> +<li><p class="first">The gap cost parameters will be ignored. The gap existence +cost will be 0 and the gap extension cost will be q + r/2.</p> +</li> +<li><p class="first">The match score (r) must be an even number.</p> +</li> +<li><p class="first">Any sequence quality data (e.g. fastq) will be ignored.</p> +</li> +</ul> +</td></tr> +<tr><td class="option-group"> <kbd><span class="option">-T <var>NUMBER</var></span></kbd></td> -<td>Type of alignment: 0 means "local alignment" and 1 means +<td><p class="first">Type of alignment: 0 means "local alignment" and 1 means "overlap alignment". Local alignments can end anywhere in the middle or at the ends of the sequences. Overlap alignments must extend to the left until they hit the end of a sequence (either query or reference), and to the right until they hit the end of -a sequence.</td></tr> -<tr><td class="option-group"> -<kbd><span class="option">-m <var>MULTIPLICITY</var></span></kbd></td> -<td><p class="first">Maximum multiplicity for initial matches. Each initial match is -lengthened until it occurs at most this many times in the -reference.</p> -<p class="last">If the reference was split into volumes by lastdb, then lastal -uses one volume at a time. The maximum multiplicity then -applies to each volume, not the whole reference. This is why -voluming changes the results.</p> +a sequence.</p> +<p class="last"><strong>Warning:</strong> it's often a bad idea to use -T1. This setting +does not change the maximum score drop allowed inside +alignments, so if an alignment cannot be extended to the end of +a sequence without exceeding this drop, it will be discarded.</p> </td></tr> <tr><td class="option-group"> -<kbd><span class="option">-l <var>LENGTH</var></span></kbd></td> -<td>Minimum length for initial matches. Length means the number of -letters spanned by the match.</td></tr> -<tr><td class="option-group"> -<kbd><span class="option">-L <var>LENGTH</var></span></kbd></td> -<td>Maximum length for initial matches.</td></tr> -<tr><td class="option-group"> <kbd><span class="option">-n <var>COUNT</var></span></kbd></td> <td>Maximum number of gapless alignments per query position. When lastal extends gapless alignments from initial matches that @@ -619,16 +651,6 @@ alignments with higher score (and on the same strand). This is a useful way to get just the top few hits to each part of each query (P Berman et al. 2000, J Comput Biol 7:293-302).</td></tr> <tr><td class="option-group"> -<kbd><span class="option">-k <var>STEP</var></span></kbd></td> -<td>Look for initial matches starting only at every STEP-th position -in the query. This makes lastal faster but less sensitive.</td></tr> -<tr><td class="option-group"> -<kbd><span class="option">-W <var>SIZE</var></span></kbd></td> -<td>Look for initial matches starting only at query positions that -are "minimum" in any window of SIZE consecutive positions (see -<a class="reference external" href="lastdb.html">lastdb.html</a>). By default, this parameter takes the same -value as was used for lastdb -W.</td></tr> -<tr><td class="option-group"> <kbd><span class="option">-i <var>BYTES</var></span></kbd></td> <td><p class="first">Search queries in batches of at most this many bytes. If a single sequence exceeds this amount, however, it is not split. @@ -671,13 +693,24 @@ applied to the DNA after translation, at the protein level.</p> </td></tr> <tr><td class="option-group"> <kbd><span class="option">-u <var>NUMBER</var></span></kbd></td> -<td>Specify treatment of lowercase letters when extending -alignments. 0 means do not mask them; 1 means mask them for -gapless extensions; 2 means mask them for gapless and gapped -extensions but not final extensions; 3 means mask them at all -stages. "Mask" means change their match/mismatch scores to -min(unmasked score, 0). This option does not affect treatment -of lowercase for initial matches.</td></tr> +<td><p class="first">Specify treatment of lowercase letters when extending +alignments:</p> +<ol class="arabic" start="0"> +<li><p class="first">Mask them for neither gapless nor gapped extensions.</p> +</li> +<li><p class="first">Mask them for gapless but not gapped extensions.</p> +</li> +<li><p class="first">Mask them for gapless but not gapped extensions, and then +discard alignments that lack any segment with score ≥ e when +lowercase is masked.</p> +</li> +<li><p class="first">Mask them for gapless and gapped extensions.</p> +</li> +</ol> +<p class="last">"Mask" means change their match/mismatch scores to min(unmasked +score, 0). This option does not affect treatment of lowercase +for initial matches.</p> +</td></tr> <tr><td class="option-group"> <kbd><span class="option">-w <var>DISTANCE</var></span></kbd></td> <td><p class="first">This option is a kludge to avoid catastrophic time and memory diff --git a/doc/lastal.txt b/doc/lastal.txt index 2d289e9..7e93663 100644 --- a/doc/lastal.txt +++ b/doc/lastal.txt @@ -20,35 +20,27 @@ You can also pipe query sequences into lastal, for example:: Steps in lastal --------------- -1) Find similar regions +1) Find initial matches. For each possible start position in the + query: find the shortest match with length ≥ l that *either* occurs + ≤ m times in the reference, *or* has length L. - a) Find initial matches. For each possible start position in the - query: find the shortest match with length ≥ l that *either* - occurs ≤ m times in the reference, *or* has length L. +2) Extend a gapless alignment from each initial match, and keep those + with score ≥ d. - b) Extend a gapless alignment from each initial match, and keep - those with score ≥ d. +3) Define cores: find the longest run of identical matches in each + gapless alignment. - c) Define cores: find the longest run of identical matches in each - gapless alignment. +4) Extend a gapped alignment from either side of each core, and keep + those with score ≥ e. - d) Extend a gapped alignment from either side of each core, and keep - those with score ≥ e. +5) Non-redundantize the alignments: if several alignments share an + endpoint (same coordinates in both sequences), remove all but one + highest-scoring one. -2) Align them +6) Estimate the ambiguity of each aligned column (OFF by default). - e) "Final" alignments: extend gapped alignments again, but with - different masking and/or maximum score drop parameters. This - step is performed only if u=2 or if z differs from x. - - f) Non-redundantize the alignments: if several alignments share an - endpoint (same coordinates in both sequences), remove all but one - highest-scoring one. - - g) Estimate the ambiguity of each aligned column (OFF by default). - - h) Redo the alignments to minimize column ambiguity, using either - gamma-centroid or LAMA (OFF by default). +7) Redo the alignments to minimize column ambiguity, using either + gamma-centroid or LAMA (OFF by default). Options ------- @@ -156,7 +148,9 @@ Score options -F COST Align DNA queries to protein reference sequences, using the specified frameshift cost. A value of 15 seems to be - reasonable. The output looks like this:: + reasonable. (As a special case, -F0 means DNA-versus-protein + alignment without frameshifts, which is faster.) The output + looks like this:: a score=108 s prot 2 40 + 649 FLLQAVKLQDP-STPHQIVPSP-VSDLIATHTLCPRMKYQDD @@ -181,7 +175,10 @@ Score options Maximum score drop for gapless alignments. -z DROP - Maximum score drop for final gapped alignments. + Maximum score drop for final gapped alignments. Setting z + different from x causes lastal to extend gapless alignments + twice: first with a maximum score drop of x, and then with a + (presumably higher) maximum score drop of z. -d SCORE Minimum score for gapless alignments. @@ -204,6 +201,36 @@ E-value options only affects the default value of -e, so if you specify -e then -E has no effect. +Initial-match options +~~~~~~~~~~~~~~~~~~~~~ + + -m MULTIPLICITY + Maximum multiplicity for initial matches. Each initial match is + lengthened until it occurs at most this many times in the + reference. + + If the reference was split into volumes by lastdb, then lastal + uses one volume at a time. The maximum multiplicity then + applies to each volume, not the whole reference. This is why + voluming changes the results. + + -l LENGTH + Minimum length for initial matches. Length means the number of + letters spanned by the match. + + -L LENGTH + Maximum length for initial matches. + + -k STEP + Look for initial matches starting only at every STEP-th position + in the query. This makes lastal faster but less sensitive. + + -W SIZE + Look for initial matches starting only at query positions that + are "minimum" in any window of SIZE consecutive positions (see + `<lastdb.html>`_). By default, this parameter takes the same + value as was used for lastdb -W. + Miscellaneous options ~~~~~~~~~~~~~~~~~~~~~ @@ -220,6 +247,18 @@ Miscellaneous options applies to either strand of the reference aligned to the forward strand of the query. + -M Find minimum-difference alignments, which is faster but cruder. + This treats all matches the same, and minimizes the number of + differences (mismatches plus gaps). + + * Any substitution score matrix will be ignored. The + substitution scores are set by the match score (r) and the + mismatch cost (q). + * The gap cost parameters will be ignored. The gap existence + cost will be 0 and the gap extension cost will be q + r/2. + * The match score (r) must be an even number. + * Any sequence quality data (e.g. fastq) will be ignored. + -T NUMBER Type of alignment: 0 means "local alignment" and 1 means "overlap alignment". Local alignments can end anywhere in the @@ -228,22 +267,10 @@ Miscellaneous options query or reference), and to the right until they hit the end of a sequence. - -m MULTIPLICITY - Maximum multiplicity for initial matches. Each initial match is - lengthened until it occurs at most this many times in the - reference. - - If the reference was split into volumes by lastdb, then lastal - uses one volume at a time. The maximum multiplicity then - applies to each volume, not the whole reference. This is why - voluming changes the results. - - -l LENGTH - Minimum length for initial matches. Length means the number of - letters spanned by the match. - - -L LENGTH - Maximum length for initial matches. + **Warning:** it's often a bad idea to use -T1. This setting + does not change the maximum score drop allowed inside + alignments, so if an alignment cannot be extended to the end of + a sequence without exceeding this drop, it will be discarded. -n COUNT Maximum number of gapless alignments per query position. When @@ -265,16 +292,6 @@ Miscellaneous options a useful way to get just the top few hits to each part of each query (P Berman et al. 2000, J Comput Biol 7:293-302). - -k STEP - Look for initial matches starting only at every STEP-th position - in the query. This makes lastal faster but less sensitive. - - -W SIZE - Look for initial matches starting only at query positions that - are "minimum" in any window of SIZE consecutive positions (see - `<lastdb.html>`_). By default, this parameter takes the same - value as was used for lastdb -W. - -i BYTES Search queries in batches of at most this many bytes. If a single sequence exceeds this amount, however, it is not split. @@ -313,12 +330,18 @@ Miscellaneous options -u NUMBER Specify treatment of lowercase letters when extending - alignments. 0 means do not mask them; 1 means mask them for - gapless extensions; 2 means mask them for gapless and gapped - extensions but not final extensions; 3 means mask them at all - stages. "Mask" means change their match/mismatch scores to - min(unmasked score, 0). This option does not affect treatment - of lowercase for initial matches. + alignments: + + 0. Mask them for neither gapless nor gapped extensions. + 1. Mask them for gapless but not gapped extensions. + 2. Mask them for gapless but not gapped extensions, and then + discard alignments that lack any segment with score ≥ e when + lowercase is masked. + 3. Mask them for gapless and gapped extensions. + + "Mask" means change their match/mismatch scores to min(unmasked + score, 0). This option does not affect treatment of lowercase + for initial matches. -w DISTANCE This option is a kludge to avoid catastrophic time and memory diff --git a/src/Alignment.cc b/src/Alignment.cc index 7fe6221..ec2db83 100644 --- a/src/Alignment.cc +++ b/src/Alignment.cc @@ -3,9 +3,9 @@ #include "Alignment.hh" #include "Alphabet.hh" #include "Centroid.hh" -#include "GappedXdropAligner.hh" #include "GeneticCode.hh" #include "GeneralizedAffineGapCosts.hh" +#include "GreedyXdropAligner.hh" #include "TwoQualityScoreMatrix.hh" #include <cassert> @@ -69,7 +69,8 @@ static bool isNext( const SegmentPair& x, const SegmentPair& y ){ return x.end1() == y.beg1() && x.end2() == y.beg2(); } -void Alignment::makeXdrop( GappedXdropAligner& aligner, Centroid& centroid, +void Alignment::makeXdrop( Centroid& centroid, + GreedyXdropAligner& greedyAligner, bool isGreedy, const uchar* seq1, const uchar* seq2, int globality, const ScoreMatrixRow* scoreMatrix, int smMax, const GeneralizedAffineGapCosts& gap, int maxDrop, @@ -96,8 +97,8 @@ void Alignment::makeXdrop( GappedXdropAligner& aligner, Centroid& centroid, // extend a gapped alignment in the left/reverse direction from the seed: std::vector<uchar>& columnAmbiguityCodes = extras.columnAmbiguityCodes; - extend( blocks, columnAmbiguityCodes, aligner, centroid, seq1, seq2, - seed.beg1(), seed.beg2(), false, globality, + extend( blocks, columnAmbiguityCodes, centroid, greedyAligner, isGreedy, + seq1, seq2, seed.beg1(), seed.beg2(), false, globality, scoreMatrix, smMax, maxDrop, gap, frameshiftCost, frameSize, pssm2, sm2qual, qual1, qual2, alph, extras, gamma, outputType ); @@ -116,8 +117,8 @@ void Alignment::makeXdrop( GappedXdropAligner& aligner, Centroid& centroid, // extend a gapped alignment in the right/forward direction from the seed: std::vector<SegmentPair> forwardBlocks; std::vector<uchar> forwardAmbiguities; - extend( forwardBlocks, forwardAmbiguities, aligner, centroid, seq1, seq2, - seed.end1(), seed.end2(), true, globality, + extend( forwardBlocks, forwardAmbiguities, centroid, greedyAligner, isGreedy, + seq1, seq2, seed.end1(), seed.end2(), true, globality, scoreMatrix, smMax, maxDrop, gap, frameshiftCost, frameSize, pssm2, sm2qual, qual1, qual2, alph, extras, gamma, outputType ); @@ -168,31 +169,34 @@ void Alignment::makeXdrop( GappedXdropAligner& aligner, Centroid& centroid, forwardAmbiguities.rend() ); } +// cost of the gap between x and y +static int gapCost(const SegmentPair &x, const SegmentPair &y, + const GeneralizedAffineGapCosts &gapCosts, + int frameshiftCost, size_t frameSize) { + size_t gapSize1 = y.beg1() - x.end1(); + size_t gapSize2, frameshift2; + sizeAndFrameshift(x.end2(), y.beg2(), frameSize, gapSize2, frameshift2); + int cost = gapCosts.cost(gapSize1, gapSize2); + if (frameshift2) cost += frameshiftCost; + return cost; +} + bool Alignment::isOptimal( const uchar* seq1, const uchar* seq2, int globality, const ScoreMatrixRow* scoreMatrix, int maxDrop, - const GeneralizedAffineGapCosts& gap, + const GeneralizedAffineGapCosts& gapCosts, int frameshiftCost, size_t frameSize, const ScoreMatrixRow* pssm2, const TwoQualityScoreMatrix& sm2qual, - const uchar* qual1, const uchar* qual2 ){ + const uchar* qual1, const uchar* qual2 ) const{ int maxScore = 0; int runningScore = 0; for( size_t i = 0; i < blocks.size(); ++i ){ const SegmentPair& y = blocks[i]; + if( i > 0 ){ // between each pair of aligned blocks: const SegmentPair& x = blocks[i - 1]; - size_t gapBeg1 = x.end1(); - size_t gapEnd1 = y.beg1(); - size_t gapSize1 = gapEnd1 - gapBeg1; - - size_t gapBeg2 = x.end2(); - size_t gapEnd2 = y.beg2(); - size_t gapSize2, frameshift2; - sizeAndFrameshift( gapBeg2, gapEnd2, frameSize, gapSize2, frameshift2 ); - if( frameshift2 ) runningScore -= frameshiftCost; - - runningScore -= gap.cost( gapSize1, gapSize2 ); + runningScore -= gapCost( x, y, gapCosts, frameshiftCost, frameSize ); if( !globality && runningScore <= 0 ) return false; if( runningScore < maxScore - maxDrop ) return false; } @@ -200,6 +204,7 @@ bool Alignment::isOptimal( const uchar* seq1, const uchar* seq2, int globality, const uchar* s1 = seq1 + y.beg1(); const uchar* s2 = seq2 + y.beg2(); const uchar* e1 = seq1 + y.end1(); + const ScoreMatrixRow* p2 = pssm2 ? pssm2 + y.beg2() : 0; const uchar* q1 = qual1 ? qual1 + y.beg1() : 0; const uchar* q2 = qual2 ? qual2 + y.beg2() : 0; @@ -219,9 +224,48 @@ bool Alignment::isOptimal( const uchar* seq1, const uchar* seq2, int globality, return true; } +bool Alignment::hasGoodSegment(const uchar *seq1, const uchar *seq2, + int minScore, const ScoreMatrixRow *scoreMatrix, + const GeneralizedAffineGapCosts &gapCosts, + int frameshiftCost, size_t frameSize, + const ScoreMatrixRow *pssm2, + const TwoQualityScoreMatrix &sm2qual, + const uchar *qual1, const uchar *qual2) const { + int score = 0; + + for (size_t i = 0; i < blocks.size(); ++i) { + const SegmentPair& y = blocks[i]; + + if (i > 0) { // between each pair of aligned blocks: + score -= gapCost(blocks[i - 1], y, gapCosts, frameshiftCost, frameSize); + if (score < 0) score = 0; + } + + const uchar *s1 = seq1 + y.beg1(); + const uchar *s2 = seq2 + y.beg2(); + const uchar *e1 = seq1 + y.end1(); + + const ScoreMatrixRow *p2 = pssm2 ? pssm2 + y.beg2() : 0; + const uchar *q1 = qual1 ? qual1 + y.beg1() : 0; + const uchar *q2 = qual2 ? qual2 + y.beg2() : 0; + + while (s1 < e1) { + /**/ if (sm2qual) score += sm2qual(*s1++, *s2++, *q1++, *q2++); + else if (pssm2) score += (*p2++)[*s1++]; + else score += scoreMatrix[*s1++][*s2++]; + + if (score >= minScore) return true; + if (score < 0) score = 0; + } + } + + return false; +} + void Alignment::extend( std::vector< SegmentPair >& chunks, std::vector< uchar >& ambiguityCodes, - GappedXdropAligner& aligner, Centroid& centroid, + Centroid& centroid, + GreedyXdropAligner& greedyAligner, bool isGreedy, const uchar* seq1, const uchar* seq2, size_t start1, size_t start2, bool isForward, int globality, @@ -233,8 +277,11 @@ void Alignment::extend( std::vector< SegmentPair >& chunks, const uchar* qual1, const uchar* qual2, const Alphabet& alph, AlignmentExtras& extras, double gamma, int outputType ){ + GappedXdropAligner& aligner = centroid.aligner(); + if( frameSize ){ assert( outputType < 4 ); + assert( !isGreedy ); assert( !globality ); assert( !pssm2 ); assert( !sm2qual ); @@ -262,22 +309,24 @@ void Alignment::extend( std::vector< SegmentPair >& chunks, } int extensionScore = - sm2qual ? aligner.align2qual( seq1 + start1, qual1 + start1, - seq2 + start2, qual2 + start2, - isForward, globality, sm2qual, - gap.delExist, gap.delExtend, - gap.insExist, gap.insExtend, - gap.pairExtend, maxDrop, smMax ) - : pssm2 ? aligner.alignPssm( seq1 + start1, pssm2 + start2, - isForward, globality, - gap.delExist, gap.delExtend, - gap.insExist, gap.insExtend, - gap.pairExtend, maxDrop, smMax ) - : aligner.align( seq1 + start1, seq2 + start2, - isForward, globality, sm, - gap.delExist, gap.delExtend, - gap.insExist, gap.insExtend, - gap.pairExtend, maxDrop, smMax ); + isGreedy ? greedyAligner.align( seq1 + start1, seq2 + start2, + isForward, sm, maxDrop, alph.size ) + : sm2qual ? aligner.align2qual( seq1 + start1, qual1 + start1, + seq2 + start2, qual2 + start2, + isForward, globality, sm2qual, + gap.delExist, gap.delExtend, + gap.insExist, gap.insExtend, + gap.pairExtend, maxDrop, smMax ) + : pssm2 ? aligner.alignPssm( seq1 + start1, pssm2 + start2, + isForward, globality, + gap.delExist, gap.delExtend, + gap.insExist, gap.insExtend, + gap.pairExtend, maxDrop, smMax ) + : aligner.align( seq1 + start1, seq2 + start2, + isForward, globality, sm, + gap.delExist, gap.delExtend, + gap.insExist, gap.insExtend, + gap.pairExtend, maxDrop, smMax ); if( extensionScore == -INF ){ score = -INF; // avoid score overflow @@ -288,14 +337,20 @@ void Alignment::extend( std::vector< SegmentPair >& chunks, if( outputType < 5 || outputType == 7 ){ // ordinary max-score alignment size_t end1, end2, size; - while( aligner.getNextChunk( end1, end2, size, - gap.delExist, gap.delExtend, - gap.insExist, gap.insExtend, - gap.pairExtend ) ) - chunks.push_back( SegmentPair( end1 - size, end2 - size, size ) ); + if( isGreedy ){ + while( greedyAligner.getNextChunk( end1, end2, size ) ) + chunks.push_back( SegmentPair( end1 - size, end2 - size, size ) ); + }else{ + while( aligner.getNextChunk( end1, end2, size, + gap.delExist, gap.delExtend, + gap.insExist, gap.insExtend, + gap.pairExtend ) ) + chunks.push_back( SegmentPair( end1 - size, end2 - size, size ) ); + } } if( outputType > 3 ){ // calculate match probabilities + assert( !isGreedy ); assert( !sm2qual ); centroid.reset(); centroid.forward( seq1, seq2, start1, start2, isForward, globality, gap ); diff --git a/src/Alignment.hh b/src/Alignment.hh index 2ebf602..4d1e9fc 100644 --- a/src/Alignment.hh +++ b/src/Alignment.hh @@ -15,8 +15,8 @@ namespace cbrc{ typedef unsigned char uchar; -class GappedXdropAligner; class GeneralizedAffineGapCosts; +class GreedyXdropAligner; class LastEvaluer; class MultiSequence; class Alphabet; @@ -77,7 +77,8 @@ struct Alignment{ // Alignment might not be "optimal" (see below). // If outputType > 3: calculates match probabilities. // If outputType > 4: does gamma-centroid alignment. - void makeXdrop( GappedXdropAligner& aligner, Centroid& centroid, + void makeXdrop( Centroid& centroid, + GreedyXdropAligner& greedyAligner, bool isGreedy, const uchar* seq1, const uchar* seq2, int globality, const ScoreMatrixRow* scoreMatrix, int smMax, const GeneralizedAffineGapCosts& gap, int maxDrop, @@ -94,11 +95,20 @@ struct Alignment{ // If "globality" is non-zero, skip the prefix and suffix checks. bool isOptimal( const uchar* seq1, const uchar* seq2, int globality, const ScoreMatrixRow* scoreMatrix, int maxDrop, - const GeneralizedAffineGapCosts& gap, + const GeneralizedAffineGapCosts& gapCosts, int frameshiftCost, size_t frameSize, const ScoreMatrixRow* pssm2, const TwoQualityScoreMatrix& sm2qual, - const uchar* qual1, const uchar* qual2 ); + const uchar* qual1, const uchar* qual2 ) const; + + // Does the Alignment have any segment with score >= minScore? + bool hasGoodSegment(const uchar *seq1, const uchar *seq2, + int minScore, const ScoreMatrixRow *scoreMatrix, + const GeneralizedAffineGapCosts &gapCosts, + int frameshiftCost, size_t frameSize, + const ScoreMatrixRow *pssm2, + const TwoQualityScoreMatrix &sm2qual, + const uchar *qual1, const uchar *qual2) const; AlignmentText write(const MultiSequence& seq1, const MultiSequence& seq2, size_t seqNum2, char strand, const uchar* seqData2, @@ -118,7 +128,8 @@ struct Alignment{ void extend( std::vector< SegmentPair >& chunks, std::vector< uchar >& ambiguityCodes, - GappedXdropAligner& aligner, Centroid& centroid, + Centroid& centroid, + GreedyXdropAligner& greedyAligner, bool isGreedy, const uchar* seq1, const uchar* seq2, size_t start1, size_t start2, bool isForward, int globality, diff --git a/src/GreedyXdropAligner.cc b/src/GreedyXdropAligner.cc new file mode 100644 index 0000000..1606038 --- /dev/null +++ b/src/GreedyXdropAligner.cc @@ -0,0 +1,203 @@ +// Copyright 2016 Martin C. Frith + +#include "GreedyXdropAligner.hh" +#include <algorithm> +#include <cassert> +//#include <iostream> // for debugging + +static int maxValue(int a, int b, int c) { + return std::max(std::max(a, b), c); +} + +static void resizeIfSmaller(std::vector<int> &v, size_t s) { + if (v.size() < s) v.resize(s); +} + +const int undefined = -9; + +static const int *definedBeg(const int *beg, const int *end) { + while (beg < end && *beg == undefined) ++beg; + return beg; +} + +static const int *definedEnd(const int *beg, const int *end) { + while (end > beg && *(end - 1) == undefined) --end; + return end; +} + +namespace cbrc { + +// i: number of seq1 letters aligned so far +// j: number of seq2 letters aligned so far +// diagonal: i - j +// antidiagonal: i + j +// distance: number of differences (mismatches + gap characters) + +// furthest: a flattened, ragged 2D array, which holds the +// furthest antidiagonal reachable on a given diagonal +// with a given distance + +// score0: score + differences * (matchScore - mismatchScore), +// or equivalently, antidiagonal * matchScore / 2 + +int GreedyXdropAligner::align(const uchar *seq1, + const uchar *seq2, + bool isForward, + const ScoreMatrixRow *scorer, + int maxScoreDrop, + uchar delimiter) { + const int matchScore = scorer[0][0]; + const int mismatchScore = scorer[0][1]; + assert(matchScore % 2 == 0); + const int halfMatchScore = matchScore / 2; + const int differenceCost = matchScore - mismatchScore; + const int lookBack = (maxScoreDrop + halfMatchScore) / differenceCost + 1; + const int minScore0gain = lookBack * differenceCost - maxScoreDrop; + + if (!isForward) { + --seq1; + --seq2; + } + + minScore0s.assign(lookBack, 0); + + rowOrigins.resize(1); + size_t furthestSize = 2; + resizeIfSmaller(furthest, furthestSize); + furthest[0] = -1; // xxx tricky initialization + furthest[1] = undefined; // add one pad cell + + int lowerStop = INT_MIN; + int upperStop = INT_MAX; + + int diagonalBeg = 0; + int diagonalEnd = 1; + + bestDistance = -1; + int bestScore0 = 0; + + int distance; + for (distance = 0; ; ++distance) { + int minScore0 = minScore0s[distance]; + + size_t rowOrigin = furthestSize - diagonalBeg; + rowOrigins.push_back(rowOrigin); + size_t furthestSizeNew = rowOrigin + diagonalEnd + 2; // + 2 pad cells + resizeIfSmaller(furthest, furthestSizeNew); + + int *to = &furthest[furthestSize]; + const int *from = sources(distance, diagonalBeg); + + *to++ = undefined; // add one pad cell + const int *toBeg = to; + + for (int diagonal = diagonalBeg; diagonal < diagonalEnd; ++diagonal) { + int antidiagonal = maxValue(from[0], from[1] + 1, from[2]) + 1; + ++from; + if (halfMatchScore * antidiagonal >= minScore0) { + int i = (antidiagonal + diagonal) / 2; + int j = (antidiagonal - diagonal) / 2; + const uchar *s1; + const uchar *s2; + if (isForward) { + s1 = seq1 + i; + s2 = seq2 + j; + while (scorer[*s1][*s2] > 0) { ++s1; ++s2; } // skip past matches + i = s1 - seq1; + j = s2 - seq2; + } else { + s1 = seq1 - i; + s2 = seq2 - j; + while (scorer[*s1][*s2] > 0) { --s1; --s2; } // skip past matches + i = seq1 - s1; + j = seq2 - s2; + } + if (*s2 == delimiter) lowerStop = diagonal; + if (*s1 == delimiter && upperStop > diagonal) upperStop = diagonal; + antidiagonal = i + j; + *to = antidiagonal; + int score0 = halfMatchScore * antidiagonal; + if (score0 > bestScore0) { + bestScore0 = score0; + bestDistance = distance; + bestDiagonal = diagonal; + } + } else { + *to = undefined; + } + ++to; + } + + lowerStop += 2; // unnecessary, but might improve speed + upperStop -= 2; // unnecessary, but might improve speed + + diagonalBeg += definedBeg(toBeg, to) - toBeg - 1; + diagonalBeg = std::max(diagonalBeg, lowerStop + 1); + + diagonalEnd += definedEnd(toBeg, to) - to + 1; + diagonalEnd = std::min(diagonalEnd, upperStop); + + if (diagonalBeg >= diagonalEnd) break; + + *to = undefined; // add one pad cell + minScore0s.push_back(bestScore0 + minScore0gain); + bestScore0 += differenceCost; + furthestSize = furthestSizeNew; + } + + return bestScore0 - differenceCost * distance; +} + +bool GreedyXdropAligner::getNextChunk(size_t &end1, + size_t &end2, + size_t &length) { + if (bestDistance < 0) return false; + + int antidiagonal = furthest[rowOrigins[bestDistance + 1] + bestDiagonal + 1]; + end1 = (antidiagonal + bestDiagonal) / 2; + end2 = (antidiagonal - bestDiagonal) / 2; + + // skip back past substitutions, until we hit an indel or the start + // of the extension: + int del, mis, ins; + do { + const int *from = sources(bestDistance, bestDiagonal); + del = from[0]; + mis = from[1] + 1; + ins = from[2]; + bestDistance -= 1; + } while (mis >= del && mis >= ins); + + int newAntidiagonal; + if (del >= ins) { + bestDiagonal -= 1; + newAntidiagonal = del; + } else { + bestDiagonal += 1; + newAntidiagonal = ins; + } + length = (antidiagonal - newAntidiagonal) / 2; // odd number / 2 + + // skip back past indels, until we hit a gapless chunk or the start + // of the extension: + while (bestDistance >= 0) { + const int *from = sources(bestDistance, bestDiagonal); + del = from[0]; + mis = from[1] + 1; + ins = from[2]; + if (mis >= del && mis >= ins) break; + --newAntidiagonal; + if (del >= ins) { + if (del < newAntidiagonal) break; + bestDiagonal -= 1; + } else { + if (ins < newAntidiagonal) break; + bestDiagonal += 1; + } + bestDistance -= 1; + } + + return true; +} + +} diff --git a/src/GreedyXdropAligner.hh b/src/GreedyXdropAligner.hh new file mode 100644 index 0000000..93f734f --- /dev/null +++ b/src/GreedyXdropAligner.hh @@ -0,0 +1,69 @@ +// Copyright 2016 Martin C. Frith + +// These routines extend an alignment in a given direction (forward or +// reverse) from given start points in two sequences. + +// The algorithm is adapted from Section 3 of "A greedy algorithm for +// aligning DNA sequences" by Z Zhang, S Schwartz, L Wagner, W Miller, +// J Comput Biol. 2000 7(1-2):203-214. + +// Forward alignments begin at the given start points, whereas reverse +// alignments begin one-before the given start points. + +// To use: first call "align", which calculates the alignment but only +// returns its score. To get the actual alignment, call +// "getNextChunk" to get each gapless chunk. + +// The sequences had better end with delimiter characters. + +// The "scorer" indicates which letter pairs are considered matches. +// The algorithm uses a match score (taken from scorer[0][0]) and a +// mismatch cost (taken from scorer[0][1]). + +// The match score must be divisible by 2. + +#ifndef GREEDY_XDROP_ALIGNER_HH +#define GREEDY_XDROP_ALIGNER_HH + +#include "ScoreMatrixRow.hh" + +#include <stddef.h> +#include <vector> + +namespace cbrc { + +typedef unsigned char uchar; + +class GreedyXdropAligner { +public: + int align(const uchar *seq1, // start point in the 1st sequence + const uchar *seq2, // start point in the 2nd sequence + bool isForward, // forward or reverse extension? + const ScoreMatrixRow *scorer, // the substitution score matrix + int maxScoreDrop, + uchar delimiter); + + // Call this repeatedly to get each gapless chunk of the alignment. + // The chunks are returned in far-to-near order. The chunk's end + // coordinates in each sequence (relative to the start of extension) + // and length are returned in the 3 out-parameters. If there are no + // more chunks, the 3 parameters are unchanged and "false" is + // returned. + bool getNextChunk(size_t &end1, size_t &end2, size_t &length); + +private: + std::vector<size_t> rowOrigins; + std::vector<int> furthest; // analogous to the R array in Zhang et al. + std::vector<int> minScore0s; // analogous to the T array in Zhang et al. + + // Our position during the trace-back: + int bestDistance; + int bestDiagonal; + + const int *sources(int distance, int diagonal) const + { return &furthest[rowOrigins[distance] + diagonal]; } +}; + +} + +#endif diff --git a/src/LastEvaluer.cc b/src/LastEvaluer.cc index 0384b95..2c7ac4c 100644 --- a/src/LastEvaluer.cc +++ b/src/LastEvaluer.cc @@ -241,6 +241,14 @@ static void shrinkCodonTable(long *codonTable, const size_t *usedLettersBeg, usedLettersEnd, codonTable[i]) - usedLettersBeg; } +static void makeTxFreqs(double *txFreqs, const double *ntFreqs, + const long *codonTable) { + for (int i = 0; i < 4; ++i) + for (int j = 0; j < 4; ++j) + for (int k = 0; k < 4; ++k) + txFreqs[*codonTable++] += ntFreqs[i] * ntFreqs[j] * ntFreqs[k]; +} + void LastEvaluer::init(const char *matrixName, int matchScore, int mismatchCost, @@ -264,7 +272,7 @@ void LastEvaluer::init(const char *matrixName, const double maxSeconds = 60.0; // seems to work nicely on my computer size_t alphabetSize = std::strlen(alphabet); - if (frameshiftCost > 0) { // DNA-versus-protein alignment with frameshifts: + if (frameshiftCost >= 0) { // DNA-versus-protein alignment: if (isGapped && insOpen == delOpen && insEpen == delEpen) { if (isProtein(alphabet) && isStandardGeneticCode) { for (size_t i = 0; i < COUNTOF(frameshiftEvalueParameters); ++i) { @@ -290,16 +298,31 @@ void LastEvaluer::init(const char *matrixName, std::vector<double> aaFreqs(matrixSize); copy(letterFreqs1, letterFreqs1 + alphabetSize, aaFreqs.begin()); - if (isGapped) { - frameshiftEvaluer.initFrameshift(4, matrixSize, codonTable, - &matrix[0], &ntFreqs[0], &aaFreqs[0], - delOpen, delEpen, insOpen, insEpen, - frameshiftCost, true, - lambdaTolerance, kTolerance, - 0, maxMegabytes, randomSeed); - } else { - frameshiftEvaluer.initGapless(4, matrixSize, codonTable, - &matrix[0], &ntFreqs[0], &aaFreqs[0]); + if (frameshiftCost > 0) { // with frameshifts: + if (isGapped) { + frameshiftEvaluer.initFrameshift(4, matrixSize, codonTable, + &matrix[0], &ntFreqs[0], &aaFreqs[0], + delOpen, delEpen, insOpen, insEpen, + frameshiftCost, true, + lambdaTolerance, kTolerance, + 0, maxMegabytes, randomSeed); + } else { + frameshiftEvaluer.initGapless(4, matrixSize, codonTable, + &matrix[0], &ntFreqs[0], &aaFreqs[0]); + } + } else { // without frameshifts: + std::vector<double> txFreqs(matrixSize); + makeTxFreqs(&txFreqs[0], &ntFreqs[0], codonTable); + + if (isGapped) { + evaluer.set_gapped_computation_parameters_simplified(maxSeconds); + evaluer.initGapped(matrixSize, &matrix[0], &txFreqs[0], &aaFreqs[0], + delOpen, delEpen, insOpen, insEpen, + true, lambdaTolerance, kTolerance, + 0, maxMegabytes, randomSeed); + } else { + evaluer.initGapless(matrixSize, &matrix[0], &txFreqs[0], &aaFreqs[0]); + } } } else { // ordinary alignment: if (isGapped && insOpen == delOpen && insEpen == delEpen) { diff --git a/src/LastEvaluer.hh b/src/LastEvaluer.hh index d8d0ccd..eee084e 100644 --- a/src/LastEvaluer.hh +++ b/src/LastEvaluer.hh @@ -33,7 +33,8 @@ public: // "bad" state and throw an Sls::error. // These arguments are only used to lookup pre-calculated cases: // matrixName, matchScore, mismatchCost, isStandardGeneticCode. - // DNA-versus-protein alignment is indicated by: frameshiftCost > 0. + // DNA-versus-protein alignment is indicated by: frameshiftCost >= 0. + // As a special case, frameshiftCost==0 means no frameshifts. // For DNA-versus-protein alignment, letterFreqs2 is not used. void init(const char *matrixName, int matchScore, diff --git a/src/LastalArguments.cc b/src/LastalArguments.cc index 375c2f7..6a568fa 100644 --- a/src/LastalArguments.cc +++ b/src/LastalArguments.cc @@ -47,6 +47,7 @@ LastalArguments::LastalArguments() : outputType(3), strand(-1), // depends on the alphabet isQueryStrandMatrix(false), + isGreedy(false), globality(0), isKeepLowercase(true), // depends on the option used with lastdb tantanSetting(-1), // depends on the option used with lastdb @@ -93,8 +94,8 @@ void LastalArguments::fromArgs( int argc, char** argv, bool optionsOnly ){ Find local sequence alignments.\n\ \n\ Score options (default settings):\n\ --r: match score (DNA: 1, 0<Q<5: 6)\n\ --q: mismatch cost (DNA: 1, 0<Q<5: 18)\n\ +-r: match score (2 if -M, else 6 if 0<Q<5, else 1 if DNA)\n\ +-q: mismatch cost (3 if -M, else 18 if 0<Q<5, else 1 if DNA)\n\ -p: match/mismatch score matrix (protein-protein: BL62, DNA-protein: BL80)\n\ -a: gap existence cost (DNA: 7, protein: 11, 0<Q<5: 21)\n\ -b: gap extension cost (DNA: 1, protein: 2, 0<Q<5: 9)\n\ @@ -119,29 +120,32 @@ Cosmetic options (default settings):\n\ -v: be verbose: write messages about what lastal is doing\n\ -f: output format: TAB, MAF, BlastTab, BlastTab+ (MAF)\n\ \n\ +Initial-match options (default settings):\n\ +-m: maximum initial matches per query position (" + + stringify(oneHitMultiplicity) + ")\n\ +-l: minimum length for initial matches (" + + stringify(minHitDepth) + ")\n\ +-L: maximum length for initial matches (infinity)\n\ +-k: use initial matches starting at every k-th position in each query (" + + stringify(queryStep) + ")\n\ +-W: use \"minimum\" positions in sliding windows of W consecutive positions\n\ +\n\ Miscellaneous options (default settings):\n\ -s: strand: 0=reverse, 1=forward, 2=both (2 for DNA, 1 for protein)\n\ -S: score matrix applies to forward strand of: 0=reference, 1=query (" + stringify(isQueryStrandMatrix) + ")\n\ +-M: find minimum-difference alignments (faster but cruder)\n\ -T: type of alignment: 0=local, 1=overlap (" + stringify(globality) + ")\n\ --m: maximum initial matches per query position (" - + stringify(oneHitMultiplicity) + ")\n\ --l: minimum length for initial matches (" - + stringify(minHitDepth) + ")\n\ --L: maximum length for initial matches (infinity)\n\ -n: maximum gapless alignments per query position (infinity if m=0, else m)\n\ -C: omit gapless alignments in >= C others with > score-per-length (off)\n\ -K: omit alignments whose query range lies in >= K others with > score (off)\n\ --k: use initial matches starting at every k-th position in each query (" - + stringify(queryStep) + ")\n\ --W: use \"minimum\" positions in sliding windows of W consecutive positions\n\ -i: query batch size (8 KiB, unless there is > 1 thread or lastdb volume)\n\ -P: number of parallel threads (" + stringify(numOfThreads) + ")\n\ -R: repeat-marking options (the same as was used for lastdb)\n\ -u: mask lowercase during extensions: 0=never, 1=gapless,\n\ - 2=gapless+gapped but not final, 3=always (2 if lastdb -c and Q<5, else 0)\n\ + 2=gapless+postmask, 3=always (2 if lastdb -c and Q<5, else 0)\n\ -w: suppress repeats inside exact matches, offset by <= this distance (" + stringify(maxRepeatDistance) + ")\n\ -G: genetic code file\n\ @@ -163,7 +167,7 @@ LAST home page: http://last.cbrc.jp/\n\ optind = 1; // allows us to scan arguments more than once(???) int c; const char optionString[] = "hVvf:" "r:q:p:a:b:A:B:c:F:x:y:z:d:e:" "D:E:" - "s:S:T:m:l:L:n:C:K:k:W:i:P:R:u:w:t:g:G:j:Q:"; + "s:S:MT:m:l:L:n:C:K:k:W:i:P:R:u:w:t:g:G:j:Q:"; while( (c = myGetopt(argc, argv, optionString)) != -1 ){ switch(c){ case 'h': @@ -215,7 +219,7 @@ LAST home page: http://last.cbrc.jp/\n\ break; case 'F': unstringify( frameshiftCost, optarg ); - if( frameshiftCost <= 0 ) badopt( c, optarg ); + if( frameshiftCost < 0 ) badopt( c, optarg ); break; case 'x': unstringify( maxDropGapped, optarg ); @@ -247,6 +251,23 @@ LAST home page: http://last.cbrc.jp/\n\ if( maxEvalue <= 0 ) badopt( c, optarg ); break; + case 'm': + unstringify( oneHitMultiplicity, optarg ); + break; + case 'l': + unstringify( minHitDepth, optarg ); + break; + case 'L': + unstringify( maxHitDepth, optarg ); + break; + case 'k': + unstringify( queryStep, optarg ); + if( queryStep <= 0 ) badopt( c, optarg ); + break; + case 'W': + unstringify( minimizerWindow, optarg ); // allow 0, meaning "default" + break; + case 's': unstringify( strand, optarg ); if( strand < 0 || strand > 2 ) badopt( c, optarg ); @@ -254,19 +275,13 @@ LAST home page: http://last.cbrc.jp/\n\ case 'S': unstringify( isQueryStrandMatrix, optarg ); break; + case 'M': + isGreedy = true; + break; case 'T': unstringify( globality, optarg ); if( globality < 0 || globality > 1 ) badopt( c, optarg ); break; - case 'm': - unstringify( oneHitMultiplicity, optarg ); - break; - case 'l': - unstringify( minHitDepth, optarg ); - break; - case 'L': - unstringify( maxHitDepth, optarg ); - break; case 'n': unstringify( maxGaplessAlignmentsPerQueryPosition, optarg ); if( maxGaplessAlignmentsPerQueryPosition <= 0 ) badopt( c, optarg ); @@ -277,13 +292,6 @@ LAST home page: http://last.cbrc.jp/\n\ case 'K': unstringify( cullingLimitForFinalAlignments, optarg ); break; - case 'k': - unstringify( queryStep, optarg ); - if( queryStep <= 0 ) badopt( c, optarg ); - break; - case 'W': - unstringify( minimizerWindow, optarg ); // allow 0, meaning "default" - break; case 'i': unstringifySize( batchSize, optarg ); if( batchSize <= 0 ) badopt( c, optarg ); // 0 means "not specified" @@ -338,11 +346,11 @@ LAST home page: http://last.cbrc.jp/\n\ if( isTranslated() && inputFormat == 5 ) ERR( "can't combine option -F with option -Q 5" ); - if( isTranslated() && outputType > 3 ) - ERR( "can't combine option -F with option -j > 3" ); + if( isFrameshift() && outputType > 3 ) + ERR( "can't combine option -F > 0 with option -j > 3" ); - if( isTranslated() && globality == 1 ) - ERR( "can't combine option -F with option -T 1" ); + if( isFrameshift() && globality == 1 ) + ERR( "can't combine option -F > 0 with option -T 1" ); if( isTranslated() && isQueryStrandMatrix ) ERR( "can't combine option -F with option -S 1" ); @@ -350,6 +358,15 @@ LAST home page: http://last.cbrc.jp/\n\ if( globality == 1 && outputType == 1 ) ERR( "can't combine option -T 1 with option -j 1" ); + if( isGreedy && outputType > 3 ) + ERR( "can't combine option -M with option -j > 3" ); + + if( isGreedy && globality == 1 ) + ERR( "can't combine option -M with option -T 1" ); + + if( isGreedy && maskLowercase == 3 ) + ERR( "can't combine option -M with option -u 3" ); + if( optionsOnly ) return; if( optind >= argc ) ERR( "please give me a database name and sequence file(s)\n\n" + usage ); @@ -396,7 +413,18 @@ void LastalArguments::setDefaultsFromAlphabet( bool isDna, bool isProtein, unsigned realNumOfThreads ){ if( strand < 0 ) strand = (isDna || isTranslated()) ? 2 : 1; - if( isProtein ){ + if( isGreedy ){ + if( matchScore < 0 ) matchScore = 2; + if( mismatchCost < 0 ) mismatchCost = 3; + gapExistCost = 0; + gapExtendCost = mismatchCost + matchScore / 2; + insExistCost = gapExistCost; + insExtendCost = gapExtendCost; + if( frameshiftCost > 0 ) frameshiftCost = 0; + if( matchScore % 2 ) + ERR( "with option -M, the match score (option -r) must be even" ); + } + else if( isProtein ){ // default match & mismatch scores: Blosum62 matrix if( matchScore < 0 && mismatchCost >= 0 ) matchScore = 1; // idiot-proof if( mismatchCost < 0 && matchScore >= 0 ) mismatchCost = 1; // idiot-proof @@ -467,12 +495,12 @@ void LastalArguments::setDefaultsFromAlphabet( bool isDna, bool isProtein, if( minimizerWindow == 0 ) minimizerWindow = refMinimizerWindow; - if( isTranslated() && frameshiftCost < gapExtendCost ) + if( isFrameshift() && frameshiftCost < gapExtendCost ) ERR( "the frameshift cost must not be less than the gap extension cost" ); if( insExistCost != gapExistCost || insExtendCost != gapExtendCost ){ - if( isTranslated() ) - ERR( "can't combine option -F with option -A or -B" ); + if( isFrameshift() ) + ERR( "can't combine option -F > 0 with option -A or -B" ); } } @@ -559,6 +587,7 @@ void LastalArguments::writeCommented( std::ostream& stream ) const{ stream << " u=" << maskLowercase; stream << " s=" << strand; stream << " S=" << isQueryStrandMatrix; + stream << " M=" << isGreedy; stream << " T=" << globality; stream << " m=" << oneHitMultiplicity; stream << " l=" << minHitDepth; diff --git a/src/LastalArguments.hh b/src/LastalArguments.hh index 9418db1..9f03b33 100644 --- a/src/LastalArguments.hh +++ b/src/LastalArguments.hh @@ -50,7 +50,10 @@ struct LastalArguments{ void writeCommented( std::ostream& stream ) const; // are we doing translated alignment (DNA versus protein)? - bool isTranslated() const{ return frameshiftCost > 0; } + bool isTranslated() const{ return frameshiftCost >= 0; } + + // are we doing translated alignment with frameshifts? + bool isFrameshift() const{ return frameshiftCost > 0; } // how many strands are we scanning (1 or 2)? int numOfStrands() const{ return (strand == 2) ? 2 : 1; } @@ -60,6 +63,7 @@ struct LastalArguments{ int outputType; int strand; bool isQueryStrandMatrix; + bool isGreedy; int globality; // type of alignment: local, semi-global, etc. bool isKeepLowercase; int tantanSetting; diff --git a/src/lastal.cc b/src/lastal.cc index 29b08d2..324833b 100644 --- a/src/lastal.cc +++ b/src/lastal.cc @@ -13,7 +13,6 @@ #include "SubsetMinimizerFinder.hh" #include "SubsetSuffixArray.hh" #include "Centroid.hh" -#include "GappedXdropAligner.hh" #include "AlignmentPot.hh" #include "Alignment.hh" #include "SegmentPairPot.hh" @@ -24,6 +23,7 @@ #include "TantanMasker.hh" #include "DiagonalTable.hh" #include "GeneralizedAffineGapCosts.hh" +#include "GreedyXdropAligner.hh" #include "gaplessXdrop.hh" #include "gaplessPssmXdrop.hh" #include "gaplessTwoQualityXdrop.hh" @@ -48,6 +48,7 @@ using namespace cbrc; struct LastAligner { // data that changes between queries Centroid centroid; + GreedyXdropAligner greedyAligner; std::vector<int> qualityPssm; std::vector<AlignmentText> textAlns; }; @@ -100,7 +101,7 @@ void complementMatrix(const ScoreMatrixRow *from, ScoreMatrixRow *to) { // Set up a scoring matrix, based on the user options void makeScoreMatrix( const std::string& matrixName, const std::string& matrixFile ){ - if( !matrixName.empty() ){ + if( !matrixName.empty() && !args.isGreedy ){ scoreMatrix.fromString( matrixFile ); } else{ @@ -131,6 +132,8 @@ void permuteComplement(const double *from, double *to) { } void makeQualityScorers(){ + if( args.isGreedy ) return; + if( args.isTranslated() ) if( isQuality( args.inputFormat ) || isQuality( referenceFormat ) ) return warn( args.programName, @@ -418,6 +421,7 @@ static const uchar *getQueryQual(size_t queryNum) { static const ScoreMatrixRow *getQueryPssm(const LastAligner &aligner, size_t queryNum) { + if (args.isGreedy) return 0; if (args.inputFormat == sequenceFormat::pssm) return query.pssmReader() + query.padBeg(queryNum); const std::vector<int> &qualityPssm = aligner.qualityPssm; @@ -428,6 +432,10 @@ static const ScoreMatrixRow *getQueryPssm(const LastAligner &aligner, namespace Phase{ enum Enum{ gapless, gapped, final }; } +static bool isMaskLowercase(Phase::Enum e) { + return (e < 1 && args.maskLowercase > 0) || args.maskLowercase > 2; +} + struct Dispatcher{ const uchar* a; // the reference sequence const uchar* b; // the query sequence @@ -446,8 +454,8 @@ struct Dispatcher{ i( text.qualityReader() ), j( getQueryQual(queryNum) ), p( getQueryPssm(aligner, queryNum) ), - m( getScoreMatrix( strand, e < args.maskLowercase ) ), - t( getTwoQualityMatrix( strand, e < args.maskLowercase ) ), + m( getScoreMatrix( strand, isMaskLowercase(e) ) ), + t( getTwoQualityMatrix( strand, isMaskLowercase(e) ) ), d( (e == Phase::gapless) ? args.maxDropGapless : (e == Phase::gapped ) ? args.maxDropGapped : args.maxDropFinal ), z( t ? 2 : p ? 1 : 0 ){} @@ -605,9 +613,8 @@ void alignGapped( LastAligner& aligner, AlignmentPot& gappedAlns, SegmentPairPot& gaplessAlns, size_t queryNum, char strand, const uchar* querySeq, Phase::Enum phase ){ - Centroid& centroid = aligner.centroid; Dispatcher dis( phase, aligner, queryNum, strand, querySeq ); - indexT frameSize = args.isTranslated() ? (query.padLen(queryNum) / 3) : 0; + indexT frameSize = args.isFrameshift() ? (query.padLen(queryNum) / 3) : 0; countT gappedExtensionCount = 0, gappedAlignmentCount = 0; // Redo the gapless extensions, using gapped score parameters. @@ -647,10 +654,10 @@ void alignGapped( LastAligner& aligner, shrinkToLongestIdenticalRun( aln.seed, dis ); // do gapped extension from each end of the seed: - aln.makeXdrop( centroid.aligner(), centroid, dis.a, dis.b, args.globality, - dis.m, scoreMatrix.maxScore, gapCosts, dis.d, - args.frameshiftCost, frameSize, dis.p, - dis.t, dis.i, dis.j, alph, extras ); + aln.makeXdrop( aligner.centroid, aligner.greedyAligner, args.isGreedy, + dis.a, dis.b, args.globality, dis.m, scoreMatrix.maxScore, + gapCosts, dis.d, args.frameshiftCost, frameSize, + dis.p, dis.t, dis.i, dis.j, alph, extras ); ++gappedExtensionCount; if( aln.score < args.minScoreGapped ) continue; @@ -682,7 +689,7 @@ void alignFinish( LastAligner& aligner, const AlignmentPot& gappedAlns, size_t queryNum, char strand, const uchar* querySeq ){ Centroid& centroid = aligner.centroid; Dispatcher dis( Phase::final, aligner, queryNum, strand, querySeq ); - indexT frameSize = args.isTranslated() ? (query.padLen(queryNum) / 3) : 0; + indexT frameSize = args.isFrameshift() ? (query.padLen(queryNum) / 3) : 0; if( args.outputType > 3 ){ if( dis.p ){ @@ -704,11 +711,11 @@ void alignFinish( LastAligner& aligner, const AlignmentPot& gappedAlns, Alignment probAln; AlignmentExtras extras; probAln.seed = aln.seed; - probAln.makeXdrop( centroid.aligner(), centroid, + probAln.makeXdrop( centroid, aligner.greedyAligner, args.isGreedy, dis.a, dis.b, args.globality, dis.m, scoreMatrix.maxScore, gapCosts, dis.d, - args.frameshiftCost, frameSize, dis.p, dis.t, - dis.i, dis.j, alph, extras, + args.frameshiftCost, frameSize, + dis.p, dis.t, dis.i, dis.j, alph, extras, args.gamma, args.outputType ); assert( aln.score != -INF ); writeAlignment( aligner, probAln, queryNum, strand, querySeq, extras ); @@ -716,6 +723,22 @@ void alignFinish( LastAligner& aligner, const AlignmentPot& gappedAlns, } } +static void eraseWeakAlignments(LastAligner &aligner, AlignmentPot &gappedAlns, + size_t queryNum, char strand, + const uchar *querySeq) { + indexT frameSize = args.isFrameshift() ? (query.padLen(queryNum) / 3) : 0; + Dispatcher dis(Phase::gapless, aligner, queryNum, strand, querySeq); + for (size_t i = 0; i < gappedAlns.size(); ++i) { + Alignment &a = gappedAlns.items[i]; + if (!a.hasGoodSegment(dis.a, dis.b, args.minScoreGapped, dis.m, gapCosts, + args.frameshiftCost, frameSize, + dis.p, dis.t, dis.i, dis.j)) { + AlignmentPot::mark(a); + } + } + erase_if(gappedAlns.items, AlignmentPot::isMarked); +} + static bool lessForCulling(const AlignmentText &x, const AlignmentText &y) { if (x.strandNum != y.strandNum) return x.strandNum < y.strandNum; if (x.queryBeg != y.queryBeg ) return x.queryBeg < y.queryBeg; @@ -768,7 +791,7 @@ void makeQualityPssm( LastAligner& aligner, size_t queryNum, char strand, const uchar* querySeq, bool isMask ){ if( !isQuality( args.inputFormat ) || isQuality( referenceFormat ) ) return; - if( args.isTranslated() ) return; + if( args.isTranslated() || args.isGreedy ) return; std::vector<int> &qualityPssm = aligner.qualityPssm; size_t queryLen = query.padLen(queryNum); @@ -804,23 +827,28 @@ void scan( LastAligner& aligner, alignGapless( aligner, gaplessAlns, queryNum, strand, querySeq ); if( args.outputType == 1 ) return; // we just want gapless alignments - if( args.maskLowercase == 1 ) + if( args.maskLowercase == 1 || args.maskLowercase == 2 ) makeQualityPssm( aligner, queryNum, strand, querySeq, false ); AlignmentPot gappedAlns; - if( args.maskLowercase == 2 || args.maxDropFinal != args.maxDropGapped ){ + if( args.maxDropFinal != args.maxDropGapped ){ alignGapped( aligner, gappedAlns, gaplessAlns, queryNum, strand, querySeq, Phase::gapped ); erase_if( gaplessAlns.items, SegmentPairPot::isNotMarkedAsGood ); } - if( args.maskLowercase == 2 ) - makeQualityPssm( aligner, queryNum, strand, querySeq, false ); - alignGapped( aligner, gappedAlns, gaplessAlns, queryNum, strand, querySeq, Phase::final ); + if (args.maskLowercase == 2) { + makeQualityPssm(aligner, queryNum, strand, querySeq, true); + eraseWeakAlignments(aligner, gappedAlns, queryNum, strand, querySeq); + LOG2("lowercase-filtered alignments=" << gappedAlns.size()); + if (args.outputType > 3) + makeQualityPssm(aligner, queryNum, strand, querySeq, false); + } + if( args.outputType > 2 ){ // we want non-redundant alignments gappedAlns.eraseSuboptimal(); LOG2( "nonredundant gapped alignments=" << gappedAlns.size() ); diff --git a/src/makefile b/src/makefile index ea6e797..efbfb2c 100644 --- a/src/makefile +++ b/src/makefile @@ -19,13 +19,14 @@ ScoreMatrix.o SubsetMinimizerFinder.o tantan.o DiagonalTable.o \ SegmentPair.o Alignment.o GappedXdropAligner.o SegmentPairPot.o \ AlignmentPot.o GeneralizedAffineGapCosts.o Centroid.o \ LambdaCalculator.o TwoQualityScoreMatrix.o OneQualityScoreMatrix.o \ -QualityPssmMaker.o GeneticCode.o LastEvaluer.o gaplessXdrop.o \ -gaplessPssmXdrop.o gaplessTwoQualityXdrop.o SubsetSuffixArraySearch.o \ -AlignmentWrite.o MultiSequenceQual.o GappedXdropAlignerPssm.o \ -GappedXdropAligner2qual.o GappedXdropAligner3frame.o lastal.o \ -alp/sls_alignment_evaluer.o alp/sls_pvalues.o alp/sls_alp_sim.o \ -alp/sls_alp_regression.o alp/sls_alp_data.o alp/sls_alp.o \ -alp/sls_basic.o alp/njn_localmaxstatmatrix.o alp/njn_localmaxstat.o \ +QualityPssmMaker.o GeneticCode.o LastEvaluer.o GreedyXdropAligner.o \ +gaplessXdrop.o gaplessPssmXdrop.o gaplessTwoQualityXdrop.o \ +SubsetSuffixArraySearch.o AlignmentWrite.o MultiSequenceQual.o \ +GappedXdropAlignerPssm.o GappedXdropAligner2qual.o \ +GappedXdropAligner3frame.o lastal.o alp/sls_alignment_evaluer.o \ +alp/sls_pvalues.o alp/sls_alp_sim.o alp/sls_alp_regression.o \ +alp/sls_alp_data.o alp/sls_alp.o alp/sls_basic.o \ +alp/njn_localmaxstatmatrix.o alp/njn_localmaxstat.o \ alp/njn_localmaxstatutil.o alp/njn_dynprogprob.o \ alp/njn_dynprogprobproto.o alp/njn_dynprogproblim.o alp/njn_ioutil.o \ alp/njn_random.o alp/sls_falp_alignment_evaluer.o \ @@ -104,7 +105,7 @@ depend: Alignment.o: Alignment.cc Alignment.hh ScoreMatrixRow.hh SegmentPair.hh \ Alphabet.hh Centroid.hh GappedXdropAligner.hh \ GeneralizedAffineGapCosts.hh OneQualityScoreMatrix.hh GeneticCode.hh \ - TwoQualityScoreMatrix.hh + GreedyXdropAligner.hh TwoQualityScoreMatrix.hh AlignmentPot.o: AlignmentPot.cc AlignmentPot.hh Alignment.hh \ ScoreMatrixRow.hh SegmentPair.hh AlignmentWrite.o: AlignmentWrite.cc Alignment.hh ScoreMatrixRow.hh \ @@ -134,6 +135,8 @@ GappedXdropAlignerPssm.o: GappedXdropAlignerPssm.cc GappedXdropAligner.hh \ GeneralizedAffineGapCosts.o: GeneralizedAffineGapCosts.cc \ GeneralizedAffineGapCosts.hh GeneticCode.o: GeneticCode.cc GeneticCode.hh Alphabet.hh +GreedyXdropAligner.o: GreedyXdropAligner.cc GreedyXdropAligner.hh \ + ScoreMatrixRow.hh LambdaCalculator.o: LambdaCalculator.cc LambdaCalculator.hh LastEvaluer.o: LastEvaluer.cc LastEvaluer.hh ScoreMatrixRow.hh \ alp/sls_alignment_evaluer.hpp alp/sls_pvalues.hpp alp/sls_basic.hpp \ @@ -191,9 +194,9 @@ lastal.o: lastal.cc LastalArguments.hh SequenceFormat.hh \ fileMap.hh Centroid.hh GappedXdropAligner.hh \ GeneralizedAffineGapCosts.hh SegmentPair.hh AlignmentPot.hh Alignment.hh \ SegmentPairPot.hh ScoreMatrix.hh Alphabet.hh MultiSequence.hh \ - TantanMasker.hh tantan.hh DiagonalTable.hh gaplessXdrop.hh \ - gaplessPssmXdrop.hh gaplessTwoQualityXdrop.hh io.hh threadUtil.hh \ - version.hh + TantanMasker.hh tantan.hh DiagonalTable.hh GreedyXdropAligner.hh \ + gaplessXdrop.hh gaplessPssmXdrop.hh gaplessTwoQualityXdrop.hh io.hh \ + threadUtil.hh version.hh lastdb.o: lastdb.cc LastdbArguments.hh SequenceFormat.hh \ SubsetSuffixArray.hh CyclicSubsetSeed.hh VectorOrMmap.hh Mmap.hh \ fileMap.hh stringify.hh Alphabet.hh MultiSequence.hh ScoreMatrixRow.hh \ diff --git a/src/version.hh b/src/version.hh index 5d61001..aefb9eb 100644 --- a/src/version.hh +++ b/src/version.hh @@ -1 +1 @@ -"731" +"737" -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/last-align.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
