Reading CRAM files in pysam
3
1
Entering edit mode
9.0 years ago

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 • 5.6k views
ADD COMMENT
1
Entering edit mode
9.0 years ago
matted 7.8k

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
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

What version of pysam are you using?

ADD REPLY
0
Entering edit mode

most recent

ADD REPLY
1
Entering edit mode
9.0 years ago
duxan ▴ 70

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!

ADD COMMENT
0
Entering edit mode
9.0 years ago
ValueError: fetch called on bamfile without index

Does the input cram file actually have an index?

ADD COMMENT
0
Entering edit mode

yes the .cram.crai is in the same directory

ADD REPLY

Login before adding your answer.

Traffic: 913 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