Mutect2_callPONs.sh

#!/bin/sh
#SBATCH -J NormMutect2
#SBATCH -o /fast/users/a1092098/launch/slurm-%j.out
#SBATCH -A robinson
#SBATCH -p batch
#SBATCH -N 1
#SBATCH -n 2
#SBATCH --time=48:00:00
#SBATCH --mem=8GB

# Notification configuration
#SBATCH --mail-type=END
#SBATCH --mail-type=FAIL
#SBATCH --mail-user=clare.vaneyk@adelaide.edu.au

# load modules
module load Java/1.8.0_121
module load GATK
module load SAMtools

# run the executable
# A script to generate PONs with gatk Mutect2, designed for the Phoenix supercomputer
# Variables that usually don't need changing once set for your system

usage()
{
echo "# A script to generate PONs with gatk Mutect2, designed for the Phoenix supercomputer
# Requires: GATK and a list of samples
#
# Usage sbatch --array 0-(nSamples-1) $0  -b /path/to/bam/files [-o /path/to/output] | [ - h | --help ]
#
# Options
# -S    REQUIRED. List of sample ID in a text file
# -b    REQUIRED. /path/to/bamfiles.bam. Path to where you want to find your bam files. Every file matching a sample ID will be used.
# -O    OPTIONAL. Path to where you want to find your file output (if not specified current directory is used)
# -h or --help  Prints this message.  Or if you got one of the options above wrong you'll be reading this too!
#
#
# Original: Derived from GATK.HC.Phoenix by Mark Corbett, 16/11/2017
# Modified: (Date; Name; Description)
# 21/06/2018; Mark Corbett; Modify for Haloplex
# 06/07/2018; Clare van Eyk; remove configuration and fix reference and dbSNP
# 09/07/2018; Clare van Eyk; modify for use with Mutect2 command from GATK

"
}

## Set Variables ##
while [ "$1" != "" ]; do
        case $1 in
                -S )                    shift
                                        SAMPLE=$1
                                        ;;
                -b )                    shift
                                        bamDir=$1
                                        ;;
                -O )                    shift
                                        VcfFolder=$1
                                        ;;
                -h | --help )           usage
                                        exit 0
                                        ;;
                * )                     usage
                                        exit 1
        esac
        shift
done

if [ -z "$SAMPLE" ]; then # If no SAMPLE name specified then do not proceed
        usage
        echo "#ERROR: You need to specify a list of sample ID in a text file eg.                                                                                                           /path/to/SAMPLEID.txt"
        exit 1
fi

# Define batch jobs based on samples
sampleID=($(cat $SAMPLE))

if [ -z "$bamDir" ]; then # If no bamDir name specified then do not proceed
        usage
        echo "#ERROR: You need to tell me where to find the bam files."
        exit 1
fi
if [ -z "$VcfFolder" ]; then # If no output directory then use default directory
        VcfFolder=$FASTDIR/WGS/Mosaic
        echo "Using $FASTDIR/WGS/Mosaic as the output directory"
fi
if [ ! -d $VcfFolder ]; then
        mkdir -p $VcfFolder
fi

tmpDir=$FASTDIR/tmp/${sampleID[$SLURM_ARRAY_TASK_ID]} # Use a tmp directory for all of the GATK and samtools temp files
if [ ! -d $tmpDir ]; then
        mkdir -p $tmpDir
fi

## Start of the script ##
##To create panel of normals, run:##
###On each control sample###
cd $tmpDir
gatk Mutect2 \
-I $bamDir/${sampleID[$SLURM_ARRAY_TASK_ID]}.dedup.realigned.recalibrated.bam \
-tumor ${sampleID[$SLURM_ARRAY_TASK_ID]} \
-R /data/biohub/Refs/human/gatk_bundle/2.8/b37/human_g1k_v37_decoy.fasta \
--germline-resource /data/neurogenetics/RefSeq/GATK/b37/af-only-gnomad.hg38.vcf.gz \
-O $VcfFolder/${sampleID[$SLURM_ARRAY_TASK_ID]}.mosaic.vcf >> $tmpDir/${sampleID[$SLURM_ARRAY_TASK_ID]}.pipeline.log 2>&1