This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Reading CRAM files in pysam

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
python cram samtools pysam htslib
• 7,173 views • link
updated by Ram Ram • written by QVINTVS_FABIVS_MAXIMVS QVINTVS_FABIVS_MAXIMVS

3 answers

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
• 0 views • link
updated by Ram Ram • written by matted matted

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

• 0 views • link
written by QVINTVS_FABIVS_MAXIMVS QVINTVS_FABIVS_MAXIMVS

What version of pysam are you using?

• 0 views • link
written by matted matted

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!

• 0 views • link
updated by Ram Ram • written by duxan • duxan 
ValueError: fetch called on bamfile without index

Does the input cram file actually have an index?

• 0 views • link
updated by Ram Ram • written by dariober dariober

Log in to answer this question.

More posts like this