#!/usr/bin/perl -w

# usage : perl sm_gatk_2.pl <BAM file> <BAM file> <FASTA file>  
#sarah maman - 23/04/2013

use strict;
use File::Basename;


my $output_bam_intervals  = $ARGV[0];
my $input_fasta           = $ARGV[1];
my $chromosome_name       = $ARGV[2];
my $inputs_bam            = @ARGV[3..$#ARGV];


#pour eviter que la commande ne fonctionne plus si le $PATH est perdu suite a un restart de galaxy
my $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my $cmd = '';
my $GATK = '/usr/local/bioinfo/src/galaxy-dev/galaxy-dist/tool-data/shared/jars/GenomeAnalysisTK.jar';
my $DICT ='/usr/local/bioinfo/src/galaxy-dev/galaxy-dist/tool-data/shared/jars/CreateSequenceDictionary.jar';
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
my $ALL = $mday."-".($mon+1)."_".$hour."h".$min."mn".$sec."sec";



####################################################################
# ETAPE : Verifier la disponibilite du fichier dictionnaire du fasta 
#
# Input :   fichier FASTA
# Output :  fichier FASTA.dict
#
####################################################################
#Pour creer le dictionnaire :
#java -jar -Xmx4G /usr/local/bioinfo/src/picard-tools/current/CreateSequenceDictionary.jar R=NC_012125.1.fasta O=NC_012125.1.dict
#POur recuperer le dictionnaire deja cree :
#copier le dictionnaire de la racine vers le repertoire de travail
#my $workdir = dirname($inputs_bam); #bien nommer le BAM en cours d'utilisation
my $workdir = "/galaxydata/database/files/workspace/"; 
#input_fasta de galaxy est au format .dat donc il faut une copie au format fasta dans la work:
`cp $input_fasta $workdir/$chromosome_name.fasta`;
my $fasta ="$workdir/$chromosome_name.fasta";


my $dict = "/galaxydata/database/files/workspace/$chromosome_name.dict";
my $dicttmp = "$workdir/$chromosome_name.dict";
if (! -e "$dict"){
system ("(java -jar -Xmx4G $DICT R=$fasta O=$dicttmp) >& ./gatk_dict.log 2>&1");
#Apres la commande java de creation du dictionnaire
`cp -a $dicttmp $dict; chmod -R 777 $dicttmp; chmod -R 777 $dict`; #Copie (seulement, pas cmde java) a faire dans chaque module GATK
#Info pour les biologistes
print STDOUT "\nCreate FASTA dictionary \n\n java -jar -Xmx4G $DICT R=$fasta O=$dicttmp \n\n";
}
else{
`cp -a $dict $workdir/.; chmod -R 777 $dicttmp; chmod -R 777 $dict`; #Copie (seulement, pas cmde java) a faire dans chaque module GATK
}


####################################################################
#
#  Pour recuperer chaque bam en argument -I
#
####################################################################
my $config = "config_".$ALL.".txt";
	open FILE,">$config";
	my $cpt=1;
        foreach my $file (@ARGV[3..$#ARGV]){
                `cp $file "$file.bam";`;
		`($SAMTOOLS index "$file.bam") >& ./gatkindexbams.log 2>&1`; #Verifier l indexation des fichiers BAM
		print FILE "-I ","$file.bam"," "; 
		$cpt++;
	}
my $code = `cat $config`;
#print "\n\n******VOIR MON CONFIG : $code \n\n";

####################################################################
# ETAPE 5 : 1ère étape du réalignement
#
# Input :   fichier FASTA
#           fichiers BAM
# Output :  fichier texte bam.intervals
####################################################################
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T RealignerTargetCreator -I SRR007327
#_sorted_rg_rmdup.bam -I ERR_merge_sorted_rg_rmdup.bam -R NC_012125.1.fasta -o bam.intervals
my $T = "RealignerTargetCreator";

system ("(java -jar -Xmx4G $GATK -T $T  $code -R $fasta -o $workdir/$ALL.bam.intervals) >& ./gatk-realign1.log 2>&1");

#Info pour les biologistes
print STDOUT "Pre-realignement : \n\n java -jar -Xmx4G $GATK -T $T  $code -R $fasta -o $workdir/$ALL.bam.intervals\n\n ";


close FILE;

#Tests
#perl sm_gatk_2.pl /work/galaxy/database/files/006/25-4_13h41mn34sec.bam.intervals  /work/galaxy/database/files/006/NC_012125.1.fasta NC_012125.1 /work/galaxy/database/files/006/dataset_6989.dat.bam /work/galaxy/database/files/006/dataset_6985.dat.bam

#Recuperation des fichiers dans Galaxy
if (! -e "$workdir/$ALL.bam.intervals"){print STDERR "Intervals TXT FILE NOT FOUND\n";}
else{`mv "$workdir/$ALL.bam.intervals" $output_bam_intervals ;`}
