#!/usr/bin/perl -w

# usage : perl bwasw.pl <FASTA ref file> -b1 <FASTQ file> <output SAM file> 
#sarah maman - 18/04/2013
# Copyright (C) 2013 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_fastq        = $ARGV[1];
my $output_sam         = $ARGV[2];
my $tailleFichier      = (stat($input_ref_genome))[7];
#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';
my $cmd = '';

#Comme les users qui envoient les jobs ne peuvent pas ecrire dans /galaxydata/database/files/00X, les donnees sont placees dans le workspace
my ($nb) = ($output_sam=~/galaxy_dataset_(\d+)\.\S+$/);
`cp $input_ref_genome /galaxydata/database/files/workspace/ref_$nb.fasta; cp $input_fastq /galaxydata/database/files/workspace/input_$nb.fastq; `;
my $ref ="/galaxydata/database/files/workspace/ref_$nb.fasta";
my $fastq ="/galaxydata/database/files/workspace/input_$nb.fastq";

#Indexation du genome de reference (ne pas placer cette commande en background)
#si le genome de référence est gros , plus de 10 MB, alors indexer avec l'option -a bwtsw
#si le genome est petit, indexer avec -a is 
if ($tailleFichier > 10485760) {
        $cmd = "$BWA index -a bwtsw $ref  >> ./bwaindex.log 2>&1";
      }
if ($tailleFichier < 2147483648){
        $cmd = "$BWA index -a is $ref  >> ./bwaindex.log 2>&1";
      }
else {
       $cmd ="$BWA index $ref  >> ./bwaindex.log 2>&1";
     }
system $cmd;

#Information pour les biologistes
print STDOUT "Etape 1 - Indexation : $cmd \n\n";

#Envoyer un message d'erreur et arreter l envoi de lignes de commandes si l indexation s est mal deroulee
if ($? != 0){print STDERR "Incident lors de l indexation. Fin du job. Veuillez relancer votre outil"; exit (1);}

#Tester la commande :
# perl sm_bwasw.pl tests/cow.fasta tests/long.fastq toto
#Alignement des longs reads
my $sam = "/galaxydata/database/files/workspace/output_$nb.sam";
system ("($BWA bwasw $ref $fastq > $sam) >& ./bwasw.log 2>&1");


#Information pour les biologistes
print STDOUT "Etape 2 - BWA : $BWA bwasw $ref $fastq > $sam \n";

if (! -e "$sam")
{
print STDERR "SAM FILE NOT FOUND\n";
}
else
{
`cp -a "$sam" $output_sam ;`
}

