Hfrxdjt

I am used to gzip/biopython solutions when dealing with sequencing data, but now I wish to switch to more elegant pysam. So I looked at the manual, but run into quite bizarre troubles with the first couple of lines using my bam file

import pysam
samfile = pysam.AlignmentFile("3_Tms_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tms_b3v08_scaf000159'):
 print(read)
samfile.close()

returns ValueError: fetching by region is not available for SAM files. Well, the file is bam. I tried to google the error but the only hits I found are the lines in the source code of pysam that check if the file is bam/cram or sam, so somehow pysam thinks that my bam is a sam. How can I convince him otherwise? I also have noticed that the manual is for python 2.7, that's maybe where the problem comes from...

If you wish to reproduce my problem, here is my bam file.

You don't have the index file for the bam file you're using. I used this script on the file you linked to (changing the name so that it corresponds to the right file and a contig in that file):

#!/usr/bin/env python3

import pysam
samfile = pysam.AlignmentFile("3_Tce_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
 print(read)
samfile.close()

The directory I ran it in had:

$ ls 3*
3_Tce_1_mapped.bam

And I got the error you described:

$ foo.py
Traceback (most recent call last):
 File "/home/terdon/scripts/foo.py", line 5, in <module>
 for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
 File "pysam/libcalignmentfile.pyx", line 1107, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetching by region is not available for SAM files

However, indexing the bam file fixed it:

$ samtools sort 3_Tce_1_mapped.bam > 3_Tce_1_mapped.sorted.bam 
$ mv 3_Tce_1_mapped.sorted.bam 3_Tce_1_mapped.bam
$ samtools index 3_Tce_1_mapped.bam 
$ ls 3*
3_Tce_1_mapped.bam 3_Tce_1_mapped.bam.bai

$ foo.py | wc
227 2724 16725

So just sort and index your files before attempting to seek in them. Which makes sense since the index's job is primarily to allow seeking.

Your 3_Tms_1_mapped.bam file, despite its filename extension, is in fact a bgzippped SAM file. You can verify this using htsfile, which is a small utility packaged with HTSlib:

$ htsfile 3_Tms_1_mapped.bam 
3_Tms_1_mapped.bam: SAM version 1.3 BGZF-compressed sequence data

(For files that really are in BAM format, it reports BAM version 1 compressed sequence data.)

So the error message is accurate in this case.