Add read group to header after samtools merge

Question

[dariober]

Hello,

I merged bam files using samtools merge -r out.bam in1.bam in2.bam ...

With -r option I got the RG tag for each read, as expected, good. However, the header of the merged bam does not have the @RG lines. So my question: Is there any off-the-shelf tools to add the @RG header lines once the reads have been tagged?

If not, the pipeline I have in mind goes along these lines:

• Scan the bam file and collect all the different RG tags
• Output the header of the bam file.
• Append the RG tags to the header the
• Reheader the original bam file.

Does it make sense?

Thanks

Dario

Answer 1

[dariober]

I ended up writing a little script for this, if anyone is interested it is here:

• https://code.google.com/p/bioinformatics-misc/source/browse/trunk/addRGtoSAMHeader.py
• https://github.com/dariober/bioinformatics-cafe/blob/master/addRGtoSAMHeader.py

Answer 2

[Pierre Lindenbaum]

Don't use samtools merge:

Note: Samtools' merge does not reconstruct the @RG dictionary in the header. Users
      must provide the correct header with -h, or uses Picard which properly maintains
      the header dictionary in merging.

but picard MergeSamFiles

[dariober]

Thanks! I'll keep it mind and use picard in the future.

[Pierre Lindenbaum]

and you can use http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups

[dariober]

Hi again Pierre. Sorry maybe I'm being thick. I'm not sure if/how to use AddOrReplaceReadGroup for this purpose as it complains that RG ID on SAMRecord not found in header. But that's what I want to do: Put the RG ID from the reads to the header.

java -jar ~/bin/picard.jar AddOrReplaceReadGroups I=rhh05082014_300ng.bam O=rhh05082014_300ng.new.bam LB=rhh05082014_300ng PL=illumina PU=mixed SM=rhh05082014 RGID=null
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 1, Read name NS500222:44:H15HDBGXX:1:11105:2110:3420, RG ID on SAMRecord not found in header: rhh088_A9_Pull1_300ng_Inp05082014.141215_NS500222_0044_H15HDBGXX.S1_L001.hg19.clean

[Jordan]

Did you find a solution for this? I'm having the same issue.

[dariober]

@Jordan My solution was to write the script linked in my answer below.