On the exact distribution of the numbers of alleles in DNA mixtures

In DNAtools: Tools for Analysing Forensic Genetic DNA Data

library(knitr)
opts_chunk$set(
  collapse = TRUE,
  comment = ""
)
options(knitr.kable.NA = '')

Introduction

Please refer to the dna_vignette ("STR markers in forensic genetics") for an explanation of forensic DNA profiles.

When more than one individual contributes biological material to a forensic stain, the resulting DNA type is termed a DNA mixture. DNA mixtures occur frequently in forensic genetic casework, and in recent years much research has been devoted to this subject.

This vignette shows how to compute the exact distribution of the number of alleles for any number of profiles and investigated loci using the DNAtools package in R.

library(DNAtools)

The per locus number of observed alleles is of interest as it indicates the plausible range of the number of contributors. If a DNA mixture has $m$ contributors, it is possible to observe anything between $1$ and $2m$ alleles at a given marker. Of particular interest is the probability at which an $m+1$ person DNA mixture is percieved as an $m$ person DNA mixture. That is, what is the probability that $m$ persons have at most $2(m-1)$ distinct alleles among them?

DNA mixtures

In forensic genetics, a DNA mixture refers to the situation where more than one individual's DNA is present in an analysed trace. DNA mixtures are regularly encountered in forensic genetics.

However, an exact expression of the number of alleles has not been presented up to now. We show how to derive the probability that $N$ alleles are observed given $m$ contributors and $L$ typed loci.

The prevailing technology for identity assessment and DNA mixture analysis in forensic genetics uses so-called Short Tandem Repeat (STR) markers or loci. The number of alleles observed at these markers are typically in the range of 10 to 30 alleles.

In the DNAtools there is a database with 1,000 simulated DNA profiles on 10 autosomal STR markers, which we will use to compute the allele frequencies.

data(dbExample, package = "DNAtools")
kable(head(dbExample)[,2:7])

We extract the relevant columns and compute the allele frequencies by the allele count and total number of alleles in the data.

allele_freqs <- lapply(1:10, function(x){
  al_freq <- table(c(dbExample[[x*2]], dbExample[[1+x*2]]))/(2*nrow(dbExample))
  al_freq[sort.list(as.numeric(names(al_freq)))]
})
names(allele_freqs) <- sub("\\.1", "", names(dbExample)[(1:10)*2])

for(l in head(names(allele_freqs))){
  # cat(paste0("**",l,"**\n\n"))
  cat("\n")
  cat(kable(t(c(setNames(NA, paste0(l,":")), allele_freqs[[l]]))), sep = "\n")
}

Single locus

The total number of alleles observed for $m$ contributors depends on the number of loci, $L$, through the locus specific allele counts, $N_l$, by $N(m) = \sum_{l=1}^L N_l(m)$. Hence, we can focus on computing the distribution of $N_l(m)$ below.

The number of alleles at locus $l$, $N_l(m)$, follows a distribution, where the number of alleles range from 1 (all $m$ profiles homozygous with the same allele) to $2m$ (all $m$ profiles heterozygous sharing no allele). For each value that $N_l(m)$ can attain there will typically be several ways $m$ profiles can produce $N_l(m)$ distinct alleles. For a given locus $l$, we use the function Pnm_locus() to compute $P(N_l(m) = n)$, for $n = 1,\dots,2m$, where we need to specify $m$ together with the allele frequencies.

m <- 3
P3_D16 <- Pnm_locus(m = m, theta = 0, alleleProbs = allele_freqs$D16S539)

cat(kable(t(setNames(P3_D16, paste0("P(n=",1:(2*m), "|m=3)")))), sep = "\n")

The theta argument refers to the $\theta$ coefficient used to account for genetic dependence of alleles in the Balding-Nicholls population genetic model. Typical values ranges from $0$ (independent alleles, no subpopulation structure) to $0.03$ (mild correlation, some subpopulation structure).

We can also compute the probabilities for several $m$. Here we let $m$ range between 1 and 6 contributors.

ms <- 1:6
Ps_D16 <- matrix(NA, ncol = length(ms), nrow = max(ms)*2, 
                 dimnames = list(alleles = 1:(max(ms)*2), noContrib = ms))
for(m in seq_along(ms)){
  Ps_D16[1:(ms[m]*2), m] <- Pnm_locus(m = ms[m], theta = 0, alleleProbs = allele_freqs$D16S539)
}

dimnames(Ps_D16) <- list(paste0("n=", rownames(Ps_D16),""), 
                         paste0("P(n|m=", colnames(Ps_D16),")"))
kable(Ps_D16, row.names = TRUE)

