From c41857813111e6eda71db8e224543aa49b9f4ad9 Mon Sep 17 00:00:00 2001 From: Nigel Delaney Date: Tue, 4 Dec 2018 23:36:40 -0800 Subject: [PATCH 1/2] Fix Segmentation Fault Under Certain Conditions This PR was motivated by a SIGSEGV that STAR produced when aligning a particular read against the Rabbit genome. Explanation of the problem and fix below. The issue was as follows. In trying to map the read, the aligner routine within `ReadAlign::maxMappableLength2strands` was trying to find exact matches for the end of the read. The first step in this process is looking up a given prefix length in the suffix array index. This index contains minimum (and maximum) locations for prefixes up to size L_max (14 in this case). For this particular read segment, no prefix of size 14 was found, so the program re-attempted to search the index with a prefix of size 13, which it did find. Because a prefix of size 13 and not 14 was found, the program knew there was no match of size 14 and so tried to set the end of the search range in the suffix array (`indStartEnd[1]`) using a short-cut that avoids the binary search. In particular, it set the End of the bounded range to the value of one minus the next element in the index. However, it appears not all elements prior to the start of the next prefix are guaranteed to match the current prefix. In this particular case, that next element pointed to a position 10 bases before the end of the genome, which meant it could not produce an alignment 13 bases in length. Later on inside `ReadAlign::stitchPieces` and its call to `sjSplitAlign` this led to the SIGSEGV when while calculating the alignment start position it inferred to be 3 basepairs past the end of the genome, giving an offset position of `-3` which when converted to an unsigned integer has value `18446744073709551613`, and consequentally the calculated value for `isj` in the code line `P->sjDstart[isj]` referred to unmapped memory, leading to the SIGSEGV. I believe this can happen as the suffix array can be sorted as follows: ACGTACGTACGTAA # Matching prefix of size 13 ACGTACGTACGTAN # Suffix array element with padded-bases/Ns, but not a valid 13-mer alignment ACGTACGTACGTCA # Next prefix of size 13 Since it seemed the short-cut avoiding the binary-search was not reliable, I removed it with this commit. Since it also appeared that this was a general problem where one could not guarantee the initial range guessed held matching prefixes up to the length tested (but rather only the length tested minus 1), I also changed the binary search to account for this. Testing showed for this case it restored the correct behavior in the inferred bounds, such that the alignment length was consistent across the start and end of the array element range set. --- source/ReadAlign_maxMappableLength2strands.cpp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/source/ReadAlign_maxMappableLength2strands.cpp b/source/ReadAlign_maxMappableLength2strands.cpp index 9db6618c..f63b87da 100644 --- a/source/ReadAlign_maxMappableLength2strands.cpp +++ b/source/ReadAlign_maxMappableLength2strands.cpp @@ -58,12 +58,7 @@ uint ReadAlign::maxMappableLength2strands(uint pieceStartIn, uint pieceLengthIn, maxL=0; Nrep = maxMappableLength(Read1, pieceStart, pieceLength, G, SA, iSA1 & P->SAiMarkNmask, iSA2, dirR, maxL, indStartEnd, P); #else - if (Lind < P->genomeSAindexNbases && (iSA1 & P->SAiMarkNmaskC)==0 ) {//no need for SA search - indStartEnd[0]=iSA1; - indStartEnd[1]=iSA2; - Nrep=indStartEnd[1]-indStartEnd[0]+1; - maxL=Lind; - } else if (iSA1==iSA2) {//unique align already, just find maxL + if (iSA1==iSA2) {//unique align already, just find maxL if ((iSA1 & P->SAiMarkNmaskC)!=0) { ostringstream errOut; errOut << "BUG: in ReadAlign::maxMappableLength2strands"; @@ -75,7 +70,7 @@ uint ReadAlign::maxMappableLength2strands(uint pieceStartIn, uint pieceLengthIn, maxL=compareSeqToGenome(Read1,pieceStart, pieceLength, Lind, G, SA, iSA1, dirR, comparRes, P); } else {//SA search, pieceLength>maxL if ( (iSA1 & P->SAiMarkNmaskC)==0 ) {//no N in the prefix - maxL=Lind; + maxL=Lind - 1; // There can be a suffix "N" base between this prefix and the next highest } else { maxL=0; }; From 927e0c7e5b5def9b487750090586b3e4342b9412 Mon Sep 17 00:00:00 2001 From: "nigel.delaney" Date: Wed, 5 Dec 2018 01:04:39 -0800 Subject: [PATCH 2/2] Fix assumption that difference was only 1 For sparse suffix arrays, presumably it could be more. --- source/ReadAlign_maxMappableLength2strands.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/source/ReadAlign_maxMappableLength2strands.cpp b/source/ReadAlign_maxMappableLength2strands.cpp index f63b87da..c0ddee78 100644 --- a/source/ReadAlign_maxMappableLength2strands.cpp +++ b/source/ReadAlign_maxMappableLength2strands.cpp @@ -70,7 +70,14 @@ uint ReadAlign::maxMappableLength2strands(uint pieceStartIn, uint pieceLengthIn, maxL=compareSeqToGenome(Read1,pieceStart, pieceLength, Lind, G, SA, iSA1, dirR, comparRes, P); } else {//SA search, pieceLength>maxL if ( (iSA1 & P->SAiMarkNmaskC)==0 ) {//no N in the prefix - maxL=Lind - 1; // There can be a suffix "N" base between this prefix and the next highest + // iSA2 can now point to a suffix with multiple "N" bases, eg. this sort order + // ACTGAAAA <- iSA1 + // ACTGNNNN <- iSA2 + // ACTTACTG <- next prefix + // As a result, we can only assume the match length is as good up to + // the length that iSA2 hits, and need to compare everything past that length in the binary search + bool cres; + maxL=compareSeqToGenome(Read1,pieceStart, pieceLength, 0, G, SA, iSA2, dirR, cres, P); } else { maxL=0; };