Adding Read Group info to a set of Bam files

[wjeck]

So I have a group of 4 or 5 bam files created with BWA sitting nice and spiffy in a directory. I need to merge them for later analysis, but I want to make sure that the reads are appropriately tagged with readgroup information that gets preserved as I go forward. Does anyone know an easy way to add readgroup info to these files individually before the merge? As I understand it there is more than just the @RG tag...

What I'd like to do is run a command line utility like this:

java -jar addReadGroup.jar -i fileIN.bam -o fileOUT.bam -r "cellLineX_readgroup1"

It didn't look like Picard has such a program, but maybe I'm missing something.

[lh3]

You may try "samtools merge", using options -r and -h. You write your @RG header lines in a file provided to -h; -r will add RG:Z: tag to each of the alignment, based on file names.

EDIT: for an example:

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_protocol#Generate_merged_alignment

[wjeck]

Perfect! I should have read the readme more closely. Many thanks!

PS: I note the -r requires the latest SAMTOOLS 1.07. Just for any future readers of the thread. Great addition, though.

[wjeck]

Follow up question: Is there a way to edit the information in the @RG tag after the files have been merged in BAM format? I'd like to add and subtract information from these lines downstream, and I can't figure out an elegant way to get into them without writing out an entire SAM file and translating it back to BAM.

[javijevi]

You may try "samtools merge", using options -r and -h. You write your @RG header lines in a file provided to -h; -r will add RG:Z: tag to each of the alignment, based on file names.

EDIT: for an example:

http://sourceforge.net/apps/mediawik...rged_alignment

In this wike, one can found the following commands:

perl -e 'print "@RG\tID:ga\tSM:hs\tLB:ga\tPL:Illumina\n@RG\tID:454\tSM:hs\tLB:454\tPL:454\n"' > rg.txt
samtools merge -rh rg.txt - ga.bam 454.bam | samtools rmdup - - | samtools rmdup -s - aln.bam

Does anybody could tell what does exactly mean the ID, SM, LB and PL tags at
"@RG\tID:ga\tSM:hs\tLB:ga\tPL:Illumina\n@RG\tID:454\tSM:hs\tLB:454\tPL:454\n", so that we could adapt them to our files?

[wjeck]

I think, tough this is not with any authority, that

ID = id name for the readgroup
SM = sample name
LB = label? dunno about this one
PL = platform

These are not currently standardized (I think) but ARE used by the Broad GATK, which means getting them right may be important for your pipeline

[javijevi]

I think, tough this is not with any authority, that

ID = id name for the readgroup
SM = sample name
LB = label? dunno about this one
PL = platform

These are not currently standardized (I think) but ARE used by the Broad GATK, which means getting them right may be important for your pipeline

Thanks a lot. By "these are not currenly standardized" do you mean that values for this tags are not pre-defined (e.g., '454', 'Illumina', and 'solid' for PL tag) so that they can be freely selected?

[wjeck]

Yes, that's what I meant. You can find the definitions of the two letter tag codes themselves in the sam specification:

http://samtools.sourceforge.net/SAM1.pdf

[Michael.James.Clark]

Trying to do something similar here.

Is there a way to add RG to a single BAM file without merging? Like, to add the RG to a BAM file that doesn't have a RG flag at all?

I'm trying to run samtools merge with just one file, but it won't work, so I made an empty dummy file, but that also didn't work (gave a segmentation fault).

Thanks!

[Brugger]

Bunked into the same problem, and after googling for day or so it is clear that no one has posted a solution.

I have a perl script that will take a sam-file/stream that will add the RG information if anyone wants it.

[allPower]

Adding readgroups to one .bam file

I also tried and failed to add a RG using both samtools and picard.

Samtools throws an error

x="sample4"
perl -e 'print "\@RG\tID:Disc1\tSM:'$x'\tPL:Illumina\n"' >rg.txt
samtools merge -rh rg.txt - index10.aln.sort.bam

Usage:   samtools merge [-nr] [-h inh.sam] <out.bam> <in1.bam> <in2.bam> [...]

Options: -n       sort by read names
         -r       attach RG tag (inferred from file names)
         -h FILE  copy the header in FILE to <out.bam> [in1.bam]

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.

from as suggested here: http://sourceforge.net/apps/mediawik...rged_alignment (note the additional "\" in my version at the beginning otherwise it loses @RG...)

Picard's MergeBamAlignment wants at least two bams and ReplaceSamHeader does not add the RG to the reads.

Does anyone have a good solution for this (Brugger) ?

[freeseek]

Re: Adding readgroups to one .bam file

@Michael.James.Clark the following two lines of bash code:

echo -e "@RG\tID:ga\tSM:hs\tLB:ga\tPL:Illumina" > rg.txt
samtools view -h ga.bam | cat rg.txt - | awk '{ if (substr($1,1,1)=="@") print; else printf "%s\tRG:Z:ga\n",$0; }' | samtools view -uS - | samtools rmdup - - | samtools rmdup -s - aln.bam

should add to the bam file the read group information in the same way samtools merge adds the read group information to the two bam files as described by javijevi. The idea is to unpack the bam file, add the read group header, add the read group information to every read, repack the file, and remove duplicates. Again, remove duplicates only if the coverage is not too deep.

[Brugger]

My more complicated perl script that adds readgroup information, but can also replace read group information. It takes sam as input and gives the same back. run the program without any parameters to get the usage message.

