Seqanswers Leaderboard Ad
Collapse
Announcement
Collapse
No announcement yet.
X
-
biaoluo, thanks for this! I feel like a dummy, but I cannot figure out how to apply the patch. Can you share how?
-
The GATK team at Broad Institute have just provided a BWA patch that can generate read group information during BWA alignment:
Leave a comment:
-
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; }
Leave a comment:
-
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
Leave a comment:
-
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) ?
Leave a comment:
-
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.
Leave a comment:
-
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!
Leave a comment:
-
Yes, that's what I meant. You can find the definitions of the two letter tag codes themselves in the sam specification:
Leave a comment:
-
Originally posted by wjeck View PostI 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
Leave a comment:
-
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
Leave a comment:
-
Originally posted by lh3 View PostYou 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
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?
Leave a comment:
-
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.
Leave a comment:
-
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:
Leave a comment:
-
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.Tags: None
Latest Articles
Collapse
-
by seqadmin
Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...-
Channel: Articles
04-04-2024, 04:25 PM -
-
by seqadmin
Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...-
Channel: Articles
03-22-2024, 06:39 AM -
ad_right_rmr
Collapse
News
Collapse
Topics | Statistics | Last Post | ||
---|---|---|---|---|
Started by seqadmin, 04-11-2024, 12:08 PM
|
0 responses
24 views
0 likes
|
Last Post
by seqadmin
04-11-2024, 12:08 PM
|
||
Started by seqadmin, 04-10-2024, 10:19 PM
|
0 responses
25 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 10:19 PM
|
||
Started by seqadmin, 04-10-2024, 09:21 AM
|
0 responses
21 views
0 likes
|
Last Post
by seqadmin
04-10-2024, 09:21 AM
|
||
Started by seqadmin, 04-04-2024, 09:00 AM
|
0 responses
52 views
0 likes
|
Last Post
by seqadmin
04-04-2024, 09:00 AM
|
Leave a comment: