Skip to content

Commit

Permalink
avoid hanging in super low entropy sequence
Browse files Browse the repository at this point in the history
We should never be running the entropy calculation past the read end. It won't
ever support a haplotype call to do so.

This change may very slightly change behavior of the algorithm in low entropy
sequence regions.
  • Loading branch information
ekg committed Jun 3, 2019
1 parent 4b8845a commit ffcee0d
Showing 1 changed file with 5 additions and 8 deletions.
13 changes: 5 additions & 8 deletions src/AlleleParser.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1266,19 +1266,16 @@ Allele AlleleParser::makeAllele(RegisteredAlignment& ra,
// a dangerous game
int start = pos - currentSequenceStart;
double minEntropy = parameters.minRepeatEntropy;
// check first that' wer'e actually ina repeat... TODO
//cerr << "entropy of " << entropy(currentSequence.substr(start, repeatRightBoundary - pos)) << " is too low, " << endl;
while (minEntropy > 0 && // ignore if turned off
repeatRightBoundary - currentSequenceStart < currentSequence.size() && //guard
// don't run off the end of the current sequence
repeatRightBoundary - currentSequenceStart < currentSequence.size() &&
// there is no point in going past the alignment end
// because we won't make a haplotype call unless we have a covering observation from a read
repeatRightBoundary < alignment.ENDPOSITION &&
entropy(currentSequence.substr(start, repeatRightBoundary - pos)) < minEntropy) {
//cerr << "entropy of " << entropy(currentSequence.substr(start, repeatRightBoundary - pos)) << " is too low, ";
//cerr << "increasing rought boundary to ";
++repeatRightBoundary;
//cerr << repeatRightBoundary << endl;
}

// now we
//cachedRepeatCounts[pos] = repeatCounts(pos - currentSequenceStart, currentSequence, 12);
// edge case, the indel is an insertion and matches the reference to the right
// this means there is a repeat structure in the read, but not the ref
if (currentSequence.substr(pos - currentSequenceStart, length) == readSequence) {
Expand Down

0 comments on commit ffcee0d

Please sign in to comment.