call_var_mutect2.sh

#!/user/bin/env bash
set -o errexit
set -o pipefail
set -o nounset

# default arg
check='T'
callgermline='true'

# help message
help_message="
Wrapper to run Mutect2 (GATK v4.2.2).

usage:
    bash $(basename "$0") [-options] -t <TUMOR_BAM> -f <FASTA>
required arguments:
    -t|--tumor : tumor sample BAM file [BAM]
    -f|--fasta : whole genome FASTA file [FASTA]
optional arguments:
    -n|--normal : matched normal sample [BAM] (default = NULL)
    --intervals : intervals file [BED,interval] (default = NULL)
    --chr : chromosome to restrict analysis to [string] (default = NULL)
    --germline : germline AF resource [VCF] (default = NULL)
    --pon : PON file [VCF] (default = NULL)
    --callgermline : whether to call germline sites [true|false] (default = true)
    -o|--outdir : output directory [string] (default = pwd)
    -l|--logdir : output directory for log files [string] (default = --outdir)
    --name : prefix to give output files [string] (default = extracted from bam)
    --check : whether to check input files [T|F] (default=T)
additional info:
    # check argument is useful for scheduling jobs where input files 
      may not presently exist
    # control, intervals, dbsnp, cosmic and pon arguments are all optional but
      if provided, they will be passed to mutect2
output:
    # creates logs and variants directories in outdir containing 
      run log and output VCF files respectively
"

# parse arg
while [[ $# -gt 1 ]]; do
    key=$1
    case $key in
        -t|--tumor)
            tumor=$2
            shift
            ;;
        -f|--fasta)
            fasta=$2
            shift
            ;;
        -n|--normal)
            normal=$2
            shift
            ;;
        -o|--outdir)
            outdir=$2
            shift
            ;;
        -l|--logdir)
            logdir=$2
            shift
            ;;
        --name)
            name=$2
            shift
            ;;
        --check)
            check=$2
            shift
            ;;
        --intervals)
            intervals=$2
            shift
            ;;
        --chr)
            chr=$2
            shift
            ;;
        --germline)
            germline=$2
            shift
            ;;
        --pon)
            pon=$2
            shift
            ;;
        --callgermline)
            callgermline=$2
            shift
            ;;
        *)
            printf "\nERROR: Undefined argument provided: %s %s\n" $1 $2
            echo "$help_message"; exit 2
            ;;
    esac
    shift
done

# set logdir
if [[ -z ${outdir:-} ]]; then
    outdir=$PWD
fi
if [[ -z ${logdir:-} ]]; then
    logdir=$outdir
fi

# check required arg
if [[ -z ${tumor:-} ]]; then
	printf "\nERROR: no tumor argument provided\n"
	echo "$help_message"; exit 2
elif [[ -z ${fasta:-} ]]; then
    printf "\nERROR: no FASTA argument provided\n"
    echo "$help_message"; exit 2
fi

# check files if flagged
if [[ $check = T ]]; then
    if [[ ! -r ${tumor} ]]; then
        printf "\nERROR: tumor file is not readable: %s\n" $tumor
        echo "$help_message"; exit 2
    elif [[ ! -r ${fasta} ]]; then
        printf "\nERROR: FASTA file is not readable: %s\n" $fasta
        echo "$help_message"; exit 2
    fi
fi

# get sample name if not provided
if [[ -z ${name:-} ]]; then
	name=$(basename ${tumor})
    name=${name%%.*}
fi

# setup output directories
mkdir -p ${outdir}
mkdir -p ${logdir}

# set optional arguments
if [[ ! -z ${normal:-} ]]; then
   if [[ $check = T ]] && [[ ! -r ${normal} ]]; then
        printf "\nERROR: normal file is not readable: %s\n" ${normal}
        echo "$help_message"; exit 2
    else
        normal_input="-I ${normal}"
        module load SAMtools/1.11-GCCcore-10.2.0
        normal_sm=$(samtools view -H ${normal} | grep "@RG" | grep -Po 'SM:\K[^\s]+')
        normal_name="-normal ${normal_sm}"
    fi
fi
if [[ ! -z ${intervals:-} ]]; then
   if [[ $check = T ]] && [[ ! -r ${intervals} ]]; then
        printf "\nERROR: intervals file is not readable: %s\n" ${intervals}
        echo "$help_message"; exit 2
    else
        intervals_arg="-L ${intervals}"
    fi
fi
if [[ ! -z ${germline:-} ]]; then
    if [[ $check = T ]] && [[ ! -r ${germline} ]]; then
        printf "\nERROR: germline file is not readable: %s\n" ${germline}
        echo "$help_message"; exit 2
    else
        germline_arg="--germline-resource ${germline}"
    fi
fi
if [[ ! -z ${pon:-} ]]; then
    if [[ $check = T ]] && [[ ! -r $pon ]]; then
        printf "\nERROR: PON file is not readable: %s\n" ${pon}
        echo "$help_message"; exit 2
    else
        pon_arg="--panel-of-normals ${pon}"
    fi
fi
if [[ ! -z ${chr:-} ]]; then
    chr_arg="-L $chr"
fi
if [[ ${callgermline} = "true" ]]; then
    callgermarg="--genotype-germline-sites true"
fi

# call variants
std_log="${logdir}/${name}.call_var_mutect2.log"
scr_log="${logdir}/${name}.call_var_mutect2.scr"
script=$(cat <<- EOS
		#!/bin/bash
		#SBATCH --time=24:00:00
		#SBATCH -n 8
		#SBATCH -N 1
		#SBATCH --mem=32gb
		#SBATCH --job-name=mutect2
		#SBATCH --output=$std_log
		${depend:-}
        
		# load modules
		module load Java/11.0.16
		module load SAMtools/1.11-GCCcore-10.2.0
		module load picard/2.25.6-Java-11
		module load GATK/4.2.2.0-GCCcore-10.2.0-Java-11

		# call
		gatk Mutect2 \
			${normal_input:-} \
			${normal_name:-} \
			${intervals_arg:-} \
			${chr_arg:-} \
			${germline_arg:-} \
			${pon_arg:-} \
			${callgermarg:-} \
			--bam-output ${name}.mutect2.bam \
			-R ${fasta} \
			-I ${tumor} \
			-O ${name}.mutect2.vcf.gz
	EOS
)
echo "$script" > ${scr_log}
subm=$(sbatch "$scr_log")
echo $subm
exit 0