samtools, convert BAM to CRAM

Question

[nrk_02]

Hello!

I have a BAM file that was aligned to Feb. 2009 GRCh37/hg19 Assembly, and I would like to convert it to a CRAM file. I use this command:

samtools view -C -T reference/hs37d5.fa -o file.cram file.bam

However, I get an error as:

[E::cram_get_ref] Failed to populate reference for id 0
[main_samview] failed to write the SAM header

Also, I know that my BAM file is sorted as the header returns -> SO:coordinate

What do you think can be the issue here?

[ATpoint]

What is the output of grep '^>' reference/hs37d5.fa and samtools view -H file.bam | grep '^@'SQ?

[nrk_02]

BAM file header format:

@HD     VN:1  SO:coordinate
@SQ     SN:chrM LN:1234
@SQ     SN:chr1 LN:1234
@SQ     SN:chr2 LN:1234
@SQ     SN:chr3 LN:1234
....
...
... 

Reference header format:

>1 dna:chromosome chromosome:GRCh37:1:1:1234:1
>2 dna:chromosome chromosome:GRCh37:2:1:1234:1
>3 dna:chromosome chromosome:GRCh37:3:1:1234:1
>4 dna:chromosome chromosome:GRCh37:4:1:1234:1
...
..

[ATpoint]

There is your answer, a mismatch in chromosome identifiers.

[nrk_02]

Thank you very much! I know that I can change the BAM file chr notation, but then I am still not sure how it should look after editing, should it just be "chromosome" instead of "chr" ?

[ATpoint]

I would rather try to find the exact reference that was used for alignment. It should match the bam file, so it must be >chr1, >chr2 etc. Don't manipulate the bam, this is probably way more prone to chaos than altering the reference file names.

[nrk_02]

Thank you very much for your replies so far! I am just curious and want ro learn: If I were to manipulate BAM to make a match, how would I do it?

[ATpoint]

You would need a conversion table that lists the chromosome name in the BAM file, and the matched name in the reference and then do some gsub operation to find replace them. I do not have code for it at hand, sorry.