#!/usr/bin/perl
#
# Add or replace a readgroup/sample/library/platform information to a samfile.
#
#
# Kim Brugger (22 Jul 2010), contact: [email protected]

use strict;
use warnings;
use Data::Dumper;
use Getopt::Std;
use File::Temp qw/ tempfile /;


my %opts;
getopts('i:o:r:s:l:p:c:a:A:R:', \%opts);
usage() if ( $opts{h});

my $infile    = $opts{i};
my $outfile   = $opts{o};
my $readgroup = $opts{r} || usage();
my $sample    = $opts{s} || $readgroup;
my $library   = $opts{l} || $readgroup;
my $platform  = $opts{p} || usage();
my $center    = $opts{c} || "";

my $aligner   = $opts{a};
my $a_line    = $opts{A};

my $replace   = $opts{R};
$infile = $replace if ( ! $infile && $replace );
open (*STDIN, $infile) || die "Could not open '$infile': $!\n" if ( $infile );

my ($tmp_fh, $tmp_file);
if ( $replace) {
  ($tmp_fh, $tmp_file) =  tempfile( DIR => './');
  *STDOUT = *$tmp_fh;
}
elsif ( $outfile ) {
  open (*STDOUT, " > $outfile " ) || die "Could not open '$outfile': $!\n" ;
}


my ($written_readgroup, $written_cmd_flags) = (0, 0);

while(<STDIN>) {
  
  if ( /^\@RG/) {
    my ( $same_readgroup, $same_sample, $same_library, $same_platform, $seq_center) = (0, 0, 0, 0);
    foreach my $field ( split("\t", $_)) {
      next if ($field =~ /^\@/);
      my ($key, $value) = split(":", $field);
      $same_readgroup++ if ( $field eq 'ID' && $value eq $readgroup);
      $same_sample++    if ( $field eq 'SM' && $value eq $sample);
      $same_library++   if ( $field eq 'LB' && $value eq $library);
      $same_platform++  if ( $field eq 'PL' && $value eq $platform);
      $seq_center++     if ( $field eq 'CN' && $value eq $center);
    }

#    $written_readgroup++ if ( $same_readgroup && $same_sample && $same_library && $same_platform );
    if ( $same_readgroup && $same_sample && $same_library && $same_platform && $seq_center ) {
      $written_readgroup++;
    }
    else {
      next;
    }
  }

  if ( ! /^\@/ ) {

    if ( $aligner && ! $written_cmd_flags ) {
      my @fields = ("\@PG", "ID:$aligner");
      push @fields, "CL:$a_line" if ( $a_line );
      print join("\t", @fields) . "\n";
      $written_cmd_flags++;
    }

    if (! $written_readgroup ) {
      print join("\t", "\@RG", "ID:$readgroup","SM:$sample","LB:$library","PL:$platform", "CN:$center") . "\n";
      $written_readgroup++;
    }

  
    if ( (/\tRG:Z:(\w+)\t/ || /\tRG:Z:(\w+)\Z/) &&
     (/\tSM:Z:(\w+)\t/ || /\tSM:Z:(\w+)\Z/)) {
      
      s/(.*\tRG:Z:)(.*?)(\t.*)/$1$readgroup$3/;
      s/(.*\tSM:Z:)(.*?)(\t.*)/$1$sample$3/;
      
      s/(.*\tRG:Z:)(.*?)\Z/$1$readgroup/;
      s/(.*\tSM:Z:)(.*?)\Z/$1$sample/;
    }
    elsif ( /\tRG:Z:(\w+)\t/ || /\tRG:Z:(\w+)\Z/ ) {
      chomp($_);
      s/(.*\tRG:Z:)(.*?)(\t.*)/$1$readgroup$3/;
      s/(.*\tRG:Z:)(.*?)\Z/$1$readgroup/;
      $_ .= "\tSM:$sample\n";
    }
    elsif ( /\tSM:Z:(\w+)\t/ || /\tSM:Z:(\w+)\Z/ ) {
      chomp($_);
      s/(.*\tSM:Z:)(.*?)(\t.*)/$1$sample$3/;
      s/(.*\tSM:Z:)(.*?)\Z/$1$sample/;
      $_ .= "\tRG:$readgroup\n";
    }
    else {
      chomp($_);
      $_ .= "\tRG:Z:$readgroup\tSM:Z:$sample\n";
    }
  }

  print;
  
}


if ( $replace) {
  close($tmp_fh);
  system "mv $tmp_file $infile";
}




#
#
#
# Kim Brugger (22 Jul 2010)
sub usage {

  $0 =~ s/.*\///;
  print "$0 adds/replaces the readgroup/library/sample/platform/center tags in a sam file/stream\n";
  print "$0 -i[nfile (or stdin] -o[utfile (or stdout)] -r[eadgroup] -s[ample] -l[ibrary] -p[latform] -c[enter] -a[ligner] -A[ligner param] -R[eplace infile with fixed file]\n";

  exit 1;
}

[biaoluo]

The GATK team at Broad Institute have just provided a BWA patch that can generate read group information during BWA alignment:

http://www.broadinstitute.org/gsa/wiki/index.php/BWA_patch_to_generate_read_group

[caddymob]

biaoluo, thanks for this! I feel like a dummy, but I cannot figure out how to apply the patch. Can you share how?