Commit 513a047b authored by Mathieu Giraud's avatar Mathieu Giraud

core/segment.{cpp,h}: computes only half of the DP matrix for VJ segments

Starting from the full read, we can not limit the DP computations
to a k-band around the diagonal without first knowing where is exactly the junction.

Nevertheless, we can avoid computing about one half of the DP matrix,
as the end of the V / the start of the J (minus some deletions) must be matched.

The BOTTOM_TRIANGLE_SHIFT is now set to 20, and this should be large enough
to handle V/J deletions until ~30 bp (see comment in segment.h).
(The current tests were even passing with BOTTOM_TRIANGLE_SHIFT set to 10.)

Now the FineSegmenter (as launched by 'make shouldvdj_with_rc_merged') is about 35% faster.
parent 5a23f9d6
......@@ -586,7 +586,7 @@ bool comp_pair (pair<int,int> i,pair<int,int> j)
* @param read: the read
* @param rep: a collection of reference sequences
* @param reverse_both: if true, reverse both the read and the reference sequences (J segment)
* @param local: if true, Local alignment (D segment), otherwise LocalEndWithSomeDeletions (V and J segments)
* @param local: if true, Local alignment (D segment), otherwise LocalEndWithSomeDeletions and onlyBottomTriangle (V and J segments)
* @param segment_cost: the cost used by the dynamic programing
* @return score: the maximized score, together with the following values:
* @return tag: - the name of the sequence
......@@ -620,7 +620,9 @@ int align_against_collection(string &read, Fasta &rep, bool reverse_both, bool l
dpMode, // DynProg::SemiGlobalTrans,
segment_cost, // DNA
reverse_both, reverse_both);
int score = dp.compute();
bool onlyBottomTriangle = !local ;
int score = dp.compute(onlyBottomTriangle, BOTTOM_TRIANGLE_SHIFT);
if (local==true){
dp.backtrack();
......
......@@ -34,6 +34,10 @@
#define THRESHOLD_NB_EXPECTED 1.0 /* Threshold of the accepted expected value for number of found k-mers */
#define THRESHOLD_NB_EXPECTED_D .05 /* e-value threshold, D-REGION */
#define BOTTOM_TRIANGLE_SHIFT 20 /* Should equal to (max allowed 'k-band') + (max allowed number of V/J deletions) - (min size to recognize facing J/V)
As we need ~10 bp to recognize the facing V/J, this value should be large enough to handle V/J deletions until ~30 bp,
(and even larger V/J deletions if there is a large facing J or V in the read). */
using namespace std;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment