Skip to contents

wgsLR: Shotgun sequencing for human identification: Dynamic SNP marker sets and likelihood ratio calculations accounting for errors

Please refer to the online documentation at https://mikldk.github.io/wgsLR, including the vignettes.

Scientific publication

The research associated with this software is described in

Andersen, M. M., Kampmann, M.-L., Jepsen, A. H., Morling, N., Eriksen, P. S., Børsting, C., & Andersen, J. D. (2025). Shotgun DNA sequencing for human identification: Dynamic SNP selection and likelihood ratio calculations accounting for errors. Forensic Science International: Genetics, 74, 103146. doi:10.1016/j.fsigen.2024.103146.

Installation

To install from Github with vignettes run this command from within R (please install remotes first if not already installed):

# install.packages('remotes')
remotes::install_github("mikldk/wgsLR", build_vignettes = TRUE)

You can also install the package without vignettes if needed as follows:

remotes::install_github("mikldk/wgsLR")

A few small examples

Estimating the genotype error probability, ww 

cases <- wgsLR::sample_data_Hp_w(n = 1000, w = 0.1, p = c(0.25, 0.25, 0.5))
tab <- table(cases$xT, cases$xR)
tab
##    
##       0   1   2
##   0 171  62   4
##   1  65 193 100
##   2  12  72 321
w_mle <- wgsLR::estimate_w(tab)
w_mle
## [1] 0.1012883

Cautionary note: not just standard VCF files

It is necessary to obtain sequencing results for all bases in the selected segments (including read depth and genotype quality). Thus, it is not sufficient to just use information from “confirmed”/high probability variants from the reference genome (variants identified in standard vcf-file format), as this can introduce bias in the results. Information from all bases in the chosen genomic areas of interest is needed. One way to achieve this by using GATK HaplotypeCaller with the additional argument --emit-ref-confidence BP_RESOLUTION for the genomic areas of interest (using -L areas.interval_list).

Calculating likelihood ratios (LRLR’s)

Same genotyping error probability for both trace sample (xt=x_t =xT) and reference sample (xr=x_r =xR):

tab <- wgsLR::d_LR_w[, c("xT", "xR")]
tab$LR <- paste0("$", wgsLR::d_LR_w$expr_tex, "$")
tab |> 
  kbl(format = "markdown")
