SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
GATK sample/library/lane meaning in BAM read group @RG Sylphide Bioinformatics 6 05-27-2014 08:20 AM
converting bam files to non-normalized read counts lpn Bioinformatics 4 10-09-2012 07:52 PM
adding platform information to a bam file efoss Bioinformatics 2 09-12-2011 10:22 AM
Consensus part from sequence read(fastq) and align(BAM) files culmen Bioinformatics 5 12-21-2010 03:57 AM
Adding read group and Platform information sbaheti Bioinformatics 2 09-25-2010 10:47 AM

Reply
 
Thread Tools
Old 02-26-2010, 06:57 AM   #1
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Question Adding Read Group info to a set of Bam files

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.
wjeck is offline   Reply With Quote
Old 02-26-2010, 07:15 AM   #2
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

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
lh3 is offline   Reply With Quote
Old 02-26-2010, 07:20 AM   #3
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

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.

Last edited by wjeck; 02-26-2010 at 07:28 AM.
wjeck is offline   Reply With Quote
Old 03-09-2010, 05:26 AM   #4
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

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.
wjeck is offline   Reply With Quote
Old 03-10-2010, 04:33 AM   #5
javijevi
Member
 
Location: Spain

Join Date: Jan 2010
Posts: 38
Default

Quote:
Originally Posted by lh3 View Post
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?
javijevi is offline   Reply With Quote
Old 03-10-2010, 07:04 AM   #6
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

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
wjeck is offline   Reply With Quote
Old 03-10-2010, 07:34 AM   #7
javijevi
Member
 
Location: Spain

Join Date: Jan 2010
Posts: 38
Default

Quote:
Originally Posted by wjeck View Post
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?
javijevi is offline   Reply With Quote
Old 03-10-2010, 12:00 PM   #8
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

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
wjeck is offline   Reply With Quote
Old 06-03-2010, 10:21 PM   #9
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

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!
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Old 07-22-2010, 06:22 AM   #10
Brugger
Member
 
Location: Cambridge, UK

Join Date: Mar 2010
Posts: 21
Default

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.
Brugger is offline   Reply With Quote
Old 08-02-2010, 12:29 AM   #11
allPower
Junior Member
 
Location: Australia

Join Date: Aug 2010
Posts: 1
Default Adding readgroups to one .bam file

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

Samtools throws an error

Code:
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) ?
allPower is offline   Reply With Quote
Old 09-15-2010, 05:58 AM   #12
freeseek
Junior Member
 
Location: Cambridge MA

Join Date: Feb 2010
Posts: 8
Default Re: Adding readgroups to one .bam file

@Michael.James.Clark the following two lines of bash code:
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.
freeseek is offline   Reply With Quote
Old 09-15-2010, 06:12 AM   #13
Brugger
Member
 
Location: Cambridge, UK

Join Date: Mar 2010
Posts: 21
Default

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.

Code:
#!/usr/bin/perl
#
# Add or replace a readgroup/sample/library/platform information to a samfile.
#
#
# Kim Brugger (22 Jul 2010), contact: kim.brugger@easih.ac.uk

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;
}
Brugger is offline   Reply With Quote
Old 10-03-2010, 12:56 AM   #14
biaoluo
Junior Member
 
Location: Philadelphia, PA

Join Date: Jan 2009
Posts: 2
Default

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/wi...ate_read_group
biaoluo is offline   Reply With Quote
Old 10-06-2010, 10:09 PM   #15
caddymob
Member
 
Location: USA

Join Date: Apr 2009
Posts: 36
Default

biaoluo, thanks for this! I feel like a dummy, but I cannot figure out how to apply the patch. Can you share how?
caddymob is offline   Reply With Quote
Old 10-06-2010, 10:23 PM   #16
caddymob
Member
 
Location: USA

Join Date: Apr 2009
Posts: 36
Default How to patch BWA

Ok, I figured it out! Here is how I did it...

go into the SVN checkout of bio-bwa/trunk/bwa, then run this:

Code:
patch bwape.c BWA_read_group_patch.diff
then:

Code:
make
then test it:

Code:
./bwa sampe

Usage:   bwa sampe [options] <prefix> <in1.sai> <in2.sai> <in1.fq> <in2.fq>

Options: -a INT   maximum insert size [500]
         -o INT   maximum occurrences for one end [100000]
         -n INT   maximum hits to output for paired reads [3]
         -N INT   maximum hits to output for discordant pairs [10]
         -c FLOAT prior of chimeric rate (lower bound) [1.0e-05]
         -f FILE sam file to output results to [stdout]

         -P       preload index into memory (for base-space reads only)
         -s       disable Smith-Waterman for the unmapped mate
         -A       disable insert size estimate (force -s)

         -i       read group identifier (ID)
         -m       read group sample (SM), required if ID is given
         -l       read group library (LB)
         -p       read group platform (PL)
Notes: 1. For SOLiD reads, <in1.fq> corresponds R3 reads and <in2.fq> to F3.
       2. For reads shorter than 30bp, applying a smaller -o is recommended to
          to get a sensible speed at the cost of pairing accuracy.
the -i, -m, -l and -p options are the ticket!
caddymob is offline   Reply With Quote
Old 04-14-2011, 10:14 AM   #17
jyli
Member
 
Location: North Carolina

Join Date: Nov 2008
Posts: 21
Default

Quote:
Originally Posted by wjeck View Post
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.
Were you ever be able to figure it out (with the already a merged bam file)?
jyli is offline   Reply With Quote
Old 04-14-2011, 10:37 AM   #18
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

The way to do this now is to use the Picard command line tool, in the latest picard version.
wjeck is offline   Reply With Quote
Old 04-15-2011, 04:53 AM   #19
Seq84
Member
 
Location: Italy

Join Date: Feb 2011
Posts: 19
Default

Quote:
Usage: bwa sampe [options] <prefix> <in1.sai> <in2.sai> <in1.fq> <in2.fq>

Options: -a INT maximum insert size [500]
-o INT maximum occurrences for one end [100000]
-n INT maximum hits to output for paired reads [3]
-N INT maximum hits to output for discordant pairs [10]
-c FLOAT prior of chimeric rate (lower bound) [1.0e-05]
-f FILE sam file to output results to [stdout]
-r STR read group header line such as `@RG\tID:foo\tSM:bar'[null]
-P preload index into memory (for base-space reads only)
-s disable Smith-Waterman for the unmapped mate
-A disable insert size estimate (force -s)

Notes: 1. For SOLiD reads, <in1.fq> corresponds R3 reads and <in2.fq> to F3.
2. For reads shorter than 30bp, applying a smaller -o is recommended to
to get a sensible speed at the cost of pairing accuracy.
This is BWA sampe in 0.5.9-r16 version.
With -r option followed by such kind of string you could insert RG directly during mapping.

For editing RG lines use AddOrReplaceReadGroup in Picard.
Seq84 is offline   Reply With Quote
Old 04-15-2011, 04:59 AM   #20
wjeck
Member
 
Location: Chapel Hill, NC

Join Date: Mar 2009
Posts: 39
Default

I can attest that both of these tools work well. The PICARD tool in particular is vastly quicker than my previous workaround.
wjeck is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:04 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO