#!/usr/bin/perl -w

# usage : perl sm_gatk_5.pl <BAM file> <BAM file> <FASTA file>  
# Sarah Maman - 23/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 $output_variants_vcf   = $ARGV[0];
my $output_idx            = $ARGV[1];
my $input_fasta           = $ARGV[2];
my $chromosome_name       = $ARGV[3];
my $glm                   = $ARGV[4];
my $inputs_recal_bam      = @ARGV[5..$#ARGV];
	

#pour eviter que la commande ne fonctionne plus si le $PATH est perdu suite a un restart de galaxy
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 $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my ($nb) = ($output_variants_vcf=~/galaxy_dataset_(\d+)\.\S+$/);
my $UG ="UnifiedGenotyper";


#pour recuperer chaque bam en argument -I

	my $config = "config_".$nb.".txt";
	open FILE,">$config";
	my $cpt=1;
        foreach my $file (@ARGV[5..$#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 : 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_recal_bam); #bien nommer le BAM en cours d'utilisation
#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
}


####################################################################
# ETAPE 9 : Recherche de variants
#
# Input :   fichier FASTA
#           fichiers BAM
# Output :  fichier VCF
#           fichier .idx (index du VCF)
####################################################################
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T UnifiedGenotyper -I ERR_merge_sorte
#d_rg_rmdup_realign_recal.bam -I SRR007327_sorted_rg_rmdup_realign_recal.bam -R NC_012125.1.fasta -o variants.vcf -glm BOTH
system ("(java -jar -Xmx4G $GATK -T $UG $code -R $fasta -o $workdir/$nb.vcf -glm $glm) >& ./gatk-variants.log 2>&1");

#Info pour les biologistes
print STDOUT "Search variants : \n\n java -jar -Xmx4G $GATK -T $UG $code -R $fasta -o $workdir/$nb.vcf -glm $glm \n\n";

close FILE;

#Tests
#perl sm_gatk_5.pl tests/TP/toto tests/TP/toto2 tests/TP/NC_012125.1.fasta tests/TP/SRR007327_sorted_rg_rmdup_realign_recal.bam tests/TP/ERR_merge_sorted_rg_rmdup_realign_recal.bam


#Recuperation des fichiers dans Galaxy
if (! -e "$workdir/$nb.vcf"){print STDERR "Variants VCF FILE NOT FOUND\n";}
else{`cp -a "$workdir/$nb.vcf" $output_variants_vcf ;`
#Affichage d'un additional Files as output :
#peak_calling/macs_wrapper.py:    out_html.write( '<html><head><title>Additional output created by MACS (%s)</title></head><body><h3>Additional Files:</h3><p><ul>\n' % experiment_name )
#print("<html><head><title>Additional output created by Variant Search  $output_idx </title></head><body><h3>Additional Files:</h3><p><ul>\n");
}

if (! -e "$workdir/$nb.vcf.idx"){print STDERR "IDX FILE NOT FOUND\n";}
else{`cp -a "$workdir/$nb.vcf.idx" $output_idx ;`}
