Extract mapped genome sequence from mappy API

[marcus1487]

I am using the mappy API with the sole end goal of extracting the genomic sequence for a read mapping. Currently, I am using pyfaidx to take the mapped coordinates from the mappy.Alignment object and extract the associated genomic sequence. This seems like a bit of extra work since the sequence in theory may be extract-able from the mappy.Aligner or mappy.Alignmentobjects (though there are likely reasons this may not be true). This would also allow my code not to load the same reference fasta file twice (once for mapping in mappy and the second time for randomly extracting genome sequence). Is there access to this right now via the python API (maybe undocumented)? Or would it be possible to add this feature either directly or indirectly to the mappy.Alignment object?

[lh3]

Not possible at the moment, but it should not be hard to add this functionality.

[lh3]

Added. You can now:

import mappy as mp

a = mp.Aligner("test/MT-human.fa")
seq = a.seq("MT_human", 100, 250)
for hit in a.map(mp.revcomp(seq)):
	print(hit)

I have only tested it on toy examples. Let me know if it has issues.

[marcus1487]

I have tested this feature and it works perfectly for my test cases.

[marcus1487]

After some more extensive testing it does not look as though this feature is working correctly. When working on some larger genomes the sequence extraction appears to be broken.

As a test case extracting the following genomic location from the linked zebrafish genome returns None, but this is a valid genomic location and I can confirm that it contains only standard bases.

Genomic location: 4:70000000-70000100
genome file source: ftp://ftp.ensembl.org/pub/release-89/fasta/danio_rerio/dna/Danio_rerio.GRCz10.dna.toplevel.fa.gz

Through a binary search I've found that sequence is only able to be extracted for each chromosome up to position 58871916. This is true for every chromosome in this particular fasta file.

[lh3]

It is a bug. It is now fixed via 96b132c.