EP: mutation only possible toobserved alleles

[thoree]

The code below returns EP = 0; I thought EP > 0

library(forrel)
#> Loading required package: pedtools
z = nuclearPed(4, children = c(3, 4, 5, "POI"))
p = c(0.1, 0.2, 0.3, 0.4)
als = 1:4
m = marker(z, alleles = als, afreq = p, "3" = 1, "4" = 2, "5" = c(1,3) )
z = addMarkers(z,m)
#plot(z, 1, skip.empty.genotypes = TRUE)
r = 0.06
M = rbind(c(1-r, r/2, r/2, 0),
          c(r/2, 1-r, r/2, 0),
          c(r/2, r/2, 1-r, 0),
          c(0, 0,  0, 1))
dimnames(M) = list(1:4, 1:4)
mutmod(z,1) = list("custom", matrix = M)
missingPersonEP(z, "POI", disableMutations = FALSE, verbose = FALSE)
#> $EPperMarker
#> 1 
#> 0 
#> 
#> $EPtotal
#> [1] 0
#> 
#> $expectedMismatch
#> [1] 0
#> 
#> $distribMismatch
#> 0 
#> 1
# Thore thinks EP should be:
p[4]^2+2*p[4]*(1-p[4])
#> [1] 0.64

^{Created on 2019-10-30 by the reprex package (v0.3.0)}

[magnusdv]

I think EP = 0 is correct in this case. Given your mutation matrix it is possible that both parents carry the 4 allele; for example they could both have genotype 1/4. Therefore POI = 4/4 does not give exclusion.

[thoree]

You are right,thanks. By changing M[4,4]=0 and M[4,3]=1, EP =0.64 as it should.