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):
You can also install the package without vignettes if needed as follows:
A few small examples
Estimating the genotype error probability,
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
w_mle <- wgsLR::estimate_w(tab)
w_mle
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 (’s)
Same genotyping error probability for both trace sample (xT
) and reference sample (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 | |
0 | 1 | |
0 | 2 | |
1 | 0 | |
1 | 1 | |
1 | 2 | |
2 | 0 | |
2 | 1 | |
2 | 2 |
Sample-specific genotyping error probabilities and for trace sample (xT
) and reference sample (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 | |
0 | 1 | |
0 | 2 | |
1 | 0 | |
1 | 1 | |
1 | 2 | |
2 | 0 | |
2 | 1 | |
2 | 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 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))
and the product is 0 due to the mismatch at the first locus.
If instead we acknowledge that errors are possible, then for 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
prod(LR_contribs)
We can also consider the s for a range for plausible values of :
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))
Different error rates
Assume that the trace donor profile has and the suspect reference profile has . Then the 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)