Difference between revisions of "Analyses of Indels"
Line 34: | Line 34: | ||
You may also work with vcf.gz, just name the output as *.vcf.gz. But it will be slower with vt. | 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" | ||
==Normalization== | ==Normalization== | ||
Line 46: | Line 84: | ||
number of duplicate variants found for some of the 1000 Genomes Trio High Coverage call sets. | 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 | Although left alignment seems to be a trivial concept, it is easily overlooked and remain a common | ||
− | mistake. | + | mistake. |
− | |||
− | |||
{| class="wikitable" | {| class="wikitable" | ||
Line 121: | Line 157: | ||
| 7541 | | 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. | ||
+ | |||
+ | |||
==Coding regions== | ==Coding regions== |
Revision as of 15:38, 20 February 2014
Motivation
This wiki page details some standard Indel analyses which hopefully can help the group in understanding the issues and perform the analyses quickly without reinventing the wheel.
Tools
You can download vt and have some working knowledge of PERL to do stuff that vt does not support.
Analyses
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"
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.
Coding regions
The proportion of frameshift Indels amongst coding region indels is a potential indicator of quality.
You can obtain it by using the profile_indels analysis.
STR
Annotation of STRs is really important. Show example of a deceptive single base pair variant
Annotation of Indels
Examining Mendelian Errors
Useful to have call sets from several different callers
Concordance
Can check concordance of genotypes between callers
Overlapping percentages with known data sets
With Mills with dbSNP with exome chips with genotyping chips if available
Useful stratifying features
AF - rare versus common Indel length - computed naively versus tract length Allele frequency bins Type of Indels - homopolymer types and STR types and isolated Adjacent SNPs Adjacent MNPs Clumping variants
Other useful evaluations
genotype likelihood concordance concordance stratified by indel length or tract length mendelian concordance by tract length