#!/usr/bin/perl -w

# usage : perl bwa_samtools.pl <FASTA ref file> <FASTQ1 file> <FASTQ1 file> <output BAM file> 
#sarah maman - 11 juillet 2013 - Mis en place dans le cadre du projet PhyloFish
#perl bwa_samtools.pl tests/results_from_merged_rename.fasta toto tests/TEST_R1.fastq tests/TEST_R2.fastq galaxy_dataset_128.dat

use strict;
use File::Basename;

my $input_fasta        = $ARGV[0];
my $tailleFichier      = (stat($input_fasta))[7]; #Selon la taille du FASTA de reference, l indexation est differente
my $output_stats       = $ARGV[1];
my $input_1_fastq      = $ARGV[2];
my $input_2_fastq      = $ARGV[3];
my $output_BAM         = $ARGV[4];

my $path = dirname($input_fasta);
my $nb = ($output_BAM=~/galaxy_dataset_(\d+)\.\S+$/);



#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 $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my $cmd1 = ''; my $cmd2 = ''; my $cmd3 = ''; my $cmd4 = ''; my $cmd5 = ''; my $cmd6 = ''; my $cmd7 = '';

#fonction exit
sub control ()
{
if ($? != 0){print STDERR "Incident. Fin du job. Veuillez relancer votre outil\n"; exit (1);}
}


####################################################################################
#									           #			
#		INDEXATION FASTA REF						   #
#                                                                                  #
####################################################################################
#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 (! -e "${input_fasta}.ann"){#si le fasta est deja indexe
	if ($tailleFichier > 10485760) {
        	$cmd1 = "$BWA index -a bwtsw $input_fasta  >> ./bwaindex.log 2>&1";
     	}
	if ($tailleFichier < 2147483648){
        	$cmd1 = "$BWA index -a is $input_fasta  >> ./bwaindex.log 2>&1";
     	}
	else {
        	$cmd1 ="$BWA index $input_fasta  >> ./bwaindex.log 2>&1";
     	}
	system $cmd1;

	#Information en STDOUT pour les biologistes
	print STDOUT "Etape 1 \nIndexation : $cmd1 \n\n";
	control();
}
else {print STDOUT "Etape 1 \nIndexation  deja realisee en amont.\n\n";}


####################################################################################
#									           #			
#		BWA - OUTPUT : BAM NON TRIE | RECUPERATION SEQ ALIGNEE		   #
#                                                                                  #
####################################################################################
#Alignement
#-1 	When -b is specified, only use the first read in a read pair in mapping (skip single-end reads and the second reads).
#-2 	When -b is specified, only use the second read in a read pair in mapping. 

$cmd2 ="($BWA aln $input_fasta $input_1_fastq > ${input_1_fastq}_1.sai) >& ./bwaaln1.log 2>&1";
system $cmd2;
if (! -e "${input_1_fastq}_1.sai"){print STDERR "SAI FILE1 NOT FOUND\n";}
control();

$cmd3 ="($BWA aln $input_fasta $input_2_fastq > ${input_2_fastq}_2.sai) >& ./bwaaln2.log 2>&1";
system $cmd3;
if (! -e "${input_2_fastq}_2.sai"){print STDERR "SAI FILE2 NOT FOUND\n";}
control();

#Information en STDOUT pour les biologistes
print STDOUT "Etape 2 \nAlignement du premier fastq : $cmd2 \nAlignement du second fastq: $cmd3 \n\n";

#bwa sampe / sequences alignees uniquement
#$cmd6 = "($BWA sampe $input_fasta ${input_1_fastq}-1.sai ${input_2_fastq}-2.sai $input_1_fastq $input_2_fastq | $SAMTOOLS view -bS -F 4 - | samtools sort - ${input_fasta}-sampe ) >& ./bwa.log 2>&1";
$cmd6 = "($BWA sampe $input_fasta ${input_1_fastq}_1.sai ${input_2_fastq}_2.sai $input_1_fastq $input_2_fastq | $SAMTOOLS view -bS -F 4 - | $SAMTOOLS sort - $path/${nb}_sampe) >& ./bwa.log 2>&1";
system $cmd6;
control();

#Information en STDOUT pour les biologistes
print STDOUT "Etape 3 \nBWA : $cmd4\nConversion des SAM en BAM : $cmd5\nBWA sur les sequences alignees uniquement : $cmd6 \n\n";


####################################################################################
#									           #			
#		TRI INDEX STATS      						   #
#                                                                                  #
####################################################################################
#Creation du BAM index et generation des statistiques 
#$cmd7 = "$SAMTOOLS index ${input_fasta}-sampe.bam ; $SAMTOOLS idxstats ${input_fasta}-sampe.bam > ${input_fasta}.txt;";
$cmd7 = "$SAMTOOLS index $path/${nb}_sampe.bam ; $SAMTOOLS idxstats $path/${nb}_sampe.bam > $path/$nb.txt;";
system $cmd7;
control();
#Information en STDOUT pour les biologistes
print STDOUT "Etape 4 \nCreation du BAM index et generation des statistiques : $cmd7\n\n";



####################################################################################
#									           #			
#		ELIMINATION DES FICHIERS INTERMEDIAIRES          		   #
#                                                                                  #
####################################################################################

#`rm ${input_2_fastq}_2.sai ${input_1_fastq}_1.sai $path/${nb}_sampe.bam.bai;`;


####################################################################################
#									           #			
#		GALAXY OUTPUT         						   #
#                                                                                  #
####################################################################################
# Recuperation du fichier sortant de statistiques
#if (! -e "${input_fasta}.txt"){print STDERR "Echec lors de la recuperation du fichier de statitiques. \n";}
#else {`cp -a "${input_fasta}.txt" $output_stats `;}
if (! -e "$path/$nb.txt"){print STDERR "Echec lors de la recuperation du fichier de statistiques. \n";}
else {`cp -a "$path/$nb.txt" $output_stats `;}

# Recuperation du fichier sortant BAM
#if (! -e "${input_fasta}-sampe.bam"){print STDERR "Echec lors de la recuperation du BAM. \n";}
#else {`cp -a "${input_fasta}-sampe.bam" $output_BAM `;}

if (! -e "$path/${nb}_sampe.bam"){print STDERR "Echec lors de la recuperation du BAM. \n";}
else {`cp -a "$path/${nb}_sampe.bam" $output_BAM `;}
