SNP calling exercise part 2

Filtering

We have seen that the VCF contains some low-quality or unreliable variant calls. Before downstream analyses, we generally want to remove poor-quality sites or annotate them so they can be excluded later. In this exercise we explore how to apply hard filters and how to remove variants in regions of poor mappability.

Please use the VCF file generated in Part 1.

Hard Filtering

Soft filtering approaches (e.g. VQSR) attempt to statistically learn which variants are “true.” However, these approaches require large cohorts or population-level resources, which may not exist for many organisms or under-sampled populations. For this reason, we often fall back on hard filtering, i.e. applying fixed cutoffs.

Hard filtering is simple but may introduce bias if the filter correlates with variant type (e.g. heterozygous sites often have lower depth). Filters should be chosen thoughtfully.

We will use the following genomic mask file:

/home/databases/databases/GRCh38/mask99.bed.gz

This file is in the BED interval format, which stores genomic regions as:

chromosome   start(0-based)   end(1-based)

• 0-based: first base has coordinate 0
• 1-based: first base has coordinate 1

This mask contains genomic regions to exclude (often low-quality or repetitive regions). Because most genotypers do not recognize duplicated regions, combining hard filtering with mappability filters is best practice.

A typical hard-filtering command using GATK is:

gatk VariantFiltration \
   -V [INPUT VCF] \
   -O [OUTPUT VCF] \
   -filter "DP < 10.0"   --filter-name "DP" \
   -filter "QUAL < 30.0" --filter-name "QUAL30" \
   -filter "SOR > 3.0"   --filter-name "SOR3" \
   -filter "FS > 60.0"   --filter-name "FS60" \
   -filter "MQ < 40.0"   --filter-name "MQ40"

Explanation of filters:

Filter	Meaning
DP < 10	Remove sites with <10× coverage
QUAL < 30	Remove sites where variant quality <30 (variant QUAL ≠ genotype quality GQ — explanation: Variant QUAL vs GQ)
SOR > 3.0	Remove sites with strong strand bias (StrandOddsRatio)
FS > 60	Remove variants failing Fisher Strand bias test (FisherStrand)
MQ < 40	Remove sites where reads have low mapping quality

Note: No filter is perfect — you should progressively add filters, evaluate their impact, and ensure that you do not introduce unwanted biases.

Q1

How many sites were filtered out? Sites that pass all filters have PASS in the 7th column. Use grep to count PASS vs non-PASS entries.

Q2

The 7th column contains the name(s) of the filters that failed. Using cut, sort, and uniq -c, determine which filter removed the most sites.

---

Filtering by Mappability

Next, we remove variants that fall inside low-mappability regions, because reads cannot be uniquely mapped there and false positives are common.

Use bedtools intersect to retain only variants located in high-mappability intervals (≥99% unique mappability):

bedtools intersect -header \
   -a [INPUT VCF] \
   -b /home/databases/databases/GRCh38/filter99.bed.gz \
 | bgzip -c > [OUTPUT VCF]

Name your output:

NA24694_hf_map99.vcf.gz

The "99" refers to the proportion of synthetic reads that map uniquely at that position. Explanation: SNPable / mappability.

Q3

How many variants remain after removing low-mappability regions?

---

Annotation of Variants

Next, we examine the genomic context of variants: intronic, exonic, intergenic, UTR, etc. We use snpEff for variant annotation.

java -jar /home/ctools/snpEff/snpEff.jar eff \
   -dataDir /home/databases/databases/snpEff/ \
   -htmlStats [OUTPUT HTML] \
   GRCh38.99 \
   [INPUT VCF] \
 | bgzip -c > [OUTPUT VCF]

• -dataDir: location of snpEff databases
• GRCh38.99: genome version — must match the reference genome you used earlier

Run snpEff on your hard-filtered VCF (before mappability filtering):

• HTML report: NA24694_hf.html
• Annotated VCF: NA24694_hf_ann.vcf.gz

Open the HTML report:

firefox NA24694_hf.html

Q4

Which genomic region category contains the most variants (exon, intron, upstream, downstream, UTR, etc.)?

Q5

How many variants are predicted to cause a codon change? See explanations at: Point mutation

---

Please find answers here: <a href="SNP_calling_exercise_part_2_answers">SNP_calling_exercise_part_2_answers</a>

Congratulations — you finished the exercise!

Note: When piping bcftools view into other tools, consider specifying the output type using:

-O {b|u|z|v}

This avoids unnecessary compression/decompression and speeds up workflows.