#! /usr/bin/perl 
# delete the INDELs that convered by another one.

# author: X.Gan

use strict;
use Getopt::Std;
use File::Basename;
use List::Util qw[min max];

my %opts =( S=>undef, d=>undef);
getopts('Sd', \%opts);
my $isSam =$opts{S} ;
my $deleteRG= $opts{d};

die "Usage: split_bylib.pl bamfile prefix \n" if (@ARGV <2); 

my $bamfile =shift;
my $prefix  =shift;
my $SMname=basename($prefix);

my $k;
if ( $isSam) {
    open(BAM_F, $bamfile) or die "Can't find the file";
} else {
    my $cmd="samtools view -h $bamfile | ";
    open(BAM_F, $cmd) or die "Can't find the file";
}

my @sticky;
my %ct;
my $secline;
my %RGlist;  ## map RG->LB
my %LBlist;  ## map LB->filehand
my %RG2file; 
my @filelist;
while (  <BAM_F>  ) 
{
    $secline=$_;
    last if ( $secline !~ /^\@/ );
    if ($secline =~ /^\@RG/ ) {
	my ($LB)= $_ =~ /^.*LB:(\S+).*/;
	my ($RG)= $_ =~ /^.*ID:(\S+).*/;
	$LBlist{$LB}=1;
	$RGlist{$RG}= $LB;
    }
}

foreach my $LB (keys %LBlist) {
    #print "\t", $`;
    my $myfilename="$prefix"."_".$LB.".sam";
    my $tf;
    open($tf, ">", $myfilename);
    $LBlist{$LB}=$tf;            ##put filehandle for each library
    push(@filelist, $myfilename);
}
foreach my $RG (keys %RGlist) {
    $RG2file{$RG} = $LBlist{ $RGlist{$RG} }; #get filehand for each RGgroup (share the same for the same lib);
}


############# processing the header ###########
close(BAM_F);
###seek(BAM_F, 0, 0);
if ( $isSam) {
    open(BAM_F, $bamfile) or die "Can't find the file";
} else {
    my $cmd="samtools view -h $bamfile | ";
    open(BAM_F, $cmd) or die "Can't find the file";
}


my $printRG=0;    ##Has already printed RG header?
if ($deleteRG) {
	$printRG=1;
}
my $printCO=1;
my $printPG=0;
while ( <BAM_F>  ) 
{	
    $secline=$_;
    last if ( $secline !~ /^\@/ );
    if ($secline =~ /^\@RG/ ) {
	next if ( $printRG );
	foreach my $LB (keys %LBlist) {
	    my $fh=$LBlist{$LB};
	    my $lineLB="\@RG\tID:$LB\tCN:WTCHG\tLB:$LB\tPL:ILLUMINA\tSM:$SMname\n";
	    print {$fh} $lineLB;
	}
	$printRG=1;
	next;  ###don;t print anymore
    }

    if ($secline =~ /^\@CO/) {
	next if ($printCO ); 
	$printCO=1;
    }
    if ( $secline =~ /^\@PG/ ) {
	next if ($printPG);
	$secline="\@PG\tID:IMR_STAMPY\tVN:0.8\n";
	$printPG=1;
    }
    foreach my $LB (keys %LBlist) {
	my $fh=$LBlist{$LB};
	print {$fh} $secline;
    }
}
if ($deleteRG) {
    while ($secline) {
	chomp $secline;
	if ($secline =~ s/[[:space:]]RG:Z:(\S+)//g) {
	    my $RG=$1;
	    my $fh=$RG2file{$RG};
	    print {$fh} $secline, "\n";
	}
	$secline = <BAM_F>;
    }

} else {
    while ($secline) {
	chomp $secline;
	if ($secline =~ /RG:Z:(\S+)/) {
	    my $RG=$1;
	    my $newID=$RGlist{$RG};
	    $secline =~ s/RG:Z:(\S+)/RG:Z:$newID/g; 
	    my $fh=$RG2file{$RG};
	    print {$fh} $secline, "\n";
	}
	$secline = <BAM_F>;
    }
}



foreach (@filelist ) {
    my $samfile=$_;
    my $bamfile =$samfile;
    $bamfile =~ s/.sam$/.bam/g;
    my $cmd="samtools view -b -S $samfile >$bamfile";
    print $cmd, "\n";
    system("$cmd");

}



exit 1;
