Lecture 2

Set up directory and copy over files

Remember to load modules

---

Index the bam file

---

Run samtools flagstat on both files

---

What are some of the similarities and differences?

Download the bam file to your local computer as previously shown using scp

---

Check what assembly it was aligned to

samtools view -H HG002.sorted.bam | tail -n 3

---

Open it in IGV with the appropriate reference around 20:43200000-43300000

Assign paths to variables for simplicity

REF=/home/hgen_share/bwa
VCF=/home/hgen_share/vcf/b37
PICARD_JAR=/home/hgen_share/utils/picard.jar
GATK_JAR=/home/hgen_share/utils/GenomeAnalysisTK-3.8.1.jar
export PATH="/home/hgen_share/utils/gatk-4.1.9.0:$PATH"
module load StdEnv/2020 samtools/1.16
module load nixpkgs/16.09 java/1.8.0_192

Find regions with indels and dense in SNPs with RealignerTargetCreator

java -Xmx4G -jar $GATK_JAR \
  -T RealignerTargetCreator \
  -R $REF/b37.fa \
  -L 20:43200000-43300000 \
  --known $VCF/dbSNP.vcf \
  --known $VCF/Mills.indels.vcf \
  -o realign.intervals \
  -I HG002.sorted.bam

---

Re-align around the identified regions with IndelRealigner

java -Xmx4G -jar $GATK_JAR \
  -T IndelRealigner \
  -R $REF/b37.fa \
  -known $VCF/dbSNP.vcf \
  -known $VCF/Mills.indels.vcf \
  -targetIntervals realign.intervals \
  -o HG002.realigned.sorted.bam \
  -I HG002.sorted.bam

Retrieve re-aligned reads based with "OC" tag ("original CIGAR", implying that a new CIGAR string was generated)

samtools view HG002.realigned.sorted.bam | grep 'OC' | cut -f1 | sort > oc.txt

---

Sort reads by name using picard SortSam to allow efficient searching

java -Xmx4G -jar $PICARD_JAR SortSam \
  -I HG002.realigned.sorted.bam \
  -O HG002.realigned.qs.bam \
  -SO queryname
 
java -Xmx2g -jar $PICARD_JAR SortSam \
  -I HG002.sorted.bam \
  -O HG002.qs.bam \
  -SO queryname

---

Subset reads with name matching those found in oc.txt with picard FilterSamReads

java -Xmx4G -jar $PICARD_JAR FilterSamReads \
  -I HG002.realigned.qs.bam \
  -O HG002.new.bam \
  --FILTER includeReadList \
  -RLF oc.txt \
  -SO coordinate \
  -CREATE_INDEX true
 
java -Xmx2g -jar $PICARD_JAR FilterSamReads \
  -I HG002.qs.bam \
  -O HG002.old.bam \
  --FILTER includeReadList \
  -RLF oc.txt \
  -SO coordinate \
  -CREATE_INDEX true

---

Download the new and old bam files. What else must you download with your bam file for visualization?

---

Go to 20:43,257,030-43,257,326 to see an obvious example.

Mark PCR duplicate reads with MarkDuplicates

java -Xmx4G -jar ${PICARD_JAR} MarkDuplicates \
  -REMOVE_DUPLICATES false \
  -CREATE_INDEX true \
  -I HG002.realigned.sorted.bam \
  -O HG002.sorted.dup.bam \
  -M HG002.sorted.dup.metrics

---

Check the results

grep 'METRICS CLASS' -A 2 HG002.sorted.dup.metrics | column -t | less -S

Establish a set of covariates for correction based on context and known SNPs with BaseRecalibrator

gatk --java-options "-Xmx4G" \
  BaseRecalibrator \
  -R $REF/b37.fa \
  -known-sites $VCF/dbSNP.vcf \
  -known-sites $VCF/Mills.indels.vcf \
  -L 20:43200000-43300000 \
  -O HG002.sorted.dup.recalibration_report.grp \
  -I HG002.sorted.dup.bam

---

Recalibrate using the metrics calculated above with PrintReads

gatk --java-options "-Xmx4G" \
  ApplyBQSR \
  -R $REF/b37.fa \
  --bqsr-recal-file HG002.sorted.dup.recalibration_report.grp \
  -O HG002.sorted.dup.recal.bam \
  -I HG002.sorted.dup.bam
   

Perform variant calling with HaplotypeCaller

gatk --java-options "-Xmx4G" \
  HaplotypeCaller \
  -R $REF/b37.fa \
  -I HG002.sorted.dup.recal.bam \
  -ERC GVCF \
  -O HG002.g.vcf \
  -L 20:43200000-43300000 \
  --native-pair-hmm-threads 1

---

Download the vcf and recalibrated bam file to your local computer

Guide to viewing vcf files in IGV found here

More information, including parameter details, for GATK HaplotypeCaller found here