gatk_snpcall.scr

#!/bin/sh
#$ -cwd
# -l mem=6G,time=8::

# This step in the pipeline performs SNP calling on the inout bam file.

GLOBAL="global_config.sh"

if [[ -e $GLOBAL ]]
then
        . $GLOBAL
else
        echo "Global config file not found. Exiting."
        exit 1
fi

USAGE="Usage: $0 -I <Input bam file> -R <Reference Fasta> -D <DBSNP file> [-L \"#:#-#\"]"

while getopts I:L:R:D:P:h o
do      case "$o" in
        I)      INP="$OPTARG";;
        L)      CHR="$OPTARG";;
        R)      REF="$OPTARG";;
        D)      DBSNP="$OPTARG";;
        P)      Platform="$OPTARG";;
        h)      echo $USAGE
                exit 1;;
        esac
done

if [[ $INP == "" || $REF == "" || $DBSNP == "" || $Platform == "" ]]
then
        echo $USAGE
        exit 1
fi

CHR_NAME=`$SAMTOOLS idxstats $INP | grep chr`

if [[ $CHR != "" ]]
then
        if [[ `grep chr $CHR` == "" && -n $CHR_NAME ]]
        then
                CHR="chr${CHR}"
        fi
fi

if [[ $Platform == "illumina" ]]
then
	Platform="UNKNOWN"
fi

$GATK \
 -L $CHR \
 -R $REF \
 -D $DBSNP \
 -T UnifiedGenotyper \
 -I $INP \
 -o $INP.snps.raw.vcf \
 -stand_call_conf 10.0 \
 -stand_emit_conf 10.0 \
 -pl $Platform \
 -dcov 50

#-bm EMPIRICAL \
#-pl $Platform
#-A DepthOfCoverage \
#set call conf based on whether it's deep or low pass data