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)

We see that the probability that an \(m\) person mixture is misconceived as an \(m-1\) person mixture is

Hence, there is a probability of 0.9723436 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)

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]]

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]

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

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})\):

We see that \(m =\) 3 is the most probable after seeing 4 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