If we plot the outcome we see that for $m \in {4, 5, 6}$ the most probable allele count is 4 indicating that D16S539 is not very polymorphic in this population.

par(mar = c(4,5,0,0))
plot(c(1,max(ms)*2), range(Ps_D16, na.rm = TRUE), type = "n",
     xlab = "no. of allele", ylab = "Probability")
for (m in seq_along(ms)) {
  lines(1:(2*ms[m]), Ps_D16[1:(ms[m]*2),m], col = m)
}
legend("topright", bty = "n", title = "m", legend = ms, col = seq_along(ms), lty = 1)

The probability of masking can be calculated from the table above. Hence, we need to compute $P(N_l(m) \le 2(m-1))$, which is done below.

Ps_D16_cumsum <- apply(Ps_D16, 2, cumsum)

nm <- dim(Ps_D16_cumsum)
dimnames(Ps_D16_cumsum) <- list(paste0("n'=",1:nm[1]), 
                                paste0("P(n&le;n'|m=",1:nm[2],")"))
kable(Ps_D16_cumsum, row.names = TRUE, escape = TRUE)

We see that the probability that an $m$ person mixture is misconceived as an $m-1$ person mixture is

Pm_D16 <- setNames(rep(NA, max(ms)-1), paste0("P(n&le;",(1:(max(ms)-1))*2,"|m=",2:max(ms),")"))
for(m in 2:max(ms)) Pm_D16[m-1] <- Ps_D16_cumsum[2*(m-1),m]
kable(t(Pm_D16))

Hence, there is a probability of r Pm_D16[2] that a three-person mixture has at most four alleles, which implies it can be interpreted as a two-person mixture.

If we are interested in computing the locuswise probabilities for all loci in our dataset, we simply use the Pnm_all function with locuswise = TRUE.

no_contrib <- 3
Pnm_all_tab <- Pnm_all(m = no_contrib, theta = 0, probs = allele_freqs, locuswise = TRUE)

kable(Pnm_all_tab, row.names = TRUE, col.names = paste0("P(n=",1:(2*no_contrib),"|m=",no_contrib,")"))

Accumulating over loci

Since the loci are assumed independent, we can also compute the probability that we see at most $n_0$ alleles at all markers: $$ P(N_1(m) \le n_0, \dots, N_L(m)\le n_0) = \prod_{l=1}^L P(N_l(m) \le n_0) $$

For instance, we may compute for $m\in{1,\dots,6}$ the probability of seeing at most $n_0 \in{1,\dots,2m}$ alleles for the ten STR markers in allele_freqs.

locus_Ps <- lapply(ms, Pnm_all, theta = 0, probs = allele_freqs, locuswise = TRUE)
locus_Ps_cumsum <- lapply(locus_Ps, apply, 1, cumsum)
all_Ps_cumsum <- lapply(locus_Ps_cumsum, apply, 1, prod)
all_Ps_table <- matrix(NA, ncol = length(ms), nrow = max(ms)*2, 
                       dimnames = list(n0 = 1:(max(ms)*2), m = ms))
for(m in seq_along(ms)) all_Ps_table[1:(ms[m]*2), m] <- all_Ps_cumsum[[m]]

nm <- dim(all_Ps_table)
dimnames(all_Ps_table) <- list(paste0("n~0~=",1:nm[1]), paste0("P(n&le;n~0~|m=",1:nm[2],")"))
kable(all_Ps_table, row.names = TRUE, escape = FALSE)

Again we may investigate the risk the a $m$-person DNA mixture is misconceived as a $(m-1)$-person DNA mixture, when considering all loci:

Pm_all <- setNames(rep(NA, max(ms)-1), paste0("P(n&le;",(1:(max(ms)-1))*2,"|m=",2:max(ms),")"))
for(m in 2:max(ms)) Pm_all[m-1] <- all_Ps_table[2*(m-1),m]

kable(t(round(Pm_all, 4)))

$\theta$-correction

We can also look at the effect of the $\theta$-correction on the distribution of observed number of alleles. The stronger the subpopulation effect (measured by increasing $\theta$), the fewer alleles will be observed.

par(mar = c(4,5,0,0))
P3_D16_t0.03 <- Pnm_locus(m = 3, theta = 0.03, alleleProbs = allele_freqs$D16S539)
P3_D16_t0.1 <- Pnm_locus(m = 3, theta = 0.1, alleleProbs = allele_freqs$D16S539)
plot(P3_D16_t0.03, type = "o", pch = 16, ylab = "Probability", xlab = "no. of alleles", col = 2)
points(P3_D16, type = "o", lty = 2, col = 1)
points(P3_D16_t0.1, type = "o", lty = 3, pch = 2, col = 3)
legend("topleft", bty = "n", title = expression(theta), col = 1:3,
       legend = c("0.00", "0.03", "0.10"), lty = c(2, 1, 3), pch = c(1, 16, 2))

