Nothing
#' Pedigree utilities
#'
#' Various utility functions for `ped` objects.
#'
#' @param x A `ped` object, or (in some functions - see Details) a list of such.
#' @param what Either "max", "compMax", "indiv" or "depth" (See Value.)
#' @param chromType Either "autosomal" (default) or "x".
#'
#' @return
#'
#' * `pedsize(x)` returns the number of pedigree members in each component of
#' `x`.
#'
#' * `generations(x)` by default returns the number of generations in `x`,
#' defined as the number of individuals in the longest line of parent-child
#' links. (Note that this is well-defined also if `x` has loops and/or
#' cross-generational marriages.) For individual generation numbers, use `what =
#' "indiv"` (generation numbering as in the plot) or `what = "depth" (length of
#' the longest chain up to a founder). Finally, if `x` has multiple components,
#' and what = "compMax"`, the function returns a vector with the generation
#' count from each component.
#'
#' * `hasUnbrokenLoops(x)` returns TRUE if `x` has loops, otherwise FALSE. (No
#' computation is done here; the function simply returns the value of
#' `x$UNBROKEN_LOOPS`).
#'
#' * `hasInbredFounders(x)` returns TRUE is founder inbreeding is specified for
#' `x` and at least one founder has positive inbreeding coefficient. See
#' [founderInbreeding()] for details.
#'
#' * `hasSelfing(x)` returns TRUE if the pedigree contains selfing events. This
#' is recognised by father and mother begin equal for some child. (Note that for
#' this to be allowed, the gender code of the parent must be 0.)
#'
#' * `hasCommonAncestor(x)` computes a logical matrix `A` whose entry `A[i,j]`
#' is TRUE if pedigree members i and j have a common ancestor in `x`, and FALSE
#' otherwise. By convention, `A[i,i]` is TRUE for all i.
#'
#' * `subnucs(x)` returns a list of all nuclear sub-pedigrees of `x`, wrapped as
#' `nucleus` objects. Each nucleus is a list with entries `father`, `mother` and
#' `children`.
#'
#' * `peelingOrder(x)` calls `subnucs(x)` and extends each entry with a `link`
#' individual, indicating a member linking the nucleus to the remaining
#' pedigree. One application of this function is the fact that if _fails_ to
#' find a complete peeling order if and only if the pedigree has loops. (In fact
#' it is called each time a new `ped` object is created by [ped()] in order to
#' detect loops.) The main purpose of the function, however, is to prepare for
#' probability calculations in other packages, as e.g. in
#' `pedprobr::likelihood`.
#'
#' @examples
#' x = fullSibMating(1)
#' stopifnot(pedsize(x) == 6)
#' stopifnot(hasUnbrokenLoops(x))
#' stopifnot(generations(x) == 3)
#'
#' # All members have common ancestors except the grandparents
#' CA = hasCommonAncestor(x)
#' stopifnot(!CA[1,2], !CA[2,1], sum(CA) == length(CA) - 2)
#'
#' # Effect of breaking the loop
#' y = breakLoops(x)
#' stopifnot(!hasUnbrokenLoops(y))
#' stopifnot(pedsize(y) == 7)
#'
#' # A pedigree with selfing (note the necessary `sex = 0`)
#' z1 = singleton(1, sex = 0)
#' z2 = addChildren(z1, father = 1, mother = 1, nch = 1)
#' stopifnot(!hasSelfing(z1), hasSelfing(z2))
#'
#' # Nucleus sub-pedigrees
#' stopifnot(length(subnucs(z1)) == 0)
#' peelingOrder(cousinPed(1))
#'
#' # Plot with generation numbers as labels
#' w = cousinPed(1)
#' g = generations(w, what = "indiv")
#' labs = setNames(labels(w), g)
#' plot(w, labs = labs)
#'
#' # ... compare with
#' plot(relabel(w, "generations"))
#'
#' @name ped_utils
NULL
#' @rdname ped_utils
#' @export
pedsize = function(x) {
if(is.ped(x))
return(length(x$ID))
else if (is.pedList(x))
return(vapply(x, function(comp) length(comp$ID), 1L))
stop2("Input to `pedsize()` must be a `ped` object or a list of such")
}
#' @rdname ped_utils
#' @export
generations = function(x, what = c("max", "compMax", "indiv", "depth")) {
if(is.pedList(x)) {
what = match.arg(what)
usewhat = if(what == "compMax") "max" else what
res = lapply(x, function(comp) generations(comp, what = usewhat))
return(switch(what,
max = max(unlist(res)),
compMax = unlist(res),
indiv = , depth = res))
}
what = match.arg(what)
if(what == "indiv") {
p = .pedAlignment(x)$plist
gen = rep(seq_along(p$n), p$n)
oldIdx = unlist(lapply(seq_along(p$n), function(i) p$nid[i, 1:p$n[i]]))
dups = duplicated(oldIdx)
if(any(dups)) {
oldIdx = oldIdx[!dups]
gen = gen[!dups]
}
# Check for error
if(length(oldIdx) != length(x$ID)) {
warning("Alignment error; cannot find generation numbers")
return(NULL)
}
names(gen) = x$ID[oldIdx]
return(gen[x$ID])
}
# Depth definition: dp[i] = 1 + max(dp[parents])
# NB: Algorithm requires "parentsBeforeChildren".
# Iteration starting with depth = 1 for founders
xorig = x
x = parentsBeforeChildren(x)
N = length(x$ID)
FIDX = x$FIDX
MIDX = x$MIDX
NONFOU = which(FIDX > 0)
dp = rep(1L, N) # depth = 1 for founders
for(i in NONFOU)
dp[i] = 1L + max(dp[c(FIDX[i], MIDX[i])])
if(what == "max")
return(max(dp))
if(what == "depth") {
names(dp) = x$ID
dp[xorig$ID]
}
}
#' @rdname ped_utils
#' @export
hasUnbrokenLoops = function(x) {
if(is.pedList(x))
return(any(vapply(x, hasUnbrokenLoops, logical(1))))
unbrokenLoops = x$UNBROKEN_LOOPS
if(length(unbrokenLoops) == 1 && is.na(unbrokenLoops)) {
# Detect loop by trying to find a peeling order
nucs = peelingOrder(x)
lastnuc_link = nucs[[length(nucs)]]$link
unbrokenLoops = is.null(lastnuc_link)
}
isTRUE(unbrokenLoops)
}
#' @rdname ped_utils
#' @export
hasInbredFounders = function(x, chromType = "autosomal") {
if(is.pedList(x))
return(any(vapply(x, hasInbredFounders, logical(1), chromType = chromType)))
if(is.null(x$FOUNDER_INBREEDING))
return(FALSE)
finb = founderInbreeding(x, named = TRUE, chromType = chromType)
# If X: only females interesting (males are always 1)
if(chromType == "x")
finb = finb[getSex(x, names(finb)) == 2]
any(finb > 0)
}
#' @rdname ped_utils
#' @export
hasSelfing = function(x) {
if(is.pedList(x))
return(any(vapply(x, hasSelfing, logical(1))))
any(x$FIDX != 0 & x$FIDX == x$MIDX)
}
#' @rdname ped_utils
#' @export
hasCommonAncestor = function(x) {
if(!is.ped(x))
stop2("Input is not a `ped` object: ", x)
n = pedsize(x)
labs = labels(x)
A = matrix(FALSE, ncol = n, nrow = n, dimnames = list(labs, labs))
FOU = founders(x, internal = TRUE)
for(i in FOU) {
# vector of all descendants of i, including i
desc = c(i, descendants(x, i, internal = TRUE))
A[fast.grid(rep(list(desc), 2))] = TRUE
}
A
}
validate_sex = function(sex, nInd, zero_allowed = TRUE) {
if(length(sex) == 0)
stop2(sprintf("`%s` cannot be empty", deparse(substitute(sex))))
if(length(sex) > nInd)
stop2(sprintf("`%s` is longer than the number of individuals",
deparse(substitute(sex))))
codes = if(zero_allowed) 0:2 else 1:2
sex_int = suppressWarnings(as.integer(sex))
ok = (sex == sex_int) & sex_int %in% codes
if(!all(ok))
stop2("Illegal sex: ", unique(sex[!ok]),
". [0 = NA; 1 = male; 2 = female]")
invisible(rep_len(sex_int, nInd))
}
#' @rdname ped_utils
#' @export
subnucs = function(x) {
if (is.singleton(x))
return(list())
n = pedsize(x)
labs = labels(x)
seqn = seq_len(n)
FIDX = x[["FIDX"]]
MIDX = x[["MIDX"]]
# Indices of unique parent couples
p_pairs_idx = seqn[FIDX + MIDX > 0 & !duplicated.default(FIDX*(n+1) + MIDX)]
# List all nucs
lapply(rev.default(p_pairs_idx), function(j) {
nuc = list(father = FIDX[j], mother = MIDX[j], children = seqn[FIDX == FIDX[j] & MIDX == MIDX[j]])
class(nuc) = "nucleus"
attr(nuc, "labels") = labs
nuc
})
}
#' @rdname ped_utils
#' @export
peelingOrder = function(x) {
# output: list of nuclear subfamilies. Format for each nuc:
# list(father,mother,children,link), where link = 0 for the last nuc.
nucs = subnucs(x)
if(length(nucs) == 0)
return(nucs)
peeling = vector("list", Nnucs <- length(nucs))
i = k = 1
while (length(nucs)) {
# Start searching for a nuc with only one link to others
nuc = nucs[[i]]
# Identify links to other remaining nucs
if(length(nucs) > 1) {
nucmembers = c(nuc$father, nuc$mother, nuc$children)
links = nucmembers[nucmembers %in% unlist(nucs[-i], use.names = FALSE)]
}
else
links = 0 # if nuc is the last
# If only one link: move nuc to peeling list, and proceed
if (length(links) == 1) {
nuc$link = links
peeling[[k]] = nuc
nucs[i] = NULL
i = 1
k = k+1
}
else {
if (i == length(nucs)) { # LOOP! Include remaining nucs without 'link', and break.
peeling[k:Nnucs] = nucs
break
}
# Otherwise try next remaining nuc
i = i+1
}
}
peeling
}
#' S3 methods
#'
#' @param x An object
#' @param ... Not used
#' @export
print.nucleus = function(x, ...) {
labs = attr(x, 'labels')
fa = x$father
mo = x$mother
ch = x$children
link = x$link
if (is.null(link))
linktext = "NO LINK ASSIGNED"
else if (link == 0)
linktext = "0 (end of chain)"
else {
who = if (link == fa) "father"
else if (link == mo) "mother"
else if (link %in% ch) "child"
else stop2("Erroneous link individual (not a member of nucleus): ", link)
linktext = sprintf("%s (%s)", labs[link], who)
}
cat(sprintf("Nucleus: Fa/mo/ch = %s/%s/%s\nLink individual: %s\n",
labs[fa], labs[mo], paste(labs[ch], collapse = ","), linktext))
}
hasNumLabs = function(x) {
# Returns TRUE if the labels of x are coercible to integers
labs = labels(x)
numlabs = suppressWarnings(as.character(as.integer(labs)))
isTRUE(all(labs == numlabs))
}
# New implementation of `hasNumLabs`; takes either a character or ped(list)
.isIntegral = function(labs) {
if(!is.character(labs))
labs = labels(labs)
intlabs = suppressWarnings(as.character(as.integer(labs)))
isTRUE(all(labs == intlabs))
}
# Utility for creating new labels
generateLabs = function(x, n = 1, avoid = NULL, prefix = "") {
labs = if(is.character(x)) x else labels(x)
# Simple strategy (no need to be clever about this...)
# Generate enough candidates and drop the taken ones
taken = c(labs, avoid)
cand = as.character(seq_len(length(taken) + n))
if(nzchar(prefix))
cand = paste0(prefix, cand)
.mysetdiff(cand, taken)[seq_len(n)]
}
#' Pedigree component
#'
#' Given a list of `ped` objects (called pedigree components), and a vector of
#' ID labels, find the index of the component holding each individual.
#'
#' @param x A `ped` object, or a list of such.
#' @param ids A vector of ID labels (coercible to character).
#' @param checkUnique A logical, by default FALSE. If TRUE, an error is raised
#' if any element of `ids` occurs more than once in `x`.
#' @param errorIfUnknown A logical, by default FALSE. If TRUE, the function
#' stops with an error if not all elements of `ids` are found in `x`.
#'
#' @return An integer vector of the same length as `ids`, with NA entries where
#' the corresponding label was not found in any of the components.
#'
#' @seealso [internalID()]
#'
#' @examples
#' x = list(nuclearPed(1), singleton(id = "A"))
#' getComponent(x, c(3, "A"))
#'
#' @export
getComponent = function(x, ids, checkUnique = FALSE, errorIfUnknown = FALSE) {
# If simple ped: Handle separately and return early
if(is.ped(x)) {
comp = rep(1L, length(ids))
if(any(unkn <- !ids %in% x$ID)) {
if(errorIfUnknown)
stop2("Unknown ID label: ", ids[unkn])
comp[unkn] = NA_integer_
}
return(comp)
}
if(!is.pedList(x)) {
print(x)
stop2("Input is not a `ped` object (or a list of such)")
}
# List labels of each component
labList = labels(x, unlist = FALSE)
# A single vector with all labels
labVec = unlist(labList, recursive = FALSE, use.names = FALSE)
# Check for duplicates if indicated
if(checkUnique) {
v = labVec[labVec %in% ids]
if(dup <- anyDuplicated.default(v))
stop2("ID label is not unique: ", v[dup])
}
# Vector of same length as labVec, with component index for each member
compi = rep(seq_along(labList), times = lengths(labList))
# Match input ids against label vector
idx = match(ids, labVec)
if(anyNA(idx) && errorIfUnknown)
stop2("Unknown ID label: ", ids[is.na(idx)])
# Return comp idx of the input ids, including NA if not present
compi[idx]
}
topolOrder = function(x) {
N = pedsize(x)
env = new.env()
logvec = logical(N)
names(logvec) = x$ID
isNonFou = logvec
isNonFou[nonfounders(x)] = TRUE
env$visited = logvec
env$ord = character(N)
env$k = 1
dfs = function(x, id, env) {
env$visited[id] = TRUE
if(isNonFou[id]) {
fa = father(x, id)
if(!env$visited[fa])
dfs(x, fa, env)
mo = mother(x, id)
if(!env$visited[mo])
dfs(x, mo, env)
}
k = env$k
env$ord[k] = id
#message("Added ", id, " in position ", N-k)
env$k = k + 1
}
lvs = sort.int(leaves(x), method = "quick")
for(id in lvs)
dfs(x, id, env)
# Return ordered labels
env$ord
}
reorderSimple = function(x, neworder, ids = NULL) {
if(!is.ped(x))
stop2("`reorderQuick` only works for `ped` objects")
if(length(neworder) != length(x$ID))
stop2("Wrong length of `neworder`")
if(is.character(neworder))
neworder = match(neworder, x$ID)
ID = x$ID[neworder]
FIDX = x$FIDX[neworder]
MIDX = x$MIDX[neworder]
SEX = x$SEX[neworder]
nonf = FIDX > 0
FIDX[nonf] = match(FIDX[nonf], neworder)
MIDX[nonf] = match(MIDX[nonf], neworder)
# Mask all but `ids` indivs
if(!is.null(ids)) {
mask = !ID %in% ids
ID[mask] = "*"
SEX[mask] = "*"
}
list(ID, FIDX, MIDX, SEX)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.