samtools merge respecting original read groups

Question

[abascalfederico]

Hi,

I'm struggling to merge many bam files respecting the original read groups with samtools. Some of the input files have the same read group, other files have others. If I don't use "-c" each "duplicate" RGs are assigned a random suffix (which I don't want to, they are the same RG). If I use "-c", RGs are not altered, but then only one RG line is printed in the header.

Does any one have a clue on how to solve this?

Ideally, with "-c" samtools should print any RG header line that is different from the others in the merged bam.

Many thanks Federico

[Devon Ryan]

Which version of samtools are you using? I recall that there were modifications to handle this relatively recently.

[abascalfederico]

Thanks for the feedback, Devon. I am using Version: 1.2 (using htslib 1.2.1). I have looked at the latest version (1.3) but don't see anything about this.

I've read in other post that for this kind of situation picard-tools is the right way to go. Wish samtools change this.

Answer 1

[abascalfederico]

In case it is useful for someone, I've written a perl routine to solve this problem.

sub create_header {
my $p_input_bams = $_[0];
my $header_file  = $_[1];
my $RG_string = "";
my %rgs;
foreach my $b ( @{$p_input_bams} ) {
    open(J,"samtools view -H $b|") || die "cñakjdñakj $b\n";
    while(<J>) {
        #print STDERR "$_";
        last unless(/^\@/);
        if(/^\@RG/) {
            $_ =~ s/\#[0-9]+//;
            $rgs{$_} = 1;
        }
    }
    close(J);
}
foreach my $rg ( keys %rgs ) {
    $RG_string .= $rg;
}
print STDERR "RG_string=\n$RG_string\n";
my $main_bam = $p_input_bams->[0];
open(IB, "samtools view -h $main_bam |") || die "Cannot open main: $main_bam\n";
open(OB, ">$header_file") || die "cannot write to header: $header_file\n";
while(<IB>) {
    last unless(/^\@/);
    if(!/^\@RG/) {
        print OB;
    } else {
        print STDERR "Printing RG string...\n";
        print OB $RG_string;
    }

}
close(OB);
close(IB);

}

To use it:

system("samtools merge -cf $lib.merge-changeHeader.bam @input_bams");
&create_header(\@input_bams,$header_file);
system("samtools reheader $header_file $lib.merge-changeHeader.bam > $lib.merge.bam");

[GenoMax]

I assume this part is lost in unicode translation cñakjdñakj :-)

[abascalfederico]

hehehe my "die" error messages are usually random typings. It's fast and helps a lot debugging errors