#!/bin/bash -x

#SBATCH --job-name="GATK_pipeline"
#SBATCH --ntasks 8
##SBATCH --nodes 1
#SBATCH --partition long
#SBATCH --mail-type=ALL
#SBATCH --mail-user=raonyguimaraes@gmail.com

#$SLURM_JOB_ID
##################### DEFINE VARIABLES ###############################

EMAIL=raonyguimaraes@gmail.com

#PROGRAMS
BWA_DIR=../../bin/bwa-0.5.9
GATK_DIR=../../bin/GenomeAnalysisTK-1.1-23-g8072bd9
PIC_DIR=../../bin/picard-tools-1.52
ST_DIR=../../bin/samtools-0.1.17

#FOLDERS
TMP_DIR=/var/tmp/bio712
OUT_DIR=output
FASTQ_DIR=../../input/exome
LOG_DIR=logs/gatk
REFERENCE=../../input/b37/human_g1k_v37.fasta
INPUT=alignment/samtools/exome.sorted.bam

EXON_CAPTURE_FILE=../../input/bed/Design_Annotation_files/Target_Regions/targets.merged.bed

#Variant Score Recalibrator variables
#DBSNP=../../input/b37_resources/dbsnp_132.b37.vcf
DBSNP=../../input/dbsnp/dbsnp-134.vcf
OMNI=../../input/b37_resources/1000G_omni2.5.b37.sites.vcf
HAPMAP=../../input/b37_resources/hapmap_3.3.b37.sites.vcf
INDELS=../../input/b37_resources/1000G_biallelic.indels.b37.vcf

#create local tmp dir to speedup the processing (do not use NFS)
rm -rf /var/tmp/bio712
mkdir /var/tmp/bio712


#################### GATK ############################################

# GATK  5 steps on this order:  
#Local realignment around indels, MarkDuplicates, Base quality score recalibration,  Unifier Genotyper, Variant quality score recalibration

#include Short Read Micro re-Aligner SRMA

##########################################STEP 1: Local realignment around indels
# 
# # 1.1 RealignerTargetCreator
#java  -Xmx4g -jar -Djava.io.tmpdir=$TMP_DIR -jar $GATK_DIR/GenomeAnalysisTK.jar -T RealignerTargetCreator \
#-l INFO \
#-R $REFERENCE \
#-I $INPUT \
#-B:dbsnp,vcf $DBSNP \
#-B:indels,vcf $INDELS \
#-B:intervals,BED $EXON_CAPTURE_FILE \
#-log $LOG_DIR/RealignerTargetCreator.log \
#-o $OUT_DIR/exome.intervals

#-L $EXON_CAPTURE_FILE \


# 1.2 IndelRealigner 
#java -Xmx22g -Djava.io.tmpdir=$TMP_DIR -jar $GATK_DIR/GenomeAnalysisTK.jar -T IndelRealigner \
#-R $REFERENCE \
#-I $INPUT \
#-targetIntervals $OUT_DIR/exome.intervals \
#-log $LOG_DIR/IndelRealigner.log \
#-o $OUT_DIR/exome.real.bam


