Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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.

  • #2
    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:

    Download SAM tools for free. SAM (Sequence Alignment/Map) is a flexible generic format for storing nucleotide sequence alignment. SAMtools provide efficient utilities on manipulating alignments in the SAM format.

    Comment


    • #3
      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, 08:28 AM.

      Comment


      • #4
        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.

        Comment


        • #5
          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?

          Comment


          • #6
            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

            Comment


            • #7
              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?

              Comment


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

                Comment


                • #9
                  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]

                  Comment


                  • #10
                    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.

                    Comment


                    • #11
                      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) ?

                      Comment


                      • #12
                        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.

                        Comment


                        • #13
                          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: [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;
                          }

                          Comment


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

                            Comment


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

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Exploring the Dynamics of the Tumor Microenvironment
                                by seqadmin




                                The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                                07-08-2024, 03:19 PM
                              • seqadmin
                                Exploring Human Diversity Through Large-Scale Omics
                                by seqadmin


                                In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                                06-25-2024, 06:43 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, Today, 06:53 AM
                              0 responses
                              11 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-10-2024, 07:30 AM
                              0 responses
                              32 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-03-2024, 09:45 AM
                              0 responses
                              203 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-03-2024, 08:54 AM
                              0 responses
                              213 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X