Question: Add read group to header after samtools merge
1
gravatar for dariober
4.4 years ago by
dariober10k
WCIP | Glasgow | UK
dariober10k wrote:

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

read group samtools merge • 8.0k views
ADD COMMENTlink modified 4.4 years ago • written 4.4 years ago by dariober10k
4
gravatar for dariober
4.4 years ago by
dariober10k
WCIP | Glasgow | UK
dariober10k wrote:

I ended up writing a little script for this, if anyone is interested is here addRGtoSAMHeader.py.

ADD COMMENTlink modified 4.0 years ago • written 4.4 years ago by dariober10k
6
gravatar for Pierre Lindenbaum
4.4 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum120k wrote:

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 http://broadinstitute.github.io/picard/command-line-overview.html#MergeSamFiles

ADD COMMENTlink written 4.4 years ago by Pierre Lindenbaum120k

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

ADD REPLYlink written 4.4 years ago by dariober10k

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

ADD REPLYlink written 4.4 years ago by Pierre Lindenbaum120k
1

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

 

 

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by dariober10k

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

ADD REPLYlink written 4.0 years ago by Jordan1.1k

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

ADD REPLYlink written 4.0 years ago by dariober10k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1104 users visited in the last hour