Calling variants GATK4 - HipSci

# index the genome
$ bwa index hs37d5.fa

# create a genome dictionary file
$ picard-tools CreateSequenceDictionary R=genome.fa O=genome.dict

# align reads using 'bwa mem'
$ls -1 *.gz | cut -d. -f1 | sort | uniq | while read l; do bwa mem  -t 40 ~/00-NGS/Annotation/HipSci/hs37d5.fa $l.1.val.1.fq.gz $l.2.val.2.fq.gz > $l.PE.sam ; done &

# sam to bam
$ for s in *.bam; do samtools view -b $s -o $s.bam ; done

# Picard 'AddOrReplaceReadGroups'
$ for b in *.bam; do picard-tools AddOrReplaceReadGroups I=$b O=$b.rg.sorted.bam SO=coordinate RGID=$b RGLB=TruSeq_DNA_prep RGPL=illumina RGPU=hiseq2000 RGSM=$b & done

# Picard 'MarkDuplicates'
$ for b in *.rg.sorted.bam; do picard-tools MarkDuplicates I=$b O=$b.dedupped.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=$b.output.metrics & done
 
# GATK 'SplitNCigarReads'
$ for d in *.dedupped.bam; do gatk SplitNCigarReads -R ~/00-NGS/Annotation/HipSci/hs37d5.fa -I $d -O $d.split.bam & done

# GATK 'BaseRecalibrator' & 'ApplyBQSR' (base quality score recalibration)
$ for s in *.split.bam; do gatk BaseRecalibrator -I $s -R ~/00-NGS/Annotation/HipSci/hs37d5.fa --known-sites ~/00-NGS/Annotation/HipSci/ALL.wgs.dbsnp.build135.snps.sites.vcf.gz -O $s.recal_data.table & done
$ for s in *.split.bam; do gatk ApplyBQSR -R ~/00-NGS/Annotation/HipSci/hs37d5.fa -I $s --bqsr-recal-file $s.recal_data.table -O $s.recal.bam & done

# GATK 'HaplotypeCaller'
$ for r in *.recal.bam; do gatk --java-options "-Xmx8g" HaplotypeCaller -R ~/00-NGS/Annotation/HipSci/hs37d5.fa -I $r -O $r.raw.vcf.gz & done