xT xR LR
0 0 p0(1−w)4+p1w2(1−w)2+p2w4p02(1−w)4+2p0p1w(1−w)3+2p0p2w2(1−w)2+p12w2(1−w)2+2p1p2w3(1−w)+p22w4\frac{p_{0} \left(1 - w\right)^{4} + p_{1} w^{2} \left(1 - w\right)^{2} + p_{2} w^{4}}{p_{0}^{2} \left(1 - w\right)^{4} + 2 p_{0} p_{1} w \left(1 - w\right)^{3} + 2 p_{0} p_{2} w^{2} \left(1 - w\right)^{2} + p_{1}^{2} w^{2} \left(1 - w\right)^{2} + 2 p_{1} p_{2} w^{3} \left(1 - w\right) + p_{2}^{2} w^{4}}
0 1 2p0w(1−w)3+p1w3(1−w)+p1w(1−w)3+2p2w3(1−w)2p02w(1−w)3+3p0p1w2(1−w)2+p0p1(1−w)4+2p0p2w3(1−w)+2p0p2w(1−w)3+p12w3(1−w)+p12w(1−w)3+p1p2w4+3p1p2w2(1−w)2+2p22w3(1−w)\frac{2 p_{0} w \left(1 - w\right)^{3} + p_{1} w^{3} \left(1 - w\right) + p_{1} w \left(1 - w\right)^{3} + 2 p_{2} w^{3} \left(1 - w\right)}{2 p_{0}^{2} w \left(1 - w\right)^{3} + 3 p_{0} p_{1} w^{2} \left(1 - w\right)^{2} + p_{0} p_{1} \left(1 - w\right)^{4} + 2 p_{0} p_{2} w^{3} \left(1 - w\right) + 2 p_{0} p_{2} w \left(1 - w\right)^{3} + p_{1}^{2} w^{3} \left(1 - w\right) + p_{1}^{2} w \left(1 - w\right)^{3} + p_{1} p_{2} w^{4} + 3 p_{1} p_{2} w^{2} \left(1 - w\right)^{2} + 2 p_{2}^{2} w^{3} \left(1 - w\right)}
0 2 p0w2(1−w)2+p1w2(1−w)2+p2w2(1−w)2p02w2(1−w)2+p0p1w3(1−w)+p0p1w(1−w)3+p0p2w4+p0p2(1−w)4+p12w2(1−w)2+p1p2w3(1−w)+p1p2w(1−w)3+p22w2(1−w)2\frac{p_{0} w^{2} \left(1 - w\right)^{2} + p_{1} w^{2} \left(1 - w\right)^{2} + p_{2} w^{2} \left(1 - w\right)^{2}}{p_{0}^{2} w^{2} \left(1 - w\right)^{2} + p_{0} p_{1} w^{3} \left(1 - w\right) + p_{0} p_{1} w \left(1 - w\right)^{3} + p_{0} p_{2} w^{4} + p_{0} p_{2} \left(1 - w\right)^{4} + p_{1}^{2} w^{2} \left(1 - w\right)^{2} + p_{1} p_{2} w^{3} \left(1 - w\right) + p_{1} p_{2} w \left(1 - w\right)^{3} + p_{2}^{2} w^{2} \left(1 - w\right)^{2}}
1 0 2p0w(1−w)3+p1w3(1−w)+p1w(1−w)3+2p2w3(1−w)2p02w(1−w)3+3p0p1w2(1−w)2+p0p1(1−w)4+2p0p2w3(1−w)+2p0p2w(1−w)3+p12w3(1−w)+p12w(1−w)3+p1p2w4+3p1p2w2(1−w)2+2p22w3(1−w)\frac{2 p_{0} w \left(1 - w\right)^{3} + p_{1} w^{3} \left(1 - w\right) + p_{1} w \left(1 - w\right)^{3} + 2 p_{2} w^{3} \left(1 - w\right)}{2 p_{0}^{2} w \left(1 - w\right)^{3} + 3 p_{0} p_{1} w^{2} \left(1 - w\right)^{2} + p_{0} p_{1} \left(1 - w\right)^{4} + 2 p_{0} p_{2} w^{3} \left(1 - w\right) + 2 p_{0} p_{2} w \left(1 - w\right)^{3} + p_{1}^{2} w^{3} \left(1 - w\right) + p_{1}^{2} w \left(1 - w\right)^{3} + p_{1} p_{2} w^{4} + 3 p_{1} p_{2} w^{2} \left(1 - w\right)^{2} + 2 p_{2}^{2} w^{3} \left(1 - w\right)}
1 1 4p0w2(1−w)2+p1w4+2p1w2(1−w)2+p1(1−w)4+4p2w2(1−w)24p02w2(1−w)2+4p0p1w3(1−w)+4p0p1w(1−w)3+8p0p2w2(1−w)2+p12w4+2p12w2(1−w)2+p12(1−w)4+4p1p2w3(1−w)+4p1p2w(1−w)3+4p22w2(1−w)2\frac{4 p_{0} w^{2} \left(1 - w\right)^{2} + p_{1} w^{4} + 2 p_{1} w^{2} \left(1 - w\right)^{2} + p_{1} \left(1 - w\right)^{4} + 4 p_{2} w^{2} \left(1 - w\right)^{2}}{4 p_{0}^{2} w^{2} \left(1 - w\right)^{2} + 4 p_{0} p_{1} w^{3} \left(1 - w\right) + 4 p_{0} p_{1} w \left(1 - w\right)^{3} + 8 p_{0} p_{2} w^{2} \left(1 - w\right)^{2} + p_{1}^{2} w^{4} + 2 p_{1}^{2} w^{2} \left(1 - w\right)^{2} + p_{1}^{2} \left(1 - w\right)^{4} + 4 p_{1} p_{2} w^{3} \left(1 - w\right) + 4 p_{1} p_{2} w \left(1 - w\right)^{3} + 4 p_{2}^{2} w^{2} \left(1 - w\right)^{2}}
1 2 2p0w3(1−w)+p1w3(1−w)+p1w(1−w)3+2p2w(1−w)32p02w3(1−w)+p0p1w4+3p0p1w2(1−w)2+2p0p2w3(1−w)+2p0p2w(1−w)3+p12w3(1−w)+p12w(1−w)3+3p1p2w2(1−w)2+p1p2(1−w)4+2p22w(1−w)3\frac{2 p_{0} w^{3} \left(1 - w\right) + p_{1} w^{3} \left(1 - w\right) + p_{1} w \left(1 - w\right)^{3} + 2 p_{2} w \left(1 - w\right)^{3}}{2 p_{0}^{2} w^{3} \left(1 - w\right) + p_{0} p_{1} w^{4} + 3 p_{0} p_{1} w^{2} \left(1 - w\right)^{2} + 2 p_{0} p_{2} w^{3} \left(1 - w\right) + 2 p_{0} p_{2} w \left(1 - w\right)^{3} + p_{1}^{2} w^{3} \left(1 - w\right) + p_{1}^{2} w \left(1 - w\right)^{3} + 3 p_{1} p_{2} w^{2} \left(1 - w\right)^{2} + p_{1} p_{2} \left(1 - w\right)^{4} + 2 p_{2}^{2} w \left(1 - w\right)^{3}}
2 0 p0w2(1−w)2+p1w2(1−w)2+p2w2(1−w)2p02w2(1−w)2+p0p1w3(1−w)+p0p1w(1−w)3+p0p2w4+p0p2(1−w)4+p12w2(1−w)2+p1p2w3(1−w)+p1p2w(1−w)3+p22w2(1−w)2\frac{p_{0} w^{2} \left(1 - w\right)^{2} + p_{1} w^{2} \left(1 - w\right)^{2} + p_{2} w^{2} \left(1 - w\right)^{2}}{p_{0}^{2} w^{2} \left(1 - w\right)^{2} + p_{0} p_{1} w^{3} \left(1 - w\right) + p_{0} p_{1} w \left(1 - w\right)^{3} + p_{0} p_{2} w^{4} + p_{0} p_{2} \left(1 - w\right)^{4} + p_{1}^{2} w^{2} \left(1 - w\right)^{2} + p_{1} p_{2} w^{3} \left(1 - w\right) + p_{1} p_{2} w \left(1 - w\right)^{3} + p_{2}^{2} w^{2} \left(1 - w\right)^{2}}
2 1 2p0w3(1−w)+p1w3(1−w)+p1w(1−w)3+2p2w(1−w)32p02w3(1−w)+p0p1w4+3p0p1w2(1−w)2+2p0p2w3(1−w)+2p0p2w(1−w)3+p12w3(1−w)+p12w(1−w)3+3p1p2w2(1−w)2+p1p2(1−w)4+2p22w(1−w)3\frac{2 p_{0} w^{3} \left(1 - w\right) + p_{1} w^{3} \left(1 - w\right) + p_{1} w \left(1 - w\right)^{3} + 2 p_{2} w \left(1 - w\right)^{3}}{2 p_{0}^{2} w^{3} \left(1 - w\right) + p_{0} p_{1} w^{4} + 3 p_{0} p_{1} w^{2} \left(1 - w\right)^{2} + 2 p_{0} p_{2} w^{3} \left(1 - w\right) + 2 p_{0} p_{2} w \left(1 - w\right)^{3} + p_{1}^{2} w^{3} \left(1 - w\right) + p_{1}^{2} w \left(1 - w\right)^{3} + 3 p_{1} p_{2} w^{2} \left(1 - w\right)^{2} + p_{1} p_{2} \left(1 - w\right)^{4} + 2 p_{2}^{2} w \left(1 - w\right)^{3}}
2 2 p0w4+p1w2(1−w)2+p2(1−w)4p02w4+2p0p1w3(1−w)+2p0p2w2(1−w)2+p12w2(1−w)2+2p1p2w(1−w)3+p22(1−w)4\frac{p_{0} w^{4} + p_{1} w^{2} \left(1 - w\right)^{2} + p_{2} \left(1 - w\right)^{4}}{p_{0}^{2} w^{4} + 2 p_{0} p_{1} w^{3} \left(1 - w\right) + 2 p_{0} p_{2} w^{2} \left(1 - w\right)^{2} + p_{1}^{2} w^{2} \left(1 - w\right)^{2} + 2 p_{1} p_{2} w \left(1 - w\right)^{3} + p_{2}^{2} \left(1 - w\right)^{4}}

