Fix several bugs involving getReadCoordinateForReferenceCoordinate #6485
Conversation
davidbenjamin
force-pushed
the
db_read_position
branch
from
ed56055 to
06a2021
Compare
davidbenjamin
force-pushed
the
db_read_position
branch
from
06a2021 to
aa2f2ba
Compare
|
Since you are working on this branch right as I have noticed a bug in this code myself. I have one but to pile onto this one. I have discovered that if you run getReadCoordinateForReferenceCoordinate() with allowGoalNotReached==true and if you request exactly the base one position after the end of the read window you will have it return that base even if the read doesn't actually have that many elements in its base array (for example read spanning from 1000-1019 with cigar 20M requesting position 1020 will return a pair of <21, false> instead of <-1,false>. It doesn't look like this is one of the bugs that you described but if you are fixing this code It would be nice to fix it as well. |
|
@jamesemery This PR fixes that. |
src/main/java/org/broadinstitute/hellbender/utils/read/ReadUtils.java
Outdated
Show resolved
Hide resolved
src/main/java/org/broadinstitute/hellbender/utils/read/ReadUtils.java
Outdated
Show resolved
Hide resolved
src/main/java/org/broadinstitute/hellbender/utils/read/ReadUtils.java
Outdated
Show resolved
Hide resolved
| final byte[] myBases = getBases(); | ||
|
|
||
| // can't insert if we don't have any sequence after the inserted alt allele to span the new variant | ||
| if (haplotypeInsertLocation + refAllele.length() >= myBases.length) { |
There was a problem hiding this comment.
Is this checking that there are non-variant bases (i.e. bases that are not in the variant context) on either side of the inserted allele?
There was a problem hiding this comment.
If the ref allele has n bases, you can think of the alt allele as deleting n - 1 bases (eg ACT -> A) followed, optionally by replacing those n - 1 bases (in the case of a MNP eg ACT -> A -> AGG). Either way, if the haplotype doesn't contain those n - 1 bases it doesn't make sense to perform this transformation.
That is, suppose our haplotype is . . .ATTG. Then we can insert an ATT -> A deletion by removing the TT, but we can't insert an ATTGC -> A deletion because the C is not there.
Closes #6342.
Closes #6314.
Closes #6294.
Closes #5492.
@takutosato Fixing bugs and simplifying code.