magnusdv / forrel

Forensic pedigree analysis and relatedness inference
GNU General Public License v2.0
10 stars 0 forks source link

'ids' argument in showInTriangle #41

Closed hildekb closed 3 years ago

hildekb commented 3 years ago

Assume that IBD coefficients have been estimated between all pedigree members. The output from e.g. IBDestimate() can then be used directly as input to showInTriangle() and the set of coefficients for each pair are plotted.

My question is then: Is it possible to specify a subset of pedigree member pairs to plot when using showInTriangle()? I cannot find 'ids' as an option in showInTriangle. This feature would be helpful.

x = nuclearPed(father = "fa",mother = "mo", children = c("ch1","ch2"), sex = c(1,2)) mlist = lapply(1:4, function(i) marker(x, alleles = letters[1:10], name = paste0("M",i))) x = setMarkers(x, mlist) x = profileSim(x,N = 1) k = IBDestimate(x) IBDtriangle()

Suggestion of code:

Only want to plot coefficients between ch1 and ch2: showInTriangle(k, ids = c("ch1","ch2"))

Or if more pairs are wanted, a matrix of ids may be specified: showInTriangle(k, ids = rbind(c("ch1","ch2"),c("fa","mo")))

magnusdv commented 3 years ago

Hi, that's a good question. There are currently two ways of doing this, but neither is very well explained in the documentation:

1) subsetting the output from IBDestimate(). 1) restrict to the subset of individuals already when estimating.

The following code shows both of these. Note that I used markerSim() to simulate: It is much simpler than profileSim() in the special case when you want a single profile with several identical markers.

library(forrel)
#> Loading required package: pedtools
x = nuclearPed(fa = "fa", mo = "mo", ch = c("ch1","ch2"), sex = 1:2)

# Simulate a profile with 4 iid markers 
x = markerSim(x, N = 4, alleles = letters[1:10], seed = 123, verb = F)

# 1) Estimate kappa for all pairs, but plot sibs only
k = IBDestimate(x)
showInTriangle(k[6, ])   # select the appropriate row(s)

# Same, but more robust: use `subset()` on `k`:
kSibs = subset(k, id1 == "ch1" & id2 == "ch2")
showInTriangle(kSibs)

####

# 2) Estimate sibs only:
ks = IBDestimate(x, ids = c("ch1","ch2"))
showInTriangle(ks)

Created on 2020-10-19 by the reprex package (v0.3.0)

PS: Check out the reprex package! :)