#!/usr/bin/perl -w

# usage : bwa puis tri
# Sarah Maman - 10 juillet 2012 - Mis en place dans le cadre du projet SwanPorc
# Copyright (C) 20XX 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_fa = $ARGV[1];

#pour les options
my $flag1 = $ARGV[2];
my $flag2 = $ARGV[3];
my $mismatch = $ARGV[4];

#pour les sorties
my $output_bam = $ARGV[5];
my $output_filter1 = $ARGV[6];
my $output_filter2 = $ARGV[7];
my ($nb) = ($output_bam=~/galaxy_dataset_(\d+)\.\S+$/);

my $input_dir=dirname($input_fa);
my $input_basename=basename($input_fa,'.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 = '';


#commande
#bwa aln bank/$i 3_redundancy/ALL.matrix.filter.tfa |bwa samse bank/$i 3_redundancy/ALL.matrix.filter.tfa | samtools view bS | samtools sort 4_annotation/$i.sort


$cmd = "($BWA aln $input_ref_genome  $input_fa  | $BWA samse $input_ref_genome - $input_fa |$SAMTOOLS view -bS - | $SAMTOOLS sort - ${input_bamsorted}.sort ) >& /galaxydata/logs/sm_bwa_annot_$nb.log 2>&1";
system $cmd;

if (! -e "${input_bamsorted}.sort.bam")
{
print STDERR "BWA failed\n";
}
else
{
`(cp -a \"${input_bamsorted}.sort.bam\" $output_bam) >& /galaxydata/logs/sm_bwa_annot2_$nb.log 2>&1`
}


#Alignement

#Premier fichier issu du filtre de l'alignement
$cmd = qq/($SAMTOOLS view ${input_bamsorted}.sort.bam | grep 'NM:i:[$mismatch]' | awk '\$2==$flag1||\$2==$flag2' | cut -f1,3 | sort -k1,1 > ${input_bamsorted}.filter1) >& filtrer_aln_fichier1.log 2>&1/;
system $cmd;

print("CMD : $cmd");

if (! -e "${input_bamsorted}.filter1")
{
print STDERR "Echec lors de la generation du premier fichier de filtre des alignements. \n";
}
else
{
`(cp -a \"${input_bamsorted}.filter1\" $output_filter1) >& ./cp_filtrer_aln_fichier1.log 2>&1`
}

#Second fichier issu du filtre de l'alignement
$cmd = "(cut -f1 ${input_bamsorted}.filter1 > ${input_bamsorted}.filter2) >& /galaxydata/logs/sm_filtrer_aln_fichier1_$nb.log 2>&1";
system $cmd;
#print ("****commande 2 : $cmd \n");

if (! -e "${input_bamsorted}.filter2")
{
print STDERR "Echec lors de la generation du second fichier de filtre des alignements. \n";
}
else
{
`(cp -a \"${input_bamsorted}.filter2\" $output_filter2) >& ./cp_filtrer_aln_fichier2.log 2>&1`
}
