Extract bed regions from CRAM files

Question

[MAPK]

Hi Everyone,

I have several hundreds of cram files. I need to extract some bed regions from these files. Is there a way to extract bed regions from cram without unwrapping them to BAM? My final output should be BAM file for selected regions for each sample.

Answer 1

[ATpoint]

In general the -L option works for CRAM (according to the manual) same as it does for BAM given you indexed them with samtools index.

samtools view --reference ref.fa -L your.bed -o overlaps.bam your.cram

Answer 2

[Pierre Lindenbaum]

 samtools view -b --reference ref.fa in.cram | bedtools bamtobed -i stdin

[MAPK]

Hi Pierre, I am starting with CRAM not BAM.

[Pierre Lindenbaum]

Hi Peter, I am starting with CRAM not BAM.

fixed

[ATpoint]

Added -b to your command to avoid unheadered SAM output :)

[Pierre Lindenbaum]

ah yes ! thanks you ! :-)

[GenoMax]

How is the solution addressing this requirement?

My final output should be BAM file for selected regions for each sample.

[MAPK]

Not sure why this gives me error, but @ATpoint's answer works great.

 [an@user]$ samtools view -b --reference /home/all_sequences.fa /40/AD/XX^92XX^0319.cram | bedtools bamtobed -i genes_build37.bed-hg38
    [W::sam_read1] Parse error at line 1
    Failed to open BAM file ades_genes_build37.bed-hg38