Unable to use know genotype data in pedigrees

[knifecake]

I'm trying to get exclusionPower to run with known genotypes passed within the pedigrees. The following example fails, I think because the unrelated pedigree does not have marker data associated with it.

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))

m = marker(claim,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = c(0.5, 0.5))

claim = setMarkers(claim, m)

exclusionPower(claim, true, ids = c(1, 3), markerindex = 1, plot = FALSE)
#> Error: `partialmarker` must have length 1

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

On the contrary, if I try the same example but directly supplying the genotypes as arguments to the exclusionPower function, I get the correct result:

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))

exclusionPower(claim, true,
               ids = c(1, 3),
               known_genotypes = list(c(3, 1, 1)),
               alleles = 2,
               afreq = c(0.5, 0.5),
               plot = FALSE)
#> [1] 0.25

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

[thoree]

I run into the same problem. The below mathematically equivalent, but artifical, reformulation, works, indicating that the problem relates to singletons. Also both pedigrees need to have the marker added contradicting the help file. exclusionPower, does not mind if sex of '3' differs for claim and true (OK?):

library(forrel, quietly = TRUE)
claim = nuclearPed(1)
true = halfSibPed(1)

m = marker(claim,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = c(0.5, 0.5))
claim = setMarkers(claim, m)
m = marker(true,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = c(0.5, 0.5))
true = setMarkers(true, m)
exclusionPower(claim, true, ids = c(1, 3), markerindex = 1, plot = TRUE)

#> [1] 0.25

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

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

[knifecake]

Thanks for checking @thoree!

On the latest commit (51b9a824eb584d06382e1494c9f3f66112455ae7) I implemented known genotypes through arguments to the exclusionPower function because I could not get it to work with genotypes attached to pedigrees. However, I think that implementing genotypes inside pedigrees is non-avoidable given the current definition of the exclusionPower function if we want to support mutation models.

As this only affects the unrelated pedigree I will implement genotype attaching to pedigrees and we can think of what to do for the special unrelated case. What you propose is a viable workaround actually.

[magnusdv]

Hi, thanks for experimenting with this.

The original error is clearly a bug. So I prefer fixing it rather than going for acrobatic reformulations of the problem. :)

I have just commited (35264d2) a complete rewrite of how known genotypes are handled by exclusionPower(). The most important change is that when markerindex is given (an integer), the indicated marker is taken from claim and transferred to true using transferMarkers().

Any marker data attached to true is ignored. Other than that, the changes I've done are all internal, and shouldn't (in theory) break any existing code.

Your example should now run without problems:

library(forrel)
#> Loading required package: pedtools
claim = nuclearPed(1)
true = list(singleton(1), singleton(3))

m = marker(claim,
           `3` = c(1, 1),
           alleles = c(1, 2),
           afreq = c(0.5, 0.5))

claim = setMarkers(claim, m)

exclusionPower(claim, true, ids = c(1, 3), markerindex = 1, plot = FALSE)
#> Removed already typed individuals from internal computations: 3
#> [1] 0.25

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

There are definitely edges to be smoothed out in the new code. In particular, more tests are needed. @thoree if you have some simple exclusion power cases (paternity, sibs, etc) with known genotypes lying around (with "formula answers"), they would be much appreciated!

One last thing: I commited this to my "exclusion-power" branch, which is already a few commits ahead of master. Is this ok for you, @knifecake? Better alternatives?