CRAM support - why not?

[pd3]

It appears CRAM is not supported by Rsamtools. I wonder why is that? Samtools support CRAM natively, so I would naively expect this should be a trivial task, no?

This issue came up from other programs that rely on Rsamtools, please see here vplagnol/ExomeDepth#17

[vjcitn]

Thanks for your note. There hasn't been strong demand for CRAM functionality. Rhtslib is up to htslib 1.15, and the latter has moved to 1.18. Resources in Bioconductor core are very tight and it has not been clear whether Rsamtools should be upgraded to handle CRAM, or another htslib-dependent package should be introduced. If you see a way to modifying Rsamtools to deal with CRAM in a way that is useful to you we would surely work with a PR.

[hpages]

FWIW some work was done by @mtmorgan last year to support CRAM, see the cram branch, but I don't know how far Martin went or how much work would still be needed before this branch can be merged into devel.

[hpages]

UPDATE: The latest Rhtslib (2.99.2) now includes htslib 1.18 (was 1.15). See https://github.com/Bioconductor/Rhtslib/blob/devel/DESCRIPTION.

I don't know how this update will impact CRAM support in Rsamtools though.

[vjcitn]

@pd3 please have a look at the cram branch noted in previous comments supply information on the use case and errors encountered thanks

[vjcitn]

As noted some time ago, our core has not encountered significant demand for CRAM support. It would be logical to add such to Rsamtools, but efforts to do this have not come to closure. I followed instructions at http://www.htslib.org/workflow/cram.html to make a fresh example, but ran into errors of this nature:

[E::easy_errno] Libcurl reported error 77 (Problem with the SSL CA cert (path? access rights?))
[W::find_file_url] Failed to open reference "https://www.ebi.ac.uk/ena/cram/md5/6681ac2f62509cfc220d78751b8dc524": Input/output error

That crops up while exploring a role for pysam in providing CRAM support. I'd really like to see some community input on use cases and good example files to make a case for a solution, which I suspect could involve reticulate/basilisk with pysam.

[DaGaMs]

Just nudging this: CRAM support is going to become more important as several new sequencing machines produce CRAM by default (Ultima, Axelios, maybe others). Converting a 1TB CRAM file to a bam is a nightmare, so native cram input would be extremely useful for us...

[vjcitn]

I am looking at the update of Rhtslib to use htslib 1.23.1. When that is done we will reinvestigate issues related to CRAM in Rsamtools. Do you have references to any manageable example datasets from the platforms you mention?

[hpages]

Not familiar with the CRAM format or with the differences between BAM and CRAM, but I suspect that adding support for CRAM in Rsamtools will require adding a set of tools similar to the set of tools that is currently available for operating on BAM files (BamFile objects, scanBam(), ScanBamParam(), sortBam(), etc...). This will be a major undertaking.

@vjcitn Rhtslib has htslib 1.18 which already exposes the CRAM API via the htslib/cram.h header file. I don't see a lot of differences in htslib/cram.h between htslib 1.18 and htslib 1.23.1 so I don't know if an update is actually necessary to start supporting CRAM in Rsamtools. My impression is that the effort should probably be on Rsamtools at this point.

[mtmorgan]

FWIW from my earlier effort (4+ years ago) I got bogged down in trying to implement support for everything in current Rsamtools. I kind of remember the branch that @hpages mentioned getting most of the core functionality working on CRAM, but getting bogged down in pileup and similar. So I'd suggest creating a plan that enables CRAM support for core functionality, with subsequent incremental additions. Or simply starting again with an Rcramtools to make a cleaner break from the accumulated technical debt.

[hpages]

Thanks Martin. I like the Rcramtools proposal. Seems more appealing than trying to stuff the already big Rsamtools package with the new CRAM code. Doesn't even need to be a BioC core package. @DaGaMs @pd3 Want to give it a try? I wonder if this could be a good candidate for the BioC2026 General Hackathon.

[DaGaMs]

I'm actually somewhat surprised that you found CRAM files to cause problems with any of the htslib functions. We built rastair around htslib, which required both read streaming and pileup, and the only "change" needed for CRAMs was to pass an indexed reference fasta[.gz]. Just glancing over some code, I think it's just a matter of calling hts_set_fai_filename once on the htsFile, and after that everything else should just behave "as usual"?

[hpages]

I think it's just a matter of calling hts_set_fai_filename once on the htsFile, and after that everything else should just behave "as usual"?

Good news then. Want to submit a PR @DaGaMs ?