Lecture 3

Take a look at `HG002.g.vcf` from last week

View the file with less

Create a new directory and move the vcf file and its index to the new directory using mv

Filter variants using GATK recommanded cut-offs

Set up variables

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

Replace the hashtag (#) with the recommended cut-offs found here

gatk --java-options "-Xmx4G" \
    VariantFiltration \
    -V HG002.g.vcf \
    -filter "QD < #" --filter-name "QD" \
    -filter "FS > #" --filter-name "FS" \
    -filter "SOR > #" --filter-name "SOR" \
    -filter "MQ < #" --filter-name "MQ" \
    -filter "MQRankSum < #" --filter-name "MQRankSum" \
    -filter "ReadPosRankSum < #" --filter-name "ReadPosRankSum" \
    -O HG002.filt.g.vcf 2> >(grep -v undefined)

How many variants appear ok?

grep -v '^#' HG002.filt.g.vcf | cut -f7 | grep 'PASS' | wc -l

How many total variants are called? Hint: use grep

Set up variables for DeepVariant

module load StdEnv/2020 apptainer/1.1.3
mkdir -p "${PWD}/data"
mkdir -p "${PWD}/output"
INPUT_DIR="${PWD}/data"
OUTPUT_DIR="${PWD}/output"
TMPDIR=.
cp /home/hgen_share/lec3/* ${INPUT_DIR}

Run DeepVariant on the same sample from last week. Note: this may take a few minutes

singularity run -B $PWD,/usr/lib/locale/ \
   /home/hgen_share/deepvariant/deepvariant_1.4.0.sif \
   /opt/deepvariant/bin/run_deepvariant \
   --model_type=WGS \
   --ref="${INPUT_DIR}"/b37.chr20.fa \
   --reads="${INPUT_DIR}"/HG002.sorted.bam \
   --regions "20:43200000-43300000" \
   --output_vcf="${OUTPUT_DIR}"/HG002.dv.vcf.gz \
   --output_gvcf="${OUTPUT_DIR}"/HG002.dv.g.vcf.gz \
   --num_shards=1

Unzip the g.vcf found in the output directory using gunzip then filter it using similar parameters as above. Compare the number of variants that pass. Note: you may need to set-up your variables for GATK again

(Extra) Compare hard filtering with VQSR

Copy whole genome variant calls

cp /home/hgen_share/lec3_part2 ./

Since SNPs and indels need to be trained separately, start with SNPs

gatk --java-options "-Xmx4G" \
  VariantRecalibrator \
   -R ${REF}/b37.fa \
   -V raw.vcf \
   --resource:hapmap,known=false,training=true,truth=true,prior=15.0 $VCF/b37/hapmap.vcf \
   --resource:omni,known=false,training=true,truth=false,prior=12.0 $VCF/b37/omni.vcf \
   --resource:1000G,known=false,training=true,truth=false,prior=10.0 $VCF/b37/1KG.snps.vcf \
   --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $VCF/b37/dbSNP.vcf \
   -an DP -an QD -an FS -an SOR -an MQ -an MQRankSum -an ReadPosRankSum \
   -mode SNP \
   -O recal_SNP.recal \
   --tranches-file recal_SNP.tranches \
   --rscript-file recal_SNP.plots.R

Apply VQSR to SNPs

gatk --java-options "-Xmx4G" \
  ApplyVQSR \
  -R ${REF}/b37.fa \
  -V raw.vcf \
  -O raw.snp.vcf \
  --truth-sensitivity-filter-level 99.5 \
  --tranches-file recal_SNP.tranches \
  --recal-file recal_SNP.recal \
  -mode SNP

Now repeat for indels

gatk --java-options "-Xmx4G" \
  VariantRecalibrator \
  -R ${REF}/b37.fa \
  -V raw.snp.vcf \
  --resource:mills,known=false,training=true,truth=true,prior=12.0 $VCF/b37/Mills.indels.vcf \
  --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 $VCF/b37/dbSNP.vcf \
  --max-gaussians 4 \
  -an QD -an DP -an FS -an SOR -an MQRankSum -an ReadPosRankSum \
  -mode INDEL \
  -O recal_INDEL.recal \
  --tranches-file recal_INDEL.tranches \
  --rscript-file recal_INDEL.plots.R
 
gatk --java-options "-Xmx4G" \
  ApplyVQSR \
  -R ${REF}/b37.fa \
  -V raw.snp.vcf \
  -O recal.vcf \
  --truth-sensitivity-filter-level 99.0 \
  --tranches-file recal_INDEL.tranches \
  --recal-file recal_INDEL.recal \
  -mode INDEL

How many variants in the same region passed?

gatk SelectVariants \
  -R ${REF}/b37.fa \
  -V recal.vcf \
  -O recal.sub.vcf \
  -L 20:43200000-43300000
 
grep -v '^#' recal.sub.vcf | cut -f7 | grep 'PASS' | wc -l

(Extra) Generate a VCF file with only variants that passed

Set up bcftoolsto filter PASS variants

module load StdEnv/2020 gcc/9.3.0 bcftools/1.16

Run bcftools. The input file can be the VCF file produced from either variant caller

bcftools view -f 'PASS,.' -O v -o output_filename input_file

(Extra) Use `VariantsToTable`to get a tsv file of variants

VariantsToTable, by default, only extracts variants that passed filtering. Use the --show-filtered parameter to show all variants.

REF=/home/hgen_share/bwa
gatk VariantsToTable \
    -R $REF/b37.fa \
    -V input_vcf \
    -F CHROM -F POS -F FILTER -F TYPE -GF AD -GF DP \
    --show-filtered \
    -O output_filename.tsv

You may need to reload modulesfor GATK again

Add functional annotations with SnpEff

Copy whole genome variant calls

cp /home/hgen_share/lec3_part2/* ./

Use SnpEff with mostly default parameters

UTILS=/home/partage/utils
SNPEFF_JAR=/home/hgen_share/utils/snpEff/snpEff.jar
java -Xmx4G -jar $SNPEFF_JAR \
  -c /home/hgen_share/utils/snpEff/snpEff.config \
  -i vcf -o vcf -v hg19 \
  ex.vcf > ex.snpeff.vcf

Download .html results file to your local computer using scp

Annotate with existing resources using GEMINI

Add the directory containing gemini

export PATH=/home/hgen_share/utils/gemini/anaconda/bin:$PATH

Generate annotated database with gemini

module load mugqic/gemini/0.20.1
gemini load --cores 4 --tempdir $PWD -t snpEff -v ex.snpeff.vcf ex.snpeff.db

Dump database to plain-text csv using sqlite

db=ex.snpeff.db
t=($(sqlite3 $db ".tables"))
for i in "${t[@]}"; do
	sqlite3 $db<<- EXIT_HERE
	.mode csv
	.headers on
	.output $i.csv
	SELECT * FROM $i;
	.exit
	EXIT_HERE
  echo "$i.csv generated"
done

Download csv file

For more details regarding columns: GEMINI database

Take a look at HG002.g.vcf from last week