Quantify gene expression from CRAM file

Question

[wenbinm]

Hello Biostars folks,

Do you know any tools that can quantify gene expression from aligned CRAM files in RNASeq? In the past I used featureCounts but it doesn't accept CRAM file. I am trying to quantify gene expression from the CRAM files downloaded from the Hartwig Medical Foundation cohort. Converting CRAM to BAM and then using featureCounts works but this takes a lot of storage space.

Thank you!

Wenbin

Answer 1

[ATpoint]

HTseq apparently can deal with CRAM (https://htseq.readthedocs.io/en/master/tutorials/bam_reader.html?highlight=cram) probably because it relies on pysam which relies on htslib which again is the native way to read CRAM. So you can try and implement that. Or, and this is my recommendation, you simply do the (relatively fast) cram2bam conversion with samtools and then run featureCounts on every single file. Then once done combine the individual featureCounts outputs into a single count matrix. This is trivial as basically it's a loop over all files, so the storage requirements are limited. Or you just spend a few bucks on a larger hard drive. Storage these days is really cheap. I would not reinvent the wheel for this. Just loop (or GNU parallel, few files at a time) over the files, let run overnight and be done.

[Ram]

EDIT: Never mind, in a previous answer, you stated that the .bai file is not required. Process substitution might work here, saving OP some storage space at the cost of slower computation - unless featureCounts has a problem with streaming BAMs.

Can featureCounts work with <(samtools view -b file.cram) or does it need an actual indexed BAM file?

[wenbinm]

Thanks ATpoint. I tried HTSeq but it's much slower than featurecounts. I guess I have no choice but using cram2bam