#!/usr/bin/perl -w

# usage : perl bwa_aln_sampe.pl <FASTA ref file> -b1 <BAM file> <output BAM file> 
# sarah maman - 10 juillet 2012 - Mis en place dans le cadre du projet SwanPorc
# Copyright (C) 2012 INRA
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

use strict;
use File::Basename;

my $input_ref_genome = $ARGV[0];
my $input_bam = $ARGV[1];
my $output_bam = $ARGV[2];
my ($nb) = ($output_bam=~/galaxy_dataset_(\d+)\.\S+$/);
my $input_dir=dirname($input_bam);
my $input_basename=basename($input_bam,'.bam');
my $input_bamsorted="$input_dir/${input_basename}-sorted";
#pour eviter que la commande ne fonctionne plus si le $PATH est perdu suite a un restart de galaxy
my $BWA = '/usr/local/bioinfo/bin/bwa';
#pour connaître le path des samtools, utiliser la commande which samtools
my $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my $cmd = '';

#Ouverture des fichiers entrants
open (IN, "<$input_ref_genome") or die "Cannot open $input_ref_genome !";
open (IN, "<$input_bam") or die "Cannot open $input_bam !";

#Tri par nom du bam entrant pour travailler sur des donnees pairees (commande pas en bg)
$cmd = "($SAMTOOLS sort -n $input_bam $input_bamsorted ) >& ./samtools.log 2>&1";
system $cmd;
if (! -e "${input_bamsorted}.bam"){print STDERR "STEP 1 - BAM SORTED NOT FOUND\n";}
else{`(cp $input_bamsorted ${input_bamsorted}-sorted.bam) >& ./bwa.log 2>&1`}
#print "\n1 - TRI\n";

#Alignement des donnees pairees uniquement
system ("($BWA aln $input_ref_genome -b1 ${input_bamsorted}.bam > ${input_bamsorted}-1.sai) >& /galaxydata/logs/sm_bwaaln1_$nb.log 2>&1");
#print "2 - Alignement\n";
#print "($BWA aln $input_ref_genome -b1 ${input_bamsorted}.bam > ${input_bamsorted}-1.sai) >& /galaxydata/logs/sm_bwaaln.log 2>&1";
#print "\n\n";
#exit();

#$cmd = "($BWA aln $input_ref_genome -b1 ${input_bamsorted}.bam ) >& ${input_bamsorted}-1.sai >> ./bwaaln.log 2>&1";
#$cmd = "($BWA aln $input_ref_genome -b1 ${input_bamsorted}.bam )>& ./bwaaln.log 2>&1";
#system $cmd;
if (! -e "${input_bamsorted}-1.sai")
{
print STDERR "SAI FILE1 NOT FOUND\n";
}

$cmd = "($BWA aln $input_ref_genome -b2 ${input_bamsorted}.bam > ${input_bamsorted}-2.sai ) >& /galaxydata/logs/sm_bwaaln2_$nb.log 2>&1";
#$cmd = "$BWA aln $input_ref_genome -b2 ${input_bamsorted}.bam >> ./bwaaln.log 2>&1";
system $cmd;
if (! -e "${input_bamsorted}-2.sai")

{
print STDERR "SAI FILE2 NOT FOUND\n";
}


#bwa sampe
$cmd = "($BWA sampe $input_ref_genome ${input_bamsorted}-1.sai ${input_bamsorted}-2.sai ${input_bamsorted}.bam ${input_bamsorted}.bam | $SAM
TOOLS view -bS - > ${input_bamsorted}-sampe.bam) >& /galaxydata/logs/sm_bwa3_$nb.log 2>&1";
#print "3 - BWA SAMPE";
#print "($BWA sampe $input_ref_genome ${input_bamsorted}-1.sai ${input_bamsorted}-2.sai ${input_bamsorted}.bam ${input_bamsorted}.bam | $SAM
#TOOLS view -bS - > ${input_bamsorted}-sampe.bam) >& ./bwa.log 2>&1";
#$cmd = "$BWA sampe $input_ref_genome $output_temp_1 $output_temp_2 ${input_bamsorted}.bam ${input_bamsorted}.bam | $SAMTOOLS view -bS - >> 
#./bwa.log 2>&1";

system $cmd;

if (! -e "${input_bamsorted}-sampe.bam")
{
print STDERR "FINAL BAM FILE NOT FOUND\n";
}
else
{
#system ("`cp -a ${input_bamsorted}-sampe.bam $output_bam`");
#system ("${input_bamsorted}-sampe.bam > $output_bam");
#print "cp -a \"${input_bamsorted}-sampe.bam\" $output_bam";
`(cp -a \"${input_bamsorted}-sampe.bam\" $output_bam) >& ./bwa.log 2>&1`
}

close( IN );
