gatk4_gather_sort_vcfs.pbs

#!/bin/bash

#########################################################
# 
# Platform: NCI Gadi HPC 
# Usage: 
# Version: 2.0 
# 
# For more details see: https://github.com/Sydney-Informatics-Hub/Germline-ShortV 
# 
# If you use this script towards a publication, support us by citing: 
# 
# Suggest citation: 
# Sydney Informatics Hub, Core Research Facilities, University of Sydney, 
# 2021, The Sydney Informatics Hub Bioinformatics Repository, <date accessed>, 
# https://github.com/Sydney-Informatics-Hub/Germline-ShortV 
# 
# Please acknowledge the Sydney Informatics Hub and the facilities: 
# 
# Suggested acknowledgement: 
# The authors acknowledge the technical assistance provided by the Sydney 
# Informatics Hub, a Core Research Facility of the University of Sydney 
# and the Australian BioCommons which is enabled by NCRIS via Bioplatforms 
# Australia. The authors acknowledge the use of the National Computational 
# Infrastructure (NCI) supported by the Australian Government. 
# 
#########################################################

#PBS -P <project>
#PBS -N gather_sort
#PBS -l walltime=00:00:00,ncpus=,mem=GB,wd
#PBS -q express
#PBS -W umask=022
#PBS -l storage=
#PBS -l jobfs=8GB
#PBS -o ./Logs/gatk4_gather_sort.o
#PBS -e ./Logs/gatk4_gather_sort.e

set -e

module load gatk/4.1.8.1

config=/path/to/cohort.config
cohort=$(basename $config | cut -d'.' -f1)
INPUTS=./Inputs
vcfdir=../${cohort}_GenotypeGVCFs
scatterdir=../Reference/ShortV_intervals
scatterlist=$(ls $scatterdir/*.list)
if [[ ${#scatterlist[@]} > 1 ]]; then
        echo "$(date): ERROR - more than one scatter list file found: ${scatterlist[@]}"
        exit
fi
logdir=./Logs/GATK4_Gather_Sort_Vcfs
errdir=./Logs/GATK4_Gather_Sort_Vcfs_error_capture
log=${logdir}/${cohort}.log
args=${INPUTS}/gatk4_gather_sort_vcfs.inputs
unsorted=${vcfdir}/${cohort}.vcf.gz
sorted=${vcfdir}/${cohort}.sorted.vcf.gz
num_int=`wc -l ${scatterlist} | cut -d' ' -f 1`

mkdir -p ${logdir} ${errdir}
rm -rf ${args} ${logfile}

for interval in $(seq -f "%04g" 0 $((${num_int}-1))); do
	echo "-I" ${vcfdir}/${cohort}.${interval}.vcf.gz >> ${args}
done

# Gather per interval cohort vcf.gz files into a single vcf.gz file
# Faster to gather .g.
# Less KSU to create MD5 later
echo "$(date) : GATK4 GatherVcfs. Gather $cohort.interval.vcf.gz files into $unsorted"
gatk GatherVcfs --java-options "-XX:ParallelGCThreads=${NCPUS} -Djava.io.tmpdir=${PBS_JOBFS}" \
	--arguments_file ${args} \
	--MAX_RECORDS_IN_RAM 100000000 \
	-O ${unsorted} >> ${log} 2>&1


# Sort g.vcf.gz file - in case GatherVcfs fails to make index file (will happen if you output gzipped file)
# Ensures g.vcf.gz is sorted
echo "$(date) : GATK 4 SortVcf. Sort $unsorted and output to $sorted"
gatk SortVcf --java-options "-XX:ParallelGCThreads=${NCPUS} -Djava.io.tmpdir=${PBS_JOBFS}"\
	-I ${unsorted} \
	-O ${sorted} >> ${log} 2>&1

# Check logs for GATK errors
if grep -q -i error $log
then
        printf "Error in GATK log ${log}\n" >> $err
fi

if grep -q Exception $log
then
        printf "Exception in GATK log ${log}\n" >> $err
fi