1,102
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchSequencing Workshop Analysis of Indels (view source)
Revision as of 16:47, 16 June 2014
, 16:47, 16 June 2014→Normalization
##INFO=<ID=FIC,Number=1,Type=Float,Description="Genotype likelihood based Inbreeding Coefficient">
##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">
##INFO=<ID=AB,Number=1,Type=Float,Description="Genotype likelihood based Allele Balance">
##FILTER=<ID=TPASS,Description="Temporary pass">
##FILTER=<ID=PASS,Description="Temporary pass">
##FILTER=<ID=overlap,Description="Overlapping variant">
##FILTER=<ID=overlap,Description="Overlapping variant">
The columns are CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, Genotype fields denoted by the sample name.
The columns are CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, Genotype fields denoted by the sample name.
22 36990877 . GGT G . TPASS AC=32;AN=116;AF=0.275862;GC=32,20,6;GN=58;
22 36990878 . GGT G 455 PASS AC=32;AN=116;AF=0.275862;GC=32,20,6;GN=58;
GF=0.551724,0.344828,0.103448;NS=58;
GF=0.551724,0.344828,0.103448;NS=58;
HWEAF=0.275797;HWEGF=0.52447,0.399466,0.0760642;
HWEAF=0.275797;HWEGF=0.52447,0.399466,0.0760642;
GT:PL:DP:AD:GQ 0/0:0,9,108:9:3,0,6:10
GT:PL:DP:AD:GQ 0/0:0,9,108:9:3,0,6:10
Let's look at the first 5 records.
Let's look at the record's fields.
22 : chromosome
22 : chromosome
36990878 : genome position
. : this is the ID field that is left blank.
. : this is the ID field that is left blank.
GGT : the reference sequence that is replaced by the alternative sequence below.
GGT : the reference sequence that is replaced by the alternative sequence below.
G : so this is basically a deletion of GT
G : so this is basically a deletion of GT
455 : QUAL field denoting validity of this variant, higher the better.
PASS : a passed variant.
INFO : fields containing information about the variant.
INFO : fields containing information about the variant.
FORMAT : format field labels for the genotype columns.
FORMAT : format field labels for the genotype columns.
0/0:0,9,108:9:3,0,6:10 : genotype information.
0/0:0,9,108:9:3,0,6:10 : genotype information.
The information field are as follows:
The information field are as follows:
AC=32 : alternate allele count
AC=32 : alternate allele count
AN=116 : total number of alleles
AN=116 : total number of alleles
AF=0.275862 : allele frequency based on AC/AN
AF=0.28 : allele frequency based on AC/AN
GC=32,20,6 : genotype counts for 0/0, 0/1, 1/1
GC=32,20,6 : genotype counts for 0/0, 0/1, 1/1
GN=58; : total number of genotypes
GN=58; : total number of genotypes
GF=0.551724,0.344828,0.103448 : genotype frequencies
GF=0.55,0.34,0.10 : genotype frequencies based on GC
NS=58 : no. of samples.
NS=58 : no. of samples
HWEAF=0.275797: Genotype likelihood based estimation of the allele frequency assuming Hardy Weinberg equilibrium.
HWEAF=0.28 : genotype likelihood based estimation of the allele frequency assuming Hardy Weinberg equilibrium
HWEGF=0.52447,0.399466,0.0760642 : genotype frequency derived from HWEAF
HWEGF=0.52,0.40,0.08 : genotype frequency derived from HWEAF
MLEAF=0.27366: genotype likelihood based estimation of the genotype frequency
MLEAF=0.27 : genotype likelihood based estimation of the genotype frequency
MLEGF=0.494275,0.464129,0.0415952 : genotype frequency derived from MLEAF
MLEGF=0.49,0.46,0.04 : genotype frequency derived from MLEAF
HWE_LLR=-0.453098 : log likelihood ratio of HWE test
HWE_LLR=-0.45 : log likelihood ratio of HWE test
HWE_LPVAL=-1.0755 : log p value of HWE test
HWE_LPVAL=-1.08 : log p value of HWE test
HWE_DF=1 : degrees of freedom of the HWE test
HWE_DF=1 : degrees of freedom of the HWE test
FIC=-0.0718807;AB=0.6129
FIC=-0.07 : genotype likelihood based inbreeding coefficient
AB=0.61 : genotype likelihood based allele balance
===GENOTYPE field===
===GENOTYPE field===
0/0: homozygous reference chosen based on in PL.
The genotype fields are described as follows:
0,9,108: PHRED scaled genotype likelihoods
9: no. of reads covering this variant
0/0 : homozygous reference chosen based on PL
3,0,6: allele depth
0,9,108 : PHRED scaled genotype likelihoods
9 : no. of reads covering this variant
3,0,6 : allele depth
counts of reads supporting the reference allele,
10 : genotype quality.
the alternative allele and neither alleles respectively.
The last category might be due to insufficient
coverage of the read over the locus
or simply an allele that is not accounted for.
10 : genotype quality
= Analysis =
= Analysis =
The following section details some simple analyses we can perform.
=Analyses=
== Summary ==
First you want to know what is in the bcf file.
vt peek all.genotypes.bcf
vt peek all.genotypes.bcf
stats: no. of samples : 62
[o] output VCF file -
no. of chromosomes : 1 <br>
no. Indels : 720
2 alleles (ins/del) : 720 (0.84) [328/392] #(insertion deletion ratio) [#insertions, #deletions]
>=3 alleles (ins/del) : 0 (-nan) [0/0] <br>
no. of observed variants : 720
The variants have filter labels PASS meaning a temporary pass and overlap, meaning that the variants are overlapping with another variant, implying multiallelicity.
We can count the number of variants with the following commands.
vt peek all.genotypes.bcf -f "FILTER.PASS"
stats: no. of samples : 62
no. of chromosomes : 1 <br>
no. Indels : 584
2 alleles (ins/del) : 584 (0.69) [239/345]
>=3 alleles (ins/del) : 0 (-nan) [0/0]
vt peek all.genotypes.bcf -f "FILTER.overlap"
stats: no. of samples : 62
no. of chromosomes : 1 <br>
no. Indels : 136
2 alleles (ins/del) : 136 (1.89) [89/47] #notice the difference in insertion deletion ratios
>=3 alleles (ins/del) : 0 (-nan) [0/0]
#passed singletons only
vt peek all.genotypes.bcf -f "FILTER.PASS&&INFO.AC==1"
#passed indels of length 1 only
vt peek all.genotypes.bcf -f "FILTER.PASS&&LEN==1"
#passed indels of length >4
vt peek all.genotypes.bcf -f "FILTER.PASS&&LEN>1"
#passed singletons of length 4 or insertions of length 3
vt peek all.genotypes.bcf -f "FILTER.PASS&&(LEN==4||DLEN==3)"
== Comparison with other data sets ==
It is usually useful to examine the call sets against known data sets for the passed variants.
vt profile_indels -g indel.reference.txt -r hs37d5.fa all.genotypes.bcf -i 22:36000000-37000000 -f "PASS"
data set
No Indels : 613 [0.72] #613 passed variants with an insertion deletion ratio of 0.72
FS/NFS : 0.50 (2/2) #frame shift / non frameshift indels proportion, the bracket gives the counts of the frame shift and non frameshift indels
Low complexity : 0.46 (283/613) #fraction of indels in low complexity region, the bracket gives the counts of the indels <br>
1000G #1000 Genomes Phase 1 data set
A-B 371 [0.76] #variants found in call set only, square brackets contain insertion deletion ratio
A&B 242 [0.66] #variants found in both data sets
B-A 276 [0.89] #variants found in 1000G phase 1 data set only
Precision 39.5% #39.5% of the call set are previously known, so 60.5% are novel variants.
Sensitivity 46.7% #sensitivity of variant calling, 46,7% of known variants from 1000 Genomes were rediscovered <br>
mills #The gold standard Mills et al. indel set
A-B 542 [0.68]
A&B 71 [1.03]
B-A 31 [1.07]
Precision 11.6%
Sensitivity 69.6% <br>
dbsnp #Indels from dbSNP
A-B 405 [0.68]
A&B 208 [0.79]
B-A 494 [2.03]
Precision 33.9%
Sensitivity 29.6%
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.
Complexity region analysis: Indels in regions marked by DUST - a low complexity region masker used in the NCBI pipeline.
Overlap analysis: overlap analysis with other data sets is an indicator of sensitivity.
* 1000G: contains Indels from 1000 Genomes, represent a wide spectrum of variants from many different populations. Variants here have an allele frequency above 0.005.
* 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.
* 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.
We perform the same analysis for the failed variants again, the relatively low overlap with known data sets imply a reasonable tradeoff in sensitivity and specificity.
profile_indels v0.5
vt profile_indels -g indel.reference.txt -r hs37d5.fa all.genotypes.bcf -i 22:36000000-37000000 -f "~PASS"
data set
data set
No Indels : 720 [0.84] #720 indels, with and insertion deletion ratio of 0.84
No Indels : 107 [2.06]
FS/NFS : 0.50 (2/2) #only 4 variants overlap with coding regions, half of which are frameshift variants
FS/NFS : -nan (0/0)
Low complexity : 0.47 (335/720) #47% of the variants are in low complexity regions <br>
Low complexity : 0.79 (85/107) <br>
1000G
1000G
A-B 719 [0.83]
A-B 107 [2.06]
A&B 0 [-nan]
B-A 518 [0.77]
A&B 0 [-nan] #no variants overlaps with Mills et al. double hit variants.
B-A 102 [1.04]
Precision 0.0%
Precision 0.0%
Sensitivity 0.0% <br>
Sensitivity 0.0% <br>
mills
A-B 105 [2.09]
A&B 2 [1.00]
B-A 100 [1.04]
Precision 1.9%
Sensitivity 2.0% <br>
dbsnp
dbsnp
A-B 720 [0.84]
A-B 102 [2.00]
A&B 0 [-nan] #no variants overlaps with Mills et al. double hit variants.
A&B 5 [4.00]
B-A 702 [1.52]
B-A 697 [1.51]
Precision 0.0%
Precision 4.7%
Sensitivity 0.0%
Sensitivity 0.7%
This analysis supports filters too.
==Normalization==
A slight digression here, when analyzing indels, it is important to normalize it. While it is a simple concept,
it is hardly standardized. The call set here had already been normalized but we feel that this is an important
concept so we discuss this a bit here.
Indel representation is not unique, you should normalize them and remove duplicates.
Indel representation is not unique, you should normalize them and remove duplicates.
| 0
| 0
| 374
| 374
|
| 0
|
| 0
|-
|-
| Left aligned
| Left aligned
To normalize and remove duplicate variants:
To normalize and remove duplicate variants:
vt normalize mills.genotypes.bcf -r ~/ref/vt/grch37/hs37d5.fa | vt mergedups - -o mills.normalized.genotypes.bcf
vt normalize mills.genotypes.bcf -r 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.
and you will observe that 3994 variants had to be left aligned and 1092 variants were removed.
no. right trimmed : 0 <br>
no. right trimmed : 0 <br>
no. variants observed : 9996 <br>
no. variants observed : 9996 <br>
<br>
stats: Total number of observed variants 9996
stats: Total number of observed variants 9996
Total number of unique variants 8904 <br>
Total number of unique variants 8904 <br>
UMICH's algorithm for normalization has been adopted by Petr Danecek in bcftools and is also used in GKNO.
