#!/usr/bin/perl -w

# usage : perl sm_gatk_4.pl <BAM file> <VCF file> <FASTA file>  
#sarah maman - 24/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_realign     = $ARGV[2];
my $input_vcf		  = $ARGV[3];
my $vcf                   = $ARGV[4];
my $output_recal_bam      = $ARGV[5];

#pour eviter que la commande ne fonctionne plus si le $PATH est perdu suite a un restart de galaxy
my $cmd = '';
my $GATK = '/usr/local/bioinfo/src/galaxy-dev/galaxy-dist/tool-data/shared/jars/GenomeAnalysisTK.jar';
my $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my $BR ="BaseRecalibrator";
my $PR = "PrintReads";
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
my $ALL = $mday."-".($mon+1)."_".$hour."h".$min."mn".$sec."sec";
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 $path = dirname($input_bam_realign); #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 $path/$chromosome_name.fasta`;
my $fasta ="$path/$chromosome_name.fasta";


my $dict = "/galaxydata/database/files/workspace/$chromosome_name.dict";
my $dicttmp = "$path/$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 $path/.; chmod -R 777 $dicttmp; chmod -R 777 $dict`; #Copie (seulement, pas cmde java) a faire dans chaque module GATK
}


####################################################################
# Empty or not empty VCF ?
####################################################################

open(FILE,">$path/${ALL}_empty.vcf") or die "Unable to open file $path/${ALL}_empty.vcf";
if ($vcf eq "NO"){
print FILE "##fileformat=VCFv4.1\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tG1";
}
else
{
`mv $input_vcf "$path/${ALL}_empty.vcf"`;
}


####################################################################
# ETAPE 7 : 1ère étape de la recalibration
#
# Input :   fichier FASTA
#           fichier BAM
#           fichier VCF
# Output :  fichier texte .grp
####################################################################
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T BaseRecalibrator -I SRR007327_sorte
#d_rg_rmdup_realign.bam -R NC_012125.1.fasta -o SRR.grp -knownSites empty.vcf
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T BaseRecalibrator -I ERR_merge_sorte
#d_rg_rmdup_realign.bam -R NC_012125.1.fasta -o ERR.grp -knownSites empty.vcf

#GATK a besoin de l extension pour ses traitements
`cp $input_bam_realign "${input_bam_realign}.bam";`;
#le bam entrant doit être indexe
`($SAMTOOLS index "${input_bam_realign}.bam") >& ./gatkindexbam.log 2>&1`; #Verifier l indexation des fichiers BAM


system ("(java -jar -Xmx4G $GATK -T $BR -I ${input_bam_realign}.bam -R $fasta -o $path/${ALL}_SRR.grp -knownSites $path/${ALL}_empty.vcf) >& ./gatk_recal1.log 2>&1");

#Pour les biologistes 
print STDOUT "First step of recalibration :\n\n java -jar -Xmx4G $GATK -T $BR -I ${input_bam_realign}.bam -R $fasta -o $path/${ALL}_SRR.grp -knownSites $path/${ALL}_empty.vcf \n\n";

####################################################################
# ETAPE 8 : 2nde étape de la recalibration
#
# Input :   fichier FASTA
#           fichier BAM
#           fichier .grp
# Output :  fichier BAM
#           fichier .bai (index)
####################################################################
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T PrintReads -I SRR007327_sorted_rg_r
#mdup_realign.bam -R NC_012125.1.fasta -BQSR SRR.grp -o SRR007327_sorted_rg_rmdup_realign_recal.bam
#java -jar -Xmx4G /usr/local/bioinfo/src/GATK/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T PrintReads -I ERR_merge_sorted_rg_r
#mdup_realign.bam -R NC_012125.1.fasta -BQSR ERR.grp -o ERR_merge_sorted_rg_rmdup_realign_recal.bam

system ("(java -jar -Xmx4G $GATK -T $PR -I ${input_bam_realign}.bam -R $fasta -BQSR $path/${ALL}_SRR.grp -o ${input_bam_realign}_recal.bam) >& ./gatk_recal2.log 2>&1");

#Pour les biologistes 
print STDOUT "2nd step of recalibration :\n\n java -jar -Xmx4G $GATK -T $PR -I ${input_bam_realign}.bam -R $fasta -BQSR $path/${ALL}_SRR.grp -o ${input_bam_realign}_recal.bam \n\n";

#Tests
#perl sm_gatk_4.pl tests/TP/NC_012125.1.fasta tests/TP/ERR_merge_sorted_rg_rmdup_realign.bam tests/TP/ERR_merge_sorted.bam_empty.vcf tests/TP/output1bam tests/TP/output1bai

#Elimination des fichiers inutiles apres traitements
`rm "${input_bam_realign}.bam";`;



#Recuperation des fichiers dans Galaxy
if (! -e "${input_bam_realign}_recal.bam"){print STDERR "Recal bam FILE NOT FOUND\n";}
else{`cp -a "${input_bam_realign}_recal.bam" $output_recal_bam ;`}
