Repository navigation
Exome III: Alignment processing and variant discovery
When dealing with variant calling, an important step is to identify possible PCR duplicates, this to remove any bias when we are doing the actual calling, for this use MarkDuplicates from picard tools. Once you have your new BAM file, don't forget to index it.
java -jar $picard MarkDuplicates
Generate some metrics for the most recent BAM files from both samples either with samtools (flagstat and idxstats) or picard (CollectAlignmentSummaryMetrics and BamIndexStats). Run multiqc ., to compile the results, and have a look at the report.
Q. How many duplicates do the samples have? Which sample has more sequencing reads?
The Genome Analysis Toolkit or GATK is a software package developed to analyze
next-generation resequencing data, focusing on variant discovery and genotyping.
For the variant calling we will use the HaplotypeCaller, which is an SNP/indel
caller that uses a Bayesian genotype likelihood model to estimate simultaneously
the most likely genotypes and allele frequency in a population of N samples.
Use $GenomeAnalysisTK with these flags:
- -T HaplotypeCaller
- -R $chr11_ref
- --max_alternate_alleles 1
Besides the required ones.
java -jar $GenomeAnalysisTK -T HaplotypeCaller -h
Inspect the VCF file you just created and get familiar with its information.
Q. How do you differentiate between a homozygous mutations from an heterozygous? How many SNPs do you have in each sample? Are there any insertions/deletions?
Let's convert the VCF file into a tab delimited file. Use -T VariantsToTable from $GenomeAnalysisTK and add the following Filters:
-F CHROM -F POS -F ID -F REF -F ALT -F QUAL -F FILTER -F AC -GF GT -GF AD -GF DP
We will use this file later. Have a look and understand the output.
SNP-callers emit a number of false-positive mutations, as any prediction program. Some of these false-positives can be
detected and rejected by various statistical tests. Filtering or flagging the variant calls
in the VCF file is a useful step to detect such false-positives. For this run GATK's Variant Filtration using only this filtering expression:
- --filterExpression "DP < 10.0" --filterName "LowDP"
These are hard filters that are applied when a variant callset is too small or when training sets are not available. Try to understand the filtering done above. You may find some info at the beginning of the VCF file.
Q. How many variants have passed our filtering criteria for both files?
Previous: Exome II: QC, filtering and mapping
Created by Marcela Dávila