Sample-specific genotyping error probabilities wtw_t and wrw_r for trace sample (xt=x_t =xT) and reference sample (xr=x_r =xR), repectively:

tab <- wgsLR::d_LR_wTwR[, c("xT", "xR")]
tab$LR <- paste0("$", wgsLR::d_LR_wTwR$expr_tex, "$")
tab |> 
  kbl(format = "markdown")
xT xR LR
0 0 p0(wr−1)2(wt−1)2+p1wrwt(wr−1)(wt−1)+p2wr2wt2p02(wr−1)2(wt−1)2−p0p1wr(wr−1)(wt−1)2−p0p1wt(wr−1)2(wt−1)+p0p2wr2(wt−1)2+p0p2wt2(wr−1)2+p12wrwt(wr−1)(wt−1)−p1p2wr2wt(wt−1)−p1p2wrwt2(wr−1)+p22wr2wt2\frac{p_{0} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + p_{2} w_r^{2} w_t^{2}}{p_{0}^{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} - p_{0} p_{1} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} - p_{0} p_{1} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{0} p_{2} w_r^{2} \left(w_t - 1\right)^{2} + p_{0} p_{2} w_t^{2} \left(w_r - 1\right)^{2} + p_{1}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_r^{2} w_t \left(w_t - 1\right) - p_{1} p_{2} w_r w_t^{2} \left(w_r - 1\right) + p_{2}^{2} w_r^{2} w_t^{2}}
0 1 2p0wr(wr−1)(wt−1)2+p1wr2wt(wt−1)+p1wt(wr−1)2(wt−1)+2p2wrwt2(wr−1)2p02wr(wr−1)(wt−1)2−p0p1wr2(wt−1)2−2p0p1wrwt(wr−1)(wt−1)−p0p1(wr−1)2(wt−1)2+2p0p2wrwt2(wr−1)+2p0p2wr(wr−1)(wt−1)2+p12wr2wt(wt−1)+p12wt(wr−1)2(wt−1)−p1p2wr2wt2−2p1p2wrwt(wr−1)(wt−1)−p1p2wt2(wr−1)2+2p22wrwt2(wr−1)\frac{2 p_{0} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + p_{1} w_r^{2} w_t \left(w_t - 1\right) + p_{1} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + 2 p_{2} w_r w_t^{2} \left(w_r - 1\right)}{2 p_{0}^{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} - p_{0} p_{1} w_r^{2} \left(w_t - 1\right)^{2} - 2 p_{0} p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{0} p_{1} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + 2 p_{0} p_{2} w_r w_t^{2} \left(w_r - 1\right) + 2 p_{0} p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + p_{1}^{2} w_r^{2} w_t \left(w_t - 1\right) + p_{1}^{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) - p_{1} p_{2} w_r^{2} w_t^{2} - 2 p_{1} p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_t^{2} \left(w_r - 1\right)^{2} + 2 p_{2}^{2} w_r w_t^{2} \left(w_r - 1\right)}
0 2 p0wr2(wt−1)2+p1wrwt(wr−1)(wt−1)+p2wt2(wr−1)2p02wr2(wt−1)2−p0p1wr2wt(wt−1)−p0p1wr(wr−1)(wt−1)2+p0p2wr2wt2+p0p2(wr−1)2(wt−1)2+p12wrwt(wr−1)(wt−1)−p1p2wrwt2(wr−1)−p1p2wt(wr−1)2(wt−1)+p22wt2(wr−1)2\frac{p_{0} w_r^{2} \left(w_t - 1\right)^{2} + p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + p_{2} w_t^{2} \left(w_r - 1\right)^{2}}{p_{0}^{2} w_r^{2} \left(w_t - 1\right)^{2} - p_{0} p_{1} w_r^{2} w_t \left(w_t - 1\right) - p_{0} p_{1} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + p_{0} p_{2} w_r^{2} w_t^{2} + p_{0} p_{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + p_{1}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_r w_t^{2} \left(w_r - 1\right) - p_{1} p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{2}^{2} w_t^{2} \left(w_r - 1\right)^{2}}
1 0 2p0wt(wr−1)2(wt−1)+p1wrwt2(wr−1)+p1wr(wr−1)(wt−1)2+2p2wr2wt(wt−1)2p02wt(wr−1)2(wt−1)−2p0p1wrwt(wr−1)(wt−1)−p0p1wt2(wr−1)2−p0p1(wr−1)2(wt−1)2+2p0p2wr2wt(wt−1)+2p0p2wt(wr−1)2(wt−1)+p12wrwt2(wr−1)+p12wr(wr−1)(wt−1)2−p1p2wr2wt2−p1p2wr2(wt−1)2−2p1p2wrwt(wr−1)(wt−1)+2p22wr2wt(wt−1)\frac{2 p_{0} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{1} w_r w_t^{2} \left(w_r - 1\right) + p_{1} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + 2 p_{2} w_r^{2} w_t \left(w_t - 1\right)}{2 p_{0}^{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) - 2 p_{0} p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{0} p_{1} w_t^{2} \left(w_r - 1\right)^{2} - p_{0} p_{1} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + 2 p_{0} p_{2} w_r^{2} w_t \left(w_t - 1\right) + 2 p_{0} p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{1}^{2} w_r w_t^{2} \left(w_r - 1\right) + p_{1}^{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} - p_{1} p_{2} w_r^{2} w_t^{2} - p_{1} p_{2} w_r^{2} \left(w_t - 1\right)^{2} - 2 p_{1} p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + 2 p_{2}^{2} w_r^{2} w_t \left(w_t - 1\right)}
1 1 −4p0wrwt(wr−1)(wt−1)−p1wr2wt2−p1wr2(wt−1)2−p1wt2(wr−1)2−p1(wr−1)2(wt−1)2−4p2wrwt(wr−1)(wt−1)−4p02wrwt(wr−1)(wt−1)+2p0p1wr2wt(wt−1)+2p0p1wrwt2(wr−1)+2p0p1wr(wr−1)(wt−1)2+2p0p1wt(wr−1)2(wt−1)−8p0p2wrwt(wr−1)(wt−1)−p12wr2wt2−p12wr2(wt−1)2−p12wt2(wr−1)2−p12(wr−1)2(wt−1)2+2p1p2wr2wt(wt−1)+2p1p2wrwt2(wr−1)+2p1p2wr(wr−1)(wt−1)2+2p1p2wt(wr−1)2(wt−1)−4p22wrwt(wr−1)(wt−1)\frac{- 4 p_{0} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} w_r^{2} w_t^{2} - p_{1} w_r^{2} \left(w_t - 1\right)^{2} - p_{1} w_t^{2} \left(w_r - 1\right)^{2} - p_{1} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} - 4 p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right)}{- 4 p_{0}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + 2 p_{0} p_{1} w_r^{2} w_t \left(w_t - 1\right) + 2 p_{0} p_{1} w_r w_t^{2} \left(w_r - 1\right) + 2 p_{0} p_{1} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + 2 p_{0} p_{1} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) - 8 p_{0} p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1}^{2} w_r^{2} w_t^{2} - p_{1}^{2} w_r^{2} \left(w_t - 1\right)^{2} - p_{1}^{2} w_t^{2} \left(w_r - 1\right)^{2} - p_{1}^{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + 2 p_{1} p_{2} w_r^{2} w_t \left(w_t - 1\right) + 2 p_{1} p_{2} w_r w_t^{2} \left(w_r - 1\right) + 2 p_{1} p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + 2 p_{1} p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) - 4 p_{2}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right)}
1 2 2p0wr2wt(wt−1)+p1wrwt2(wr−1)+p1wr(wr−1)(wt−1)2+2p2wt(wr−1)2(wt−1)2p02wr2wt(wt−1)−p0p1wr2wt2−p0p1wr2(wt−1)2−2p0p1wrwt(wr−1)(wt−1)+2p0p2wr2wt(wt−1)+2p0p2wt(wr−1)2(wt−1)+p12wrwt2(wr−1)+p12wr(wr−1)(wt−1)2−2p1p2wrwt(wr−1)(wt−1)−p1p2wt2(wr−1)2−p1p2(wr−1)2(wt−1)2+2p22wt(wr−1)2(wt−1)\frac{2 p_{0} w_r^{2} w_t \left(w_t - 1\right) + p_{1} w_r w_t^{2} \left(w_r - 1\right) + p_{1} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + 2 p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right)}{2 p_{0}^{2} w_r^{2} w_t \left(w_t - 1\right) - p_{0} p_{1} w_r^{2} w_t^{2} - p_{0} p_{1} w_r^{2} \left(w_t - 1\right)^{2} - 2 p_{0} p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + 2 p_{0} p_{2} w_r^{2} w_t \left(w_t - 1\right) + 2 p_{0} p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{1}^{2} w_r w_t^{2} \left(w_r - 1\right) + p_{1}^{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} - 2 p_{1} p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_t^{2} \left(w_r - 1\right)^{2} - p_{1} p_{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + 2 p_{2}^{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right)}
2 0 p0wt2(wr−1)2+p1wrwt(wr−1)(wt−1)+p2wr2(wt−1)2p02wt2(wr−1)2−p0p1wrwt2(wr−1)−p0p1wt(wr−1)2(wt−1)+p0p2wr2wt2+p0p2(wr−1)2(wt−1)2+p12wrwt(wr−1)(wt−1)−p1p2wr2wt(wt−1)−p1p2wr(wr−1)(wt−1)2+p22wr2(wt−1)2\frac{p_{0} w_t^{2} \left(w_r - 1\right)^{2} + p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + p_{2} w_r^{2} \left(w_t - 1\right)^{2}}{p_{0}^{2} w_t^{2} \left(w_r - 1\right)^{2} - p_{0} p_{1} w_r w_t^{2} \left(w_r - 1\right) - p_{0} p_{1} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{0} p_{2} w_r^{2} w_t^{2} + p_{0} p_{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + p_{1}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_r^{2} w_t \left(w_t - 1\right) - p_{1} p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + p_{2}^{2} w_r^{2} \left(w_t - 1\right)^{2}}
2 1 2p0wrwt2(wr−1)+p1wr2wt(wt−1)+p1wt(wr−1)2(wt−1)+2p2wr(wr−1)(wt−1)22p02wrwt2(wr−1)−p0p1wr2wt2−2p0p1wrwt(wr−1)(wt−1)−p0p1wt2(wr−1)2+2p0p2wrwt2(wr−1)+2p0p2wr(wr−1)(wt−1)2+p12wr2wt(wt−1)+p12wt(wr−1)2(wt−1)−p1p2wr2(wt−1)2−2p1p2wrwt(wr−1)(wt−1)−p1p2(wr−1)2(wt−1)2+2p22wr(wr−1)(wt−1)2\frac{2 p_{0} w_r w_t^{2} \left(w_r - 1\right) + p_{1} w_r^{2} w_t \left(w_t - 1\right) + p_{1} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + 2 p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2}}{2 p_{0}^{2} w_r w_t^{2} \left(w_r - 1\right) - p_{0} p_{1} w_r^{2} w_t^{2} - 2 p_{0} p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{0} p_{1} w_t^{2} \left(w_r - 1\right)^{2} + 2 p_{0} p_{2} w_r w_t^{2} \left(w_r - 1\right) + 2 p_{0} p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} + p_{1}^{2} w_r^{2} w_t \left(w_t - 1\right) + p_{1}^{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) - p_{1} p_{2} w_r^{2} \left(w_t - 1\right)^{2} - 2 p_{1} p_{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2} + 2 p_{2}^{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2}}
2 2 p0wr2wt2+p1wrwt(wr−1)(wt−1)+p2(wr−1)2(wt−1)2p02wr2wt2−p0p1wr2wt(wt−1)−p0p1wrwt2(wr−1)+p0p2wr2(wt−1)2+p0p2wt2(wr−1)2+p12wrwt(wr−1)(wt−1)−p1p2wr(wr−1)(wt−1)2−p1p2wt(wr−1)2(wt−1)+p22(wr−1)2(wt−1)2\frac{p_{0} w_r^{2} w_t^{2} + p_{1} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) + p_{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2}}{p_{0}^{2} w_r^{2} w_t^{2} - p_{0} p_{1} w_r^{2} w_t \left(w_t - 1\right) - p_{0} p_{1} w_r w_t^{2} \left(w_r - 1\right) + p_{0} p_{2} w_r^{2} \left(w_t - 1\right)^{2} + p_{0} p_{2} w_t^{2} \left(w_r - 1\right)^{2} + p_{1}^{2} w_r w_t \left(w_r - 1\right) \left(w_t - 1\right) - p_{1} p_{2} w_r \left(w_r - 1\right) \left(w_t - 1\right)^{2} - p_{1} p_{2} w_t \left(w_r - 1\right)^{2} \left(w_t - 1\right) + p_{2}^{2} \left(w_r - 1\right)^{2} \left(w_t - 1\right)^{2}}

