Varlociraptor: Towards a unified theory of variant calling
Johannes Köster
2021
small
structural
SNVs, Indels
Inversions, Duplications, Indels, ...
?
?
?
?
somatic
germline
Variant calling
small
structural
?
?
somatic
germline
Variant calling
tumor/normal
sample, cohort, pedigree
?
?
small
structural
germline
somatic
GATK
Freebayes
Platypus
Pindel
Mutect
Strelka
Octopus
Lancet
Variant calling
Delly
Manta
Gridss
##fileformat=VCFv4.3
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27
20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60
20 1234567 microsat1 GTC G,GTCT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 ##fileformat=VCFv4.3
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 1|0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 2|1:2:0:18,2 2/2:35:4
20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:48:4:51,51 0/0:61:2
20 1234567 microsat1 GTC G,GTCT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/2:17:2 1/1:40:3⊖
Problem:
accumulation of errors
More complex: so far...
small
structural
germline
somatic
Delly
GATK
Freebayes
Platypus
Pindel
Mutect
Strelka
Lancet
Manta
Octopus
arbitrary
Gridss
Approach
Given:
- aligned NGS reads
- candidate variants
Find:
- classification of variants into events (somatic, germline, ...)
- unbiased allele frequency estimates
- while controlling FDR
Model
sample a
sample b
...
contaminates
Model
- insert size
- read sequence
- strand
- read orientation
- ...
Latent variables:
- allele frequency
- fragment from variant allele
- fragment correctly mapped
- biases/artifacts
Observed variables:
Model
ξi∼Bernoulli(αθcτ+(1−α)θhτ)
ωi∼Bernoulli(πi)
Zi∣ξi,ωi=1,β,δ∼
β,δ
- allele frequency
- contamination
- sampling bias
- typing uncertainty
- biases/artifacts (strand, orientation, softclip, homopolymer, ...)
- mapping uncertainty
Model
Pr(Zi∣θc,θh,β)=
⊕ Bayes' theorem
= posterior probability for events
Bayesian FDR control
Variant calling grammar
species:
heterozygosity: 0.001
ploidy:
male:
all: 2
X: 1
Y: 1
female:
all: 2
X: 2
Y: 0
samples:
cersei:
sex: female
events:
het: "cersei:0.5"
hom: "cersei:1.0"
image: gameofthrones.fandom.com
Variant calling grammar
species:
heterozygosity: 0.001
ploidy:
male:
all: 2
X: 1
Y: 1
female:
all: 2
X: 2
Y: 0
samples:
cersei:
sex: female
jamie:
sex: female
joffrey:
sex: male
inheritance:
mendelian:
mother: cersei
father: jamieevents:
jamie: "jamie:0.5 & jamie:1.0"
cersei: "cersei:0.5 & cersei:1.0"
denovo_joffrey: "(joffrey:0.5 | joffrey:1.0) & !$cersei & !$jamie"
images: gameofthrones.fandom.com
Variant calling grammar
species:
heterozygosity: 0.001
ploidy:
male:
all: 2
X: 1
Y: 1
female:
all: 2
X: 2
Y: 0
samples:
cersei:
sex: female
somatic-effective-mutation-rate: 1e-10
tumor:
inheritance:
clonal:
from: cersei
contamination:
by: cersei
fraction: 0.1
somatic-effective-mutation-rate: 1e-6events:
germline: "cersei:0.5 | cersei:1.0"
somatic: "cersei:]0.0,0.5["
somatic_tumor_low: "cersei:0.0 & tumor:]0.0,0.1["
somatic_tumor_high: "cersei:0.0 & tumor:[0.1,1.0]"
images: gameofthrones.fandom.com, 123rf.com
Variant calling grammar
samples:
cersei:
sex: female
somatic-effective-mutation-rate: 1e-10
tumor:
inheritance:
clonal:
from: cersei
contamination:
by: cersei
fraction: 0.1
somatic-effective-mutation-rate: 1e-6
relapse:
inheritance:
clonal:
from: cersei
contamination:
by: cersei
fraction: 0.2
somatic-effective-mutation-rate: 1e-6expressions:
somatic_tumor: "cersei:0.0 & tumor:]0.0,1.0]"
events:
germline: "cersei:0.5 | cersei:1.0"
somatic: "cersei:]0.0,0.5["
somatic_tumor_no_increase: "$somatic_tumor & l2fc(relapse,tumor) < 1"
somatic_tumor_increase: "$somatic_tumor & l2fc(relapse,tumor) >= 1"
somatic_relapse: "cersei:0.0 & tumor:0.0 & relapse:]0.0,1.0]"
images: gameofthrones.fandom.com, 123rf.com
Variant calling grammar
species:
heterozygosity: 0.001
ploidy:
male:
all: 2
X: 1
Y: 1
female:
all: 2
X: 2
Y: 0
samples:
cersei_dna_illumina:
sex: female
cersei_dna_nanopore:
sex: female
inheritance:
clonal:
from: cersei
cersei_rna_illumina:
universe: [0.0,1.0]events:
het: "cersei_dna_illumina:0.5 & cersei_dna_nanopore:0.5 & cersei_rna_illumina:]0.0,1.0]"
hom: "cersei_dna_illumina:1.0 & cersei_dna_nanopore:1.0 & cersei_rna_illumina:1.0"
White box model
chr1 1200 . GAAGCCATCTGCCCCATGGGAAAGTGACTTAAT G . .
SVLEN=-32;PROB_SOMATIC_TUMOR=0.110292;PROB_ABSENT=202.181;PROB_GERMLINE_HET=61.7704;
PROB_SOMATIC_NORMAL=16.0075;PROB_GERMLINE_HOM=4346.23;PROB_ARTIFACT=103.395
DP:AF:OBS:SB:ROB
36:0:11Np+<10Np+>8Np-<7Np->:.:.
48:0.102431:14Np+<13Np-<11Np->8Np+>1Vp-<1Vp+>:.:.- posterior allele frequency estimate
- observations (count, posterior odds, type, strand, read orientation, ...)
- bias/artifact estimates
- event probabilities (PHRED scaled)
https://varlociraptor.github.io
Varlociraptor: Towards a unified theory of variant calling Johannes Köster 2021