1
$\begingroup$

Ideally, I'd like to do something like the following:

f = pysam.AlignmentFile(filename, index_filename=index_filename)
read = list(f.fetch('Chromosome', 1, 10000))[0]
read.get_variations()

And have that output something like:

[
  (1, 'A', 'C'), # An A -> C mismatch at position 1
  (15, 'I', 'CGGG'), # An insertion of the string 'CGGG' starting at position 15 of the read
  (32, 'D', 4), # A deletion of 4 base pairs starting at position 32
]

I know that it's possible to figure out all of those values by parsing the CIGAR string and MD tags but I'm wondering if there's a more intuitive interface.

$\endgroup$
1
$\begingroup$

You're looking for pileup, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            ...handle mismatches...
$\endgroup$
2
  • $\begingroup$ Thanks for the response! I think this would work but it's less than ideal since I would be iterating over each read multiple times, no? I'm wondering if it's possible to iterate over the list of reads that span a region and retrieve that same information without going over pileup columns. $\endgroup$
    – juniper-
    Dec 16 '20 at 16:12
  • $\begingroup$ Also, this bit of code for pileupread in pileupcolumn.pileups crashes the Python kernel when running inside of an IPython notebook :-/ $\endgroup$
    – juniper-
    Dec 16 '20 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.