Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Segmentation Fault Under Certain Conditions #535

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

evolvedmicrobe
Copy link

@evolvedmicrobe evolvedmicrobe commented Dec 12, 2018

This PR was motivated by a SIGSEGV that STAR produced when aligning a particular read against the Rabbit genome. The problem is reproducible, and a file with the rabbit genome and one read that produces the error is available here.

After downloading, the following command will exit with a SIGSEGV

STAR --genomeDir rabbit --outSAMmultNmax -1 --readFilesIn temp.fastq 

I'd love to integrate a fix for this problem, and attempted one with this PR.

I diagnosed the issue 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 consequently 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:

ACTGAAAA <- iSA1 / first prefix
ACTGNNNN <- iSA2
ACTTACTG <- next prefix

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 prefix length tested, I also changed the binary search to account for this. Testing showed for particular read this had these changes restored the correct behavior in the inferred bounds (as determined by doing a full binary search).

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:

    ACTGAAAA <- iSA1 / first prefix
    ACTGNNNN <- iSA2
    ACTTACTG <- next prefix

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.
@alexdobin
Copy link
Owner

Hi Nigel,

thanks a lot for a very thorough investigation of this problem.
I am generally not happy with this part of the code as it also produces seg-faults for really small genomes. As you pointed out, the real problem is in the way this index is generated in the indexing step, the locations very close to the genome ends should be excluded. I am reluctant to remove the part of the code that you suggested as it may change its behavior unpredictably in other cases.
I will work on your test case to find a consistent solution.

Cheers
Alex

@evolvedmicrobe
Copy link
Author

Hi Alex,

Thanks for taking a look at being willing to solve the problem, I thought you might have a better fix, will look forward to it!

Thanks again,
Nigel

@evolvedmicrobe
Copy link
Author

Hi Alex,

I hope you're doing well in the New Year! Just wanted to inquire quickly if you thought you might have come up with a better and more consistent solution.

All the best,
Nigel

@alexdobin
Copy link
Owner

Hi Nigel,

sorry for the delay, I was busy with the new release and could not look into this problem. I hope to get back to it within the next week or two.

Cheers
Alex

@evolvedmicrobe
Copy link
Author

Wonderful, thank you Alex!

@evolvedmicrobe
Copy link
Author

Hi @alexdobin,

Just wanted to ping on this issue in case you might have a chance to revisit it. We're continuing to see a lot of segfaults when people align things to the rabbit genome that we'd like to avoid. I would be happy to help out with anything needed if I can as well.

Cheers,
Nigel

@alexdobin
Copy link
Owner

Hi Nigel,

have you tried it with one of the latest versions?
The seg-fault for small (or masked) genomes issue was fixed in 2.7.4a. I would hope that this issue might have been fixed as well?

Cheers
Alex

@alexdobin alexdobin added this to the 2.7.10 milestone Sep 28, 2021
alexdobin added a commit that referenced this pull request Sep 30, 2021
…seg-fault whem mapping to the rabbit genome. Issue #1223: fixed the N_unmapped value reported in ReadsPerGene.out.tab. The single-end (i.e. partially mapped alignment are not excluded from N_unmapped. dev_EoI_2.7.9a_2021-09-30
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants