Analyses of Indels
Motivation
This wiki page details some standard Indel analyses for the sequencing workshop in the example indel data set.
Tools
This walkthrough requires vt.
Analyses
The file generated from the indel calling is a binary version [BCFv2.1] of the Variant Call Format (VCF). BCFv2.1 is more efficient to process as the data is already stored in computer readable format on the hard disk. It is however not necessarily more compact than VCF4.2 especially when the format fields are rich in details.
Anatomy of all.genotypes.bcf
You can access the header by running the command:
vt view -H all.genotypes.bcf.
The header is as follows:
##fileformat=VCFv4.2 ##FILTER=<ID=PASS,Description="All filters passed"> ##contig=<ID=22,length=51304566> ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype"> ##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes"> ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Depth"> ##FORMAT=<ID=AD,Number=3,Type=Integer,Description="Allele Depth"> ##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality"> ##INFO=<ID=AC,Number=A,Type=Integer,Description="Alternate Allele Counts"> ##INFO=<ID=AN,Number=1,Type=Integer,Description="Total Number Allele Counts"> ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data"> ##INFO=<ID=AF,Number=A,Type=Float,Description="Alternate Allele Frequency"> ##INFO=<ID=GC,Number=G,Type=Integer,Description="Genotype Counts"> ##INFO=<ID=GN,Number=1,Type=Integer,Description="Total Number of Genotypes Counts"> ##INFO=<ID=GF,Number=G,Type=Float,Description="Genotype Frequency"> ##INFO=<ID=HWEAF,Number=A,Type=Float,Description="Genotype likelihood based MLE Allele Frequency assuming HWE"> ##INFO=<ID=HWEGF,Number=G,Type=Float,Description="Genotype likelihood based MLE Genotype Frequency assuming HWE"> ##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Genotype likelihood based MLE Allele Frequency"> ##INFO=<ID=MLEGF,Number=G,Type=Float,Description="Genotype likelihood based MLE Genotype Frequency"> ##INFO=<ID=HWE_LLR,Number=1,Type=Float,Description="Genotype likelihood based Hardy Weinberg ln(Likelihood Ratio)"> ##INFO=<ID=HWE_LPVAL,Number=1,Type=Float,Description="Genotype likelihood based Hardy Weinberg Likelihood Ratio Test Statistic ln(p-value)"> ##INFO=<ID=HWE_DF,Number=1,Type=Integer,Description="Degrees of freedom for Genotype likelihood based Hardy Weinberg Likelihood Ratio Test Statistic"> ##INFO=<ID=FIC,Number=1,Type=Float,Description="Genotype likelihood based Inbreeding Coefficient"> ##INFO=<ID=AB,Number=1,Type=Float,Description="Genotype likelihood based Allele Balance"> ##FILTER=<ID=TPASS,Description="Temporary pass"> ##FILTER=<ID=overlap,Description="Overlapping variant">
body
To view some of the records:
22 36990877 . GGT G . TPASS AC=32;AN=116;AF=0.275862;GC=32,20,6;GN=58;GF=0.551724,0.344828,0.103448;NS=58;
HWEAF=0.275797;HWEGF=0.52447,0.399466,0.0760642;MLEAF=0.27366;
MLEGF=0.494275,0.464129,0.0415952;HWE_LLR=-0.453098;HWE_LPVAL=-1.0755;HWE_DF=1;FIC=-0.0718807;AB=0.6129
GT:PL:DP:AD:GQ 0/0:0,9,108:9:3,0,6:10 22 36991203 . TGAG T . TPASS AC=5;AN=124;AF=0.0403226;GC=58,3,1;GN=62;GF=0.935484,0.0483871,0.016129;NS=62;HWEAF=0.0355594;HWEGF=0.930145,0.0685899,0.00126447;MLEAF=0.0353706;MLEGF=0.929259,0.0707412,5.94815e-11;HWE_LLR=-0.0443401;HWE_LPVAL=-0.266754;HWE_DF=1;FIC=-0.0109029;AB=0.562243 GT:PL:DP:AD:GQ 0/0:0,12,155:6:4,0,2:12 22 36995311 . GA G . TPASS AC=61;AN=124;AF=0.491935;GC=21,21,20;GN=62;GF=0.33871,0.33871,0.322581;NS=62;HWEAF=0.492227;HWEGF=0.257834,0.499879,0.242287;MLEAF=0.492028;MLEGF=0.298019,0.419905,0.282076;HWE_LLR=-0.605122;HWE_LPVAL=-1.30459;HWE_DF=1;FIC=0.0444598;AB=0.53981 GT:PL:DP:AD:GQ 0/1:55,0,24:3:1,2,0:24 22 36995329 . GA G . TPASS AC=2;AN=124;AF=0.016129;GC=60,2,0;GN=62;GF=0.967742,0.0322581,0;NS=62;HWEAF=0.0164696;HWEGF=0.967332,0.0323966,0.000271246;MLEAF=0.0165148;MLEGF=0.96697,0.0330296,7.28675e-44;HWE_LLR=-0.0171385;HWE_LPVAL=-0.158856;HWE_DF=1;FIC=-0.00275028;AB=0.339746 GT:PL:DP:AD:GQ 0/0:0,9,97:5:3,0,2:10
File Preparation
The VCF file you work with should preferably be BCF2.1 compatible. Here we provide an example in /net/fantasia/home/atks/indel_analysis_tutorial. \\
To convert to BCF format which will work fast with vt:
vt view mills.vcf -o mills.bcf
You will encounter an error as the header does not contain contigs. To fix this, you should construct a complete header for mills.vcf. This is done for you in mills.with.alt.hdr.*
vt view mills.with.alt.hdr.vcf -o mills.genotypes.bcf
To index:
vt index mills.genotypes.bcf
To extract just the site list which is convenient for working with if you are not analysing the genotypes of the individuals
vt view -s mills.genotypes.bcf -o mills.sites.bcf
To index:
vt index mills.sites.bcf
You may also work with vcf.gz, just name the output as *.vcf.gz. But it will be slower with vt.
Peek
You can see what you have in the file with:
vt peek mills.genotypes.bcf
You can also focus on a chromosome:
vt peek mills.genotypes.bcf -i 20
Or with just passed variants:
vt peek mills.genotypes.bcf -i 20 -f PASS
Or with failed variants:
vt peek mills.genotypes.bcf -i 20 -f ~PASS
Or with just 1bp indels:
vt peek mills.genotypes.bcf -i 20 -f "PASS&&DLEN==1"
Or with just 1bp deletions:
vt peek mills.genotypes.bcf -i 20 -f "PASS&&LEN==-1"
Or with just biallelic 1bp indels:
vt peek mills.genotypes.bcf -i 20 -f "PASS&&N_ALLELE==2&&LEN==1"
Or with just biallelic 1bp indels that are somewhat rare:
vt peek mills.sites.bcf -f "PASS&&N_ALLELE==2&&LEN==1&&INFO.AF<0.03"
Or with just biallelic 1bp indels that are somewhat rare with sanity checking:
vt peek mills.sites.bcf -f "PASS&&N_ALLELE==2&&LEN==1&&INFO.AC/INFO.AN<0.03"
which you will observe discrepancies due to rounding off in AF. So you should probably use INFO.AC/INFO.AN.
Normalization
Indel representation is not unique, you should normalize them and remove duplicates.
Variant normalization is implemented in vt and this page explains the algorithm and also provides a simple proof of correctness - Variant Normalization
The following table shows the number of variants that had to be normalized and the corresponding type of normalization performed and the ensuing number of duplicate variants found for some of the 1000 Genomes Trio High Coverage call sets. Although left alignment seems to be a trivial concept, it is easily overlooked and remain a common mistake.
| Dataset | Freebayes | Haplotyecaller | PINDEL | Platypus | RTG | Samtools | SGA |
|---|---|---|---|---|---|---|---|
| Biallelic | |||||||
| Left trim | 27069 | 1 | 0 | 0 | 0 | 0 | 15047 |
| Left aligned | 3 | 1 | 1 | 0 | 12262 | 2 | 1892 |
| Multi-allelic | |||||||
| Left trim | 40782 | 0 | 0 | 0 | 374 | ||
| Left aligned | 1892 | 0 | 0 | 0 | 1329 | 1 | 0 |
| Right trimmed | 0 | 0 | 0 | 25393 | 0 | 11 | 0 |
| Duplicate variants | 0 | 1 | 155 | 3143 | 286 | 8 | 7541 |
Another example is the Mills et al. data set which followed up with 10004 Indels for validation. Out of 9996 passed variants, it was found that after normalization, only 8904 distinct Indels remain - about a loss of 11% of variant thought distinct.
To normalize and remove duplicate variants:
vt normalize mills.genotypes.bcf -r ~/ref/vt/grch37/hs37d5.fa | vt mergedups - -o mills.normalized.genotypes.bcf
and you will observe that 3994 variants had to be left aligned and 1092 variants were removed.
stats: biallelic
no. left trimmed : 0
no. left trimmed and left aligned : 0
no. left trimmed and right trimmed : 0
no. left aligned : 3994
no. right trimmed : 0
multiallelic
no. left trimmed : 0
no. left trimmed and left aligned : 0
no. left trimmed and right trimmed : 0
no. left aligned : 0
no. right trimmed : 0
no. variants observed : 9996
Time elapsed: 0.14s
stats: Total number of observed variants 9996
Total number of unique variants 8904
Time elapsed: 0.13s
The following will be slight faster: + denotes using of uncompressed bcf stream.
vt normalize mills.genotypes.bcf -r ~/ref/vt/grch37/hs37d5.fa -o + | vt mergedups + -o mills.normalized.genotypes.bcf
Also remember to index this file and extract the sites.
Insertion/Deletion ratios, Coding Regions and Overlap analysis
You can obtain measure of insertion deletion ratios, coding region indels and sensitivity analysis by using the profile_indels analysis.
vt profile_indels -g indel.reference.txt -r ~/ref/vt/grch37/hs37d5.fa mills.normalized.sites.bcf
The indel.reference.txt file contains the required reference to perform the overlap analysis.
data set
No Indels : 8904 [0.93] //#variants in your data set [ins/del ratio]
FS/NFS : 0.66 (67/35) //Proportion of frameshift Indels. (#Frameshift Indels/#Nonframeshift Indels)
dbsnp //A represents the data set you input, B represents dbsnp
A-B 2975 [1.06] //#variants in A only [ins/del ratio]
A&B 5929 [0.86] //#variants in A and B
B-A 2059845 [1.51]
Precision 66.6% //A&B/A this represents how novel your data set is in the variants represented.
Sensitivity 0.3% //A&B/B this represents sensitivity somewhat if dbsnp is considered a high quality Indel
//set and the sample are the same in both data sets. (which they usually are not, this is still
//nonetheless a useful indicator)
mills
A-B 5705 [0.81]
A&B 3199 [1.18]
B-A 203819 [0.98]
Precision 35.9%
Sensitivity 1.5%
mills.chip
A-B 0 [-nan]
A&B 8904 [0.93]
B-A 0 [-nan]
Precision 100.0%
Sensitivity 100.0%
affy.exome.chip
A-B 8821 [0.93]
A&B 83 [0.69]
B-A 34011 [0.47]
Precision 0.9%
Sensitivity 0.2%
Ins/Del ratios: Reference alignment based methods tend to be biased towards the detection of deletions. This provides a useful measure for discovery Indel sets to show the varying degree of biasness. It also appears that as coverage increases, the ins/del ratio tends to 1.
Coding region analysis: Coding region Indels may be categorised as Frame shift Indels and Non frameshift Indels. A lower proportion of Frameshift Indels may indicate a better quality data set but this depends also on the individuals sequenced.
Overlap analysis: overlap analysis with other data sets is an indicator of sensitivity.
- dbsnp: contains Indels submitted from everywhere, I am not sure what does this represent exactly. But assuming most are real, then precision is a useful estimated quantity from this reference data set.
- Mills: contains doublehit common indels from the Mills. et al paper and is a relatively good measure of sensitivity for common variants. Because not all Indels in this set is expected to be present in your sample, this actually gives you an underestimate of sensitivity.
- Mills chip: This is a subset of the Mills data set. There are genotypes here that are useful for subsetting polymophic subsets of variants that are present in samples common with your data set, this can potentially provide a better estimate of sensitivity. In general not very useful unless you happen to be working on 1000 Genomes data or any data set who's individuals are commonly studied.
- Affy Exome Chip: This contains somewhat rare variants in exonic regions and is useful for exome chip analysis. You should subset your exome data to exome region Indels before comparing against this data set.
This analysis supports filters too.
to document
- Annotation of STRs is really important. Show example of a deceptive single base pair variant
- Mendelian analysis
- AFS
- Can check concordance of genotypes between callers - partitiion
- Type of Indels - homopolymer types and STR types and isolated, Adjacent SNPs ,Adjacent MNPs,Clumping variants
- genotype likelihood concordance
- concordance stratified by indel length or tract length
- mendelian concordance by tract length