Example

Assume that a trace sample had four loci with genotypes (0/1 = 1, 0/0 = 0, 1/1 = 2, 1/1 = 2). A person of interest is then typed for the same four loci and has genotypes (0/0 = 0, 0/0 = 0, 1/1 = 2, 1/1 = 2), i.e. a mismatch on the first locus.

For simplicity, assume that the genotype probabilites are P(0/0 = 0) = 0.25, P(0/1 = 1/0 = 1) = 0.25, and P(1/1 = 2) = 0.5.

If no errors are possible, then w=0w=0 and

wgsLR::calc_LRs_w(xT = c(1, 0, 2, 2),
                  xR = c(0, 0, 2, 2), 
                  w = 0, 
                  p = c(0.25, 0.25, 0.5))
## [1] 0 4 2 2

and the product is 0 due to the mismatch at the first locus.

If instead we acknowledge that errors are possible, then for w=0.001w = 0.001 we obtain that

LR_contribs <- wgsLR::calc_LRs_w(xT = c(1, 0, 2, 2), 
                                 xR = c(0, 0, 2, 2), 
                                 w = 0.001, 
                                 p = c(0.25, 0.25, 0.5))
LR_contribs
## [1] 0.01192834 3.99199202 1.99799850 1.99799850
prod(LR_contribs)
## [1] 0.1900904

