#!/usr/bin/perl -w

# usage : perl sm_gatk_3.pl <BAM file> <BAM intervals> <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 $input_fasta           = $ARGV[0];
my $chromosome_name       = $ARGV[1];
my $input_bam_first_GATK  = $ARGV[2];
my $input_bam_intervals   = $ARGV[3];
my $output_realign_bam    = $ARGV[4];

#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 $T ="IndelRealigner";
my $DICT ='/usr/local/bioinfo/src/galaxy-dev/galaxy-dist/tool-data/shared/jars/CreateSequenceDictionary.jar';


####################################################################
# 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($input_bam_intervals); #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
}


####################################################################
# ETAPE 6 : 2nde étape du réalignement
#
# Input :   fichier FASTA
#           fichier BAM
#           fichier bam.intervals
# Output :  fichier BAM
#           fichier .bai (index)
####################################################################
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T IndelRealigner -I SRR007327_sorted_
#rg_rmdup.bam -R NC_012125.1.fasta -targetIntervals bam.intervals -o SRR007327_sorted_rg_rmdup_realign.bam
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T IndelRealigner -I ERR_merge_sorted_
#rg_rmdup.bam -R NC_012125.1.fasta -targetIntervals bam.intervals -o ERR_merge_sorted_rg_rmdup_realign.bam

#extension .bam indispensables pour GATK
$cmd = "cp $input_bam_first_GATK \"${input_bam_first_GATK}.bam\";";
system $cmd;
$cmd = "cp $input_bam_intervals \"${input_bam_intervals}.bam.intervals\";";
system $cmd;
#le bam entrant doit être indexe
$cmd = "($SAMTOOLS index \"${input_bam_first_GATK}.bam\") >& ./gatkindexbam.log 2>&1";
system $cmd; 
#Verifier l indexation des fichiers BAM


system ("(java -jar -Xmx4G $GATK -T $T -I ${input_bam_first_GATK}.bam -R $fasta -targetIntervals ${input_bam_intervals}.bam.intervals -o ${input_bam_first_GATK}_realign.bam) >& ./gatk_realign2.log 2>&1");
#pour les biologistes
print STDOUT "Realignment : \n\n java -jar -Xmx4G $GATK -T $T -I ${input_bam_first_GATK}.bam -R $fasta -targetIntervals ${input_bam_intervals}.bam.intervals -o ${input_bam_first_GATK}_realign.bam \n\n";

#Tests
#perl sm_gatk_3.pl tests/TP/NC_012125.1.fasta tests/TP/ERR_merge_sorted_rg_rmdup.bam 23-4_17h17mn32sec.bam.intervals tests/output1


#Elimination des fichiers inutiles apres l execution de la ligne de commandes
$cmd = "rm \"${input_bam_first_GATK}.bam\"; rm \"${input_bam_intervals}.bam.intervals\";";
system $cmd;

#Recuperation des fichiers dans Galaxy
if (! -e "${input_bam_first_GATK}_realign.bam"){print STDERR "Realign bam FILE NOT FOUND\n";}
else{`cp -a "${input_bam_first_GATK}_realign.bam" $output_realign_bam ;`}

