Possible reference_ends bug?

[owensnick]

Came across this comparing pyft (v0.2.6) on a CRAM file annotated with fibertools-rs v0.6.2 to our own independent processing of same bam file, in which we found discrepant nucleosomes or MSP positions on the genome with our interpretation.

If a nucleosome or MSP end on read coincides with the end of an align block then the length of that feature in reference coordinates may be shifted to the start of the next align block making the MSP/nucleosome artificially large. Illustration:

             0         1         2
             01234567890123456789012
Read MSP:        | MSP  |
Read         ===========-------=====
Genome       =======================
Genome MSP:      |     MSP     |

| marks start and end in 0-based exclusive coordinates
= aligned base
- deletion in read 

In the above the read MSP is [4, 11) length 7 but the genome MSP becomes [4, 18) length 12. The issue is that although position 11 is not in the MSP, it maps in the subsequent alignblock. If you calculate the MSP in 1-based inclusive coords the MSP is [5, 11] on read and [5, 11] on genome both length 7.

Attached is a sam record read in which this occurs (mapped to hg19). The nucleosome at forward pos 3525, reverse pos 169, and length 103, is mapped to [3367501, 3367606) by fibertools, but I believe it should be [3367501, 3367604).

Title says possible bug as I can imagine there is some subjectivity on how these align blocks are handled. But in the above motivating example, I think the length of the deletion following a msp/nucleosome should not influence the length of that feature on the genome.

nuc_example.sam.txt

[mrvollger]

I agree this is something to work on. Internally I thought I was doing all these lift overs with inclusive coordinates to avoid this (even though I output BED style coordinates). I'll take a look and see what's going on here.

If you have your ranges against hg19 for this read as well I could use those to help develop some test cases.

[owensnick]

Attached is a table of your coordinates and my coordinates for all intervals for that read. It is only the single nucleosome where this issue occurs. I can easily generate similar tables for other discrepant reads if that is useful.

If it is of any help we have a super simple base level align map, just a vector that maps read to genome coords. The logic for mapping intervals is very straightforward with that data structure - only a few lines of code. This may be useful for testing. This takes more memory than align blocks, but preallocating and reusing over multiple reads we have found the performance is pretty good.

nuc_example_coords.tsv

[mrvollger]

I did get the chance to confirm at least one of these corner cases and I plan to release at least that fix in v0.8.

I'd be very interested in using your data structure to validate my results in test cases. Does it happen to be a rust library I could load?

[owensnick]

That is great.

Unfortunately this is not in Rust - we have this in a Julia library that uses a wrapper around hts lib. I think the Rust calling Julia is doable, but it may not be worth it in this case. We came across this as we put comparison to fibertools in our test case - although we are not calling fibertools directly from Julia, rather we got fibertools output for a test bam and compare to that. I can show you how to call our Julia library or alternatively, I can send you all the divergences for our test bam. There were more than this.

[mrvollger]

Would you be willing to check if one of my test files (tests/data/ctcf.bam) has any of these issues? I'd like to avoid brining new test datasets into the repo. Unless of course its needed.

Thanks!