We can also consider the LRLRs for a range for plausible values of ww:

ws <- c(1e-6, 1e-3, 1e-2, 1e-1)
LRs <- sapply(ws, \(w) wgsLR::calc_LRs_w(xT = c(0, 0, 2, 2), 
                                         xR = c(1, 0, 2, 2), 
                                         w = w, 
                                         p = c(0.25, 0.25, 0.5)) |> 
                prod())
data.frame(log10w = log10(ws), w = ws, 
           LR = LRs, WoElog10LR = log10(LRs))
##   log10w     w           LR WoElog10LR
## 1     -6 1e-06 0.0001919981 -3.7167031
## 2     -3 1e-03 0.1900903523 -0.7210399
## 3     -2 1e-02 1.7379421372  0.2400353
## 4     -1 1e-01 7.1409934157  0.8537586

Different error rates

Assume that the trace donor profile has wD=10−4w_D = 10^{-4} and the suspect reference profile has wS=10−8w_S = 10^{-8}. Then the LRLR is:

LR_contribs <- wgsLR::calc_LRs_wTwR(xT = c(1, 0, 2, 2), 
                                    xR = c(0, 0, 2, 2), 
                                    wT = 1e-4, 
                                    wR = 1e-8,
                                    p = c(0.25, 0.25, 0.5))
prod(LR_contribs)
## [1] 0.01279168