Pysam not finding mapped reads in BAM file
1
1
Entering edit mode
2.2 years ago
taraeicher ▴ 50

I have a BAM file with the following stats from samtools flagstat:

8022676 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
1336910 + 0 duplicates
8022676 + 0 mapped (100.00% : N/A)
8022676 + 0 paired in sequencing
4011338 + 0 read1
4011338 + 0 read2
8022676 + 0 properly paired (100.00% : N/A)
8022676 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

However, when I run the following code interactively in python:

import pysam
samfile = pysam.AlignmentFile("mybam.bam", "rb")
samfile.mapped

it's printing out 0. Interestingly, it does the same if I try to print the unmapped reads. Is there something that pysam looks for in a BAM file that is different from what samtools flagstat does? If so, is there any way I can convert my BAM file into a format that is compatible with pysam?

BAM ChIP-Seq PySam Samtools • 867 views
ADD COMMENT
1
Entering edit mode
2.2 years ago
taraeicher ▴ 50

After some investigating, I found that it was the index to the BAM file that was causing the problem. I generated this index using bamtools index. Given that this happened with several of my files which were indexed in this way, it looks like there may be a compatibility issue between my version of bamtools (2.2.2) and pysam.

I'm still not sure exactly where the bug is, but I was able to run the following instead of bamtools index as a workaround:

import pysam
samfile = pysam.index("mybam.bam")
ADD COMMENT

Login before adding your answer.

Traffic: 2842 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6