Zero length read alignments make ReadWalkers crash.

[vruano]

Doing some empirical data testing I found some instance of reads that have no single based aligned with the reference.

E.g. CIGAR: 50I2S so insertion followed by soft-clip.

That causes ReadWalker to crash when trying to create SimpleInterval on the read with a IAE.

I guess the solution is to add additional Wellformed filter.

[lbergelson]

Adding a new filter seems like the right thing. Why is this read aligned like this? It seems like a bug...

[vruano]

No clue, if not BWA (or whatever aligner) some post-alignment editing that went wrong.

I have seem occurrences of this in at least 3 WEx bams that I presume have been produced using good practices.

[vruano]

Test fail in ApplyBQSRIntegrationTest after the change so It might be some of these in the test input.

[akiezun]

what does GATK3 do?

[droazen]

In GATK3, this check is performed by the BadCigarFilter

[droazen]

Perhaps we should just port BadCigarFilter in its entirety, or would that be overkill? It is a pretty expensive read filter...

[akiezun]

yes, port. worry about speed later

On Tue, Apr 14, 2015 at 3:13 PM, droazen [email protected] wrote:

Perhaps we should just port BadCigarFilter in its entirety, or would that
be overkill? It is a pretty expensive read filter...

—
Reply to this email directly or view it on GitHub
#373 (comment)
.

[akiezun]

@vruano can you clarify? The specific CIGAR "50I2S" that you mentioned is not flagged as bad by BadCigarFilter in GATK3. Can you provide test cases for this issue?

[vruano]

I see, then I guess is another filter that is getting rid of it in GATK3... @droazen pointed out to that filter class and I guess we just assumed that BadCigarFilter was getting rid of it without actually checking whether that is true.

[droazen]

It appeared to me upon a quick inspection that the old BadCigarFilter was doing this check (via its "hasMeaningfulElements" condition), but I guess I was mistaken. It seems to be making sure that there is at least one operator that consumes read bases, but doesn't do the equivalent check for operators that consume reference bases.

[vruano]

Then should be another filter or it might be the case that GATK3 does not get rid of it but it tolerates them ; it does not blow like Hellbender/ReadWalker is doing at the moment.

[akiezun]

@vruano so a bam file with 1 read and cigar 50I2S would work in GATK3 (for print reads) and blow up on hellbender?

[droazen]

If we want the hellbender engine to tolerate these reads, we can probably do that very easily -- just need to modify ReadWalker to set readInterval to null for these reads, like it currently does with unmapped reads:

            .forEach(read -> {
                final SimpleInterval readInterval = read.isUnmapped() ? null :
                                                                        new SimpleInterval(read.getContig(), read.getStart(), read.getEnd());
                apply(read,
                      new ReferenceContext(reference, readInterval), // Will create an empty ReferenceContext if reference or readInterval == null
                      new FeatureContext(features, readInterval));   // Will create an empty FeatureContext if features or readInterval == null
            });

We should confirm that GATK3 actually tolerates these reads, though -- if it instead filters them out, I'd advocate for filtering them in hellbender as well.

[vruano]

The way I got an exception is by running a ReadWalker on one of the CEU WEx...

/humgen/gsa-hpprojects/dev/valentin/WExCNv/bams/CEUTrio.HiSeq.WEx.b37_decoy.NA12891.bam

using the broad target file
/humgen/gsa-hpprojects/GATK/bundle/current/b37/Broad.human.exome.b37.interval_list

I presume that that bams have been extensively analyzed using GATK3.X and older and so my conclusion/conjecture.

Since I'd got the same kind of problem on another file

/humgen/gsa-hpprojects/dev/valentin/WExCNv/bams/CEUTrio.HiSeq.WEx.b37_decoy.HG03006.bam

now I'm not sure which one is the one that includes that particular CIGAR but I presume is the same error mode.

[droazen]

What is the appropriate check for these reads? if ( read.getCigar().getReferenceLength() == 0 ) ?

[vruano]

I guess that an N (or P?) operation in the CIGAR may result in a non-zero reference length without guarantee that there is a single base of the read aligned anywere.

[vruano]

Perhaps an efficient form of

read.getAlignmentBlocks().stream().mapToInt(AlignmentBlock::getLength).sum() > 0

[droazen]

That is quite a check to be performing for every single read just to handle this edge case....perhaps a filter would be better?

[vruano]

That would be the filter :)
But yeah is inefficient.

[droazen]

(also, getAlignmentBlocks() has no equivalent in the GA4GH Read API, so this approach wouldn't be viable)

[vruano]

What do they have instead?

[vruano]

Can you share a pointer to that API?

[droazen]

It's just not in the API. I suppose if we needed it we could write a utility method that constructs the alignment blocks from the Cigar.

[droazen]

Sure: http://ga4gh.org/documentation/api/v0.5.1/ga4gh_api.html#/schema/org.ga4gh.GAReadAlignment

[akiezun]

pull req #380 is now merged. assigning to @vruano to decide if this is enough for this issue and what to do if it's not

[vruano]

Afaik, is not as the original CIGAR reported above "50I2S" remains a valid one and the the exception would come up.

There are at least three ways to handle this:

1. Make 0-length SimpleIntervals valid.
2. All ReadWalkers filter reads with no base aligned to the reference.
3. Consider no base aligned reads bad formed so they are generally excluded.

[droazen]

There is a 4th option, which I mentioned above -- modify ReadWalker to set readInterval to null for these reads, like it currently does with unmapped reads:

        .forEach(read -> {
            final SimpleInterval readInterval = read.isUnmapped() ? null :
                                                                    new SimpleInterval(read.getContig(), read.getStart(), read.getEnd());
            apply(read,
                  new ReferenceContext(reference, readInterval), // Will create an empty ReferenceContext if reference or readInterval == null
                  new FeatureContext(features, readInterval));   // Will create an empty FeatureContext if features or readInterval == null
        });

This approach has the advantage of allowing ReadWalkers to process these reads if they want to -- they'll just get empty (but non-null) ReferenceContexts and FeatureContexts, just like they do for unmapped reads.

[vruano]

Could you comment on what would be the problem in allowing for SimpleInterval to have 0-length.

[droazen]

0-length intervals open up a huge can of worms (not much code is designed to be able to handle them, and none of the htsjdk query interfaces we rely on work properly with them). Modifying ReadWalker to set readInterval to null for these reads is a more conservative fix for the problem at hand that allows these reads to be processed without trying to query overlapping reference/feature data for them.

[vruano]

I guess that those problems are all GenomeLoc(Parser) inheritance and since we are breaking away from GenomeLoc perhaps is a good time to lift this restriction as well; I don't think being conservative at this point is necessarily a good thing considering the longer road ahead.

[droazen]

The right way to approach this ticket is not to make a huge change that requires a massive amount of work across both our codebase and our dependencies, but instead to make a targeted fix that addresses the original bug. If we want to support zero-length intervals, that should be a separate ticket, as it's a non-trivial task to do correctly.

I'll add: unless the semantics of using zero-length intervals are precisely defined in every case, we should probably not allow them. Eg., what should happen when you query using a zero-length interval? Should you always get back nothing, or should you get back abutting records? How does overlap work? Allowing intervals to be zero-length complicates everything immensely -- what is the tangible benefit of allowing them?

[vruano]

In that case, whatever the people wanna do. 4 then?

[akiezun]

closed by #1474