# ##########################################STEP 2 MarkDuplicates
# 
#MARKDUPLICATES
#java -Xmx4g -jar $PIC_DIR/MarkDuplicates.jar \
#INPUT=$OUT_DIR/exome.real.bam \
#REMOVE_DUPLICATES=true \
#VALIDATION_STRINGENCY=LENIENT \
#AS=true \
#METRICS_FILE=alignment/picard/exome.dedup.metrics \
#OUTPUT=$OUT_DIR/exome.real.dedup.bam
#AS = ASSUME_SORTED
# 
# #Reindex BAM FIle - the'realigned'duplicates'removed'bam
#$ST_DIR/samtools index $OUT_DIR/exome.real.dedup.bam $OUT_DIR/exome.real.dedup.bam.bai
# 
# #Ignore Next 2 Steps !!!!
# #2.1 AddorReplaceGroups
# # java -jar $PIC_DIR/AddOrReplaceReadGroups.jar \
# # I=$OUT_DIR/exome.real.bam \
# # O=$OUT_DIR/exome.real.dedup.bam \
# # SORT_ORDER=coordinate \
# # RGID=EXOME RGLB=EXOME RGPL=illumina RGPU=EXOME RGSM=EXOME CREATE_INDEX=True \
# # VALIDATION_STRINGENCY=LENIENT \
# # TMP_DIR=/var/tmp/bio712
# # 
# # #2.2 FixMate --- should be used after realigment ???
# # java -jar  $PIC_DIR/FixMateInformation.jar \
# # INPUT=$OUT_DIR/exome.real.dedup.bam \
# # OUTPUT=$OUT_DIR/exome.real.dedup.matefixed.bam \
# # SORT_ORDER=coordinate \
# # VALIDATION_STRINGENCY=LENIENT \
# # TMP_DIR=/var/tmp/bio712
# 
# 
# ##########################################STEP 3 Base quality score recalibration
# 
# # #3.1 CountCovariates Before
#java -Xmx12g -jar $GATK_DIR/GenomeAnalysisTK.jar -T CountCovariates \
#-l INFO \
#-I $OUT_DIR/exome.real.dedup.bam \
#-R $REFERENCE \
#-B:dbsnp,VCF $DBSNP \
#-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate \
#-recalFile $OUT_DIR/exome.recal_data.before.csv \
#-log $LOG_DIR/CountCovariates_before.log \
#-nt 4
# # 
# # #3.2 TableRecalibration (try -baq RECALCULATE)
#java -jar $GATK_DIR/GenomeAnalysisTK.jar -T TableRecalibration \
#-l INFO \
#-I $OUT_DIR/exome.real.dedup.bam \
#-R $REFERENCE \
#-recalFile $OUT_DIR/exome.recal_data.before.csv \
#-log $LOG_DIR/TableRecalibration.log \
#-o $OUT_DIR/exome.real.dedup.recal.bam
# 
# #Reindex Bam FILE
#$ST_DIR/samtools index $OUT_DIR/exome.real.dedup.recal.bam $OUT_DIR/exome.real.dedup.recal.bam.bai
# 
# # 3.2.1 CountCovariates After
#java -Xmx12g -jar $GATK_DIR/GenomeAnalysisTK.jar -T CountCovariates \
#-l INFO \
#-I $OUT_DIR/exome.real.dedup.recal.bam \
#-R $REFERENCE \
#-B:dbsnp,VCF $DBSNP \
#-cov ReadGroupCovariate -cov QualityScoreCovariate -cov CycleCovariate -cov DinucCovariate \
#-recalFile $OUT_DIR/exome.recal_data.after.csv \
#-log $LOG_DIR/CountCovariates_after.log \
#-nt 4
# 
# 
# ##########################################STEP 4 UnifierGenotyper
# # 
# # #Standard Raw VCF
java -Xmx15g -jar $GATK_DIR/GenomeAnalysisTK.jar -T UnifiedGenotyper \
-l INFO \
-I $OUT_DIR/exome.real.dedup.recal.bam \
-R $REFERENCE \
-B:intervals,BED $EXON_CAPTURE_FILE \
-B:dbsnp,VCF $DBSNP \
-glm BOTH \
-stand_call_conf 50.0 \
-stand_emit_conf 20.0 \
-dcov 300 \
-A AlleleBalance \
-A DepthOfCoverage \
-A FisherStrand \
-o $OUT_DIR/exome.raw.vcf \
-log $LOG_DIR/UnifiedGenotyper.log \
-nt 4
# #-B:targetIntervals,BED $EXON_CAPTURE_FILE \
# 
# ##########################################STEP 5 Variant quality score recalibration
# 
# #SelectVariants.Select SNPS from the unified genotyper raw VCF
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T SelectVariants \
-R $REFERENCE \
-B:variant,VCF $OUT_DIR/exome.raw.vcf \
-snps \
-log $LOG_DIR/SelectVariants.snps.log \
-o $OUT_DIR/exome.snps.raw.vcf
# 
# 
# # #create folder
mkdir $OUT_DIR/VariantRecalibrator
# #VariantRecalibrator
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T VariantRecalibrator \
-R $REFERENCE \
-B:input,VCF $OUT_DIR/exome.snps.raw.vcf \
-B:hapmap,VCF,known=false,training=true,truth=true,prior=15.0 $HAPMAP \
-B:omni,VCF,known=false,training=true,truth=false,prior=12.0 $OMNI \
-B:dbsnp,VCF,known=true,training=false,truth=false,prior=8.0 $DBSNP \
-an QD -an HaplotypeScore -an MQRankSum -an ReadPosRankSum -an FS -an MQ \
--maxGaussians 6 \
-mode SNP \
-recalFile $OUT_DIR/VariantRecalibrator/exome.recal \
-tranchesFile $OUT_DIR/VariantRecalibrator/exome.tranches \
-rscriptFile $OUT_DIR/VariantRecalibrator/exome.plots.R \
-log $LOG_DIR/VariantRecalibrator.log \
-nt 4
# 
#  
# #Apply Recalibrator
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T ApplyRecalibration \
-R $REFERENCE \
-B:input,VCF $OUT_DIR/exome.snps.raw.vcf \
-tranchesFile $OUT_DIR/VariantRecalibrator/exome.tranches \
-recalFile $OUT_DIR/VariantRecalibrator/exome.recal \
--ts_filter_level 99.0 \
-log $LOG_DIR/ApplyRecalibration.log \
-o $OUT_DIR/exome.recal.snps.vcf

# #SelectVariants.Select Indels from the Unified Genotyper raw VCF
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T SelectVariants \
-R $REFERENCE \
-B:variant,VCF $OUT_DIR/exome.raw.vcf \
-indels \
-log $LOG_DIR/SelectVariants.indels.log \
-o $OUT_DIR/exome.indels.raw.vcf
# 
# #Variant Filtration
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T VariantFiltration \
-R $REFERENCE \
-B:variant,VCF $OUT_DIR/exome.indels.raw.vcf \
--filterExpression "QD < 2.0 || ReadPosRankSum < -20.0 || FS > 200.0"
--filterName GATK_standard \
--missingValuesInExpressionsShouldEvaluateAsFailing \
-log $LOG_DIR/VariantFiltration.indels.log \
-o $OUT_DIR/exome.indels.filtered.vcf
# 
# #CombineVariants: Combine Recalibrated SNPs and Filtered  Indels
java -Xmx4g -jar $GATK_DIR/GenomeAnalysisTK.jar -T CombineVariants \
-R $REFERENCE \
-B:variant,VCF $OUT_DIR/exome.recal.snps.vcf \
-B:variant,VCF $OUT_DIR/exome.indels.filtered.vcf \
-log $LOG_DIR/CombineVariants.log \
-o $OUT_DIR/exome.snps.indels.vcf

#Include Deph Of Coverage!!

rm -rf /var/tmp/bio712