This is a test version of Biostars. For the public version, visit https://www.biostars.org.
pysam fetch: find desired position within short read

Using pysam fetch, I can find & iterate over all short reads that map to a given position/column (1000, say) in the reference genome:

BamFile = pysam.AlignmentFile(FileName, "rb")
position=1000
ShortReadsMappingHere = BamFile.fetch("MyReference", position-1, position)
for ShortRead in ShortReadsMappingHere:
  # do something

What's the most efficient way to find where the position I'm interested in (here 1000) is within each ShortRead? i.e. ShortRead.seq[0] will give me the first element of this short read, but which element is the position I want?
This line of code must get it right:

  ThisPositionInThisRead = ShortRead.get_reference_positions(full_length=True).index(position-1)

because it finds, for every position in the short read, the position in the reference genome that it maps to (or 'None' if it doesn't map) and retrieves the index of the one mapping to the position I want. It's clearly super inefficient however. The following line

  ThisPositionInThisRead = position - ShortRead.reference_start + ShortRead.query_alignment_start - 1

almost always gets it right, by comparing coordinates, but can get it wrong if this short read has an indel not present in the consensus.
Is there a simple (and correct) way of determining this number based on coordinates?

pysam fetch mapping short reads
• 6,703 views • link
updated by Devon Ryan Devon Ryan • written by c.wymant • c.wymant

1 answer

Have a look at the mpileup functions (I've not used it with pysam, just the base C API). It'll let you iterate over reads covering a given position and indicate for you which base in a read actually covers that position (it also indicates insertions/deletions).

An example of doing this with the C API (the function names are more or less the same, since pysam is a wrapper around the htslib C API) can be found in PileOMeth (within each pileup, the qpos value indicates the actual position within the alignment).

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this