samtools merge error

Question

[Hatem Elshazly]

Hi there,

I'm trying to merge an output bam file of SHRiMP and an output bam file of bowtie using samtools merge but there is a problem because @SN is not identical between the two files, i tried to use -h to override the header but it didn't work.

Anyone knows how to work around this problem without adjusting the files manually?

In case anyone needed them, here are the files headers:

file1:

@HD     VN:1.0  SO:coordinate
@SQ     SN:Chr17        LN:75158596
@PG     ID:gmapper      VN:2.2.3        CL:/var/www/tavaxy/engines/galaxy/tool-data/SHRiMP_2_2_3/bin/gmapper-cs /var/www/tavaxy/files/inputFiles/0337d77e79d8d1268b8d728e3696adb1/step1/run1_lane1.csfasta.folder/fileContents.csfasta /var/www/tavaxy/files/inputFiles/0337d77e79d8d1268b8d728e3696adb1/bowtie_mapper/Chr17.fa.folder/fileContents.fa -N 4 -o 5 -h 80X

file2:

@HD     VN:1.3  SO:coordinate
@SQ     SN:1_2_193      LN:75
@SQ     SN:1_2_252      LN:75
@SQ     SN:1_2_1015     LN:75
@SQ     SN:1_2_1113     LN:75

Command:

samtools merge -h header.tmp sample3.bam sample1.bam sample2.bam
[bam_merge_core] different target sequence name: '1_2_193' != 'Chr17' in file 'sample2.bam'

Thanks in advance,
Shazly

Answer 1

[Devon Ryan]

These files aren't mergable with samtools or picard or any other tool that I know of. What you can do is write file 2's header to a file and then edit it to contain Chr17. You can then use process substitution to merge them into a single file:

cat header <(samtools view file1) <(samtools view file2) | samtools view -Sbo merged.bam -

Edit: Apparently using MERGE_SEQUENCE_DICTIONARIES=true in MergeSamFiles from Picard tools allows this (I just tested and confirmed this).

[Pierre Lindenbaum]

even with picard ? there is an option inMergeSamFiles called "MERGE_DICT" . (not tested)

[Devon Ryan]

It does in fact work (I've updated my answer). Nice to know.

[Devon Ryan]

Guess I've never tested that option. I'll give it a whirl.

[Hatem Elshazly]

Ok, I tried using Picard with the MERGE_SEQUENCE_DICTIONARIES=true option and it worked.

Thank you very much for replying.