Reading CRAM files in pysam

Question

[QVINTVS_FABIVS_MAXIMVS]

I'm new to CRAM format and I am wondering if I can read a CRAM file in pysam or samtools.

Or generally speaking how do I use a CRAM file like a BAM file?

#!/usr/env python
import pysam
CRAM_FH = '/home/admin/data/HG00096.alt_bwamem_GRCh38DH.20150718.GBR.low_coverage.cram'
cram = pysam.AlignmentFile(CRAM_FH,'rb')
for read in cram.fetch(region='chr1:10000000-10000100'):
        print read
cram.close()

Returned:

  File "src/test_cram.py", line 5, in <module>
    for read in cram.fetch(region='chr1:10000000-10000100'):
  File "pysam/calignmentfile.pyx", line 874, in pysam.calignmentfile.AlignmentFile.fetch (pysam/calignmentfile.c:10986)
ValueError: fetch called on bamfile without index

Answer 1

[matted]

It should. Pysam uses recent versions of htslib, and if you look at the release notes for pysam, you'll see this entry with version 0.8.2 (the latest is 0.8.4):

• Pysam now wraps htslib 1.2.1 and samtools version 1.2
• Added CRAM file support to pysam

[QVINTVS_FABIVS_MAXIMVS]

If the index is in the same directory why wont the code work then?

[matted]

What version of pysam are you using?

[QVINTVS_FABIVS_MAXIMVS]

most recent

Answer 2

[duxan]

It seems that pysam doesn't recognize your file as CRAM. You can check this by: cram.is_bam and cram.is_cram.

From the source you can see a way how you can hack it (since I haven't found this in documentation, it is safe to call it a hack) by setting mode to 'rc' instead of 'rb'. This same issue is reported here. I can confirm that changing mode this way worked in my setup (CRAM built with samtools view -C BAM > CRAM, not with pysam, and also index built with samtools index CRAM).

Hope it helps!

Answer 3

[dariober]

ValueError: fetch called on bamfile without index

Does the input cram file actually have an index?

[QVINTVS_FABIVS_MAXIMVS]

yes the .cram.crai is in the same directory