Summing over loci

In order to obtain $n$ alleles for $m$ contributors and $L$ loci, the locus counts must sum to $n$. By convolution, we obtain $$ P(N^{l{+}1}(m) = n) = \sum_{i=1}^I P(N^{l}(m) = n-i)P(N_{l+1}(m)=i), $$ where $I = \min(2m,n{-}1)$ with $P(N_{l+1}(m)=i)$ denoting the probability that locus $l{+}1$ has $i$ observed alleles, and $P(N^l(m)=n)$ is the probability to observe $n$ alleles for $l$ loci and $m$ profiles.

We compute this by setting locuswise = FALSE in Pnm_all.

par(mar = c(4,5,0,0))
P3_all <- Pnm_all(m = 3, theta = 0, probs = allele_freqs, locuswise = FALSE)
plot(P3_all, xlab = "no. of alleles", ylab = "Probability", type = "o")

Posterior probability for the number of contributors

Assume that we specify the prior, $P(m)$, on $m$ based on some general idea on the number of contributors we typically see, or based on other background information regarding the case. Then $P(m)$ reflects the belief we have about $m$ prior to seeing the DNA mixture.

Using Bayes' theorem, we can express the probability distribution for the number of contributors, $m$, given the number of observed alleles, $n_{l}$: $$ P(m\mid n_{l}) = \frac{P(n_{l}\mid m)P(m)}{P(n_{l})} = \frac{P(N_l(m) = n_{l})P(m)}{\sum_m P(N_l(m) = n_{l})P(m)}, $$ where $P(m)$ is a prior distribution on the number of contributors, $m$. We may also generalise this to including all information, $\underline{n} = (n_1, \dots, n_L)$, from the loci based on the independent assumption.

$$ P(m\mid \underline{n}) = \frac{P(\underline{n}\mid m)P(m)}{P(\underline{n})} = \frac{P(m)\prod_{l=1}^L P(N_l(m) = n_{l})}{\sum_m\left{ P(m)\prod_{l=1}^LP(N_l(m) = n_{l})\right}}, $$

Assume that $P(m)$ is given by

par(mar = c(4,5,0,0))
m_prior <- c(0.2,0.35,0.3,0.2,0.1,0.025)
m_prior <- m_prior/sum(m_prior)
barplot(m_prior, names.arg = seq_along(m_prior), ylab = "Prior probability,"~italic(P(m)), xlab = expression(italic(m)))

Let $n_{vWA} = 4$. Then we can update the belief based on this information using the expression above.

n_vWA <- 4
P_vWA <- lapply(ms, Pnm_locus, theta = 0, alleleProbs = allele_freqs$vWA)
P_vWA_n <- sapply(P_vWA, function(p) ifelse(length(p) < n_vWA, 0, p[n_vWA]))
Pn_vWA <- (m_prior * P_vWA_n)/sum(m_prior * P_vWA_n)

Hence, the change in belief about $m$ can be seen below, where the first row is the prior $P(m)$, the second $P(N_{vWA}(m) = n_{vWA})$ and the last $P(n \mid N_{vWA}(m) = n_{vWA})$:

kable(rbind("P(m)" = m_prior, "P(n~vWA~=4|m)" = P_vWA_n, "P(m|n~vWA~=4)" = Pn_vWA), col.names = paste0("m=",paste(ms)))

We see that $m =$ r which.max(Pn_vWA) is the most probable after seeing r n_vWA alleles at vWA.

Assuming that we observe $n_l = 4$ for all loci, $l=1,\dots,L$, we get

get_pn <- function(p, n){
  pn <- split(as.data.frame(p), n)
  p_n <- unlist(lapply(names(pn), function(n_l) if(ncol(pn[[n_l]])<n_l) rep(0,nrow(pn[[n_l]])) else pn[[n_l]][,n_l]))
  prod(p_n)
}

n_loci <- rep(4, length(allele_freqs))
P_all <- lapply(ms, Pnm_all, theta = 0, probs = allele_freqs, locuswise = TRUE)
P_all_n <- sapply(P_all, get_pn, n_loci)
Pn_all <- (m_prior * P_all_n)/sum(m_prior * P_all_n)

Hence, the posterior probability of $m$ after observing four alleles at all loci is given by

kable(t(setNames(Pn_all, paste0("P(m=", seq_along(Pn_all),"|n=4)"))))

DNAtools documentation built on March 18, 2022, 7:01 p.m.