# check that kinship matrices are square
# and have same row and column IDs
#
# returns vector of IDs
check_kinship <-
function(kinship, n_chr)
{
if(is_kinship_decomposed(kinship)) {
if(is_kinship_list(kinship)) {
if(length(kinship) != n_chr)
stop("length(kinship) != no. chromosomes (", n_chr, ")")
return(rownames(kinship[[1]]$vectors))
}
else {
return(rownames(kinship$vectors))
}
}
if(!is_kinship_list(kinship)) { # one kinship matrix
stopifnot(nrow(kinship) == ncol(kinship))
stopifnot( all(rownames(kinship) == colnames(kinship)) )
return(rownames(kinship))
} else {
if(length(kinship) != n_chr)
stop("length(kinship) != no. chromosomes (", n_chr, ")")
kinship_square <- vapply(kinship, function(mat) nrow(mat) == ncol(mat), TRUE)
stopifnot( all(kinship_square) )
kinship_sameIDs <- vapply(kinship, function(mat) (nrow(mat) == nrow(kinship[[1]])) &&
all((rownames(mat) == rownames(kinship[[1]])) &
(colnames(mat) == colnames(kinship[[1]])) &
(rownames(mat) == colnames(kinship[[1]]))), TRUE)
if(!all(kinship_sameIDs))
stop("All kinship matrices should be the same size ",
"and have the same row and column names")
return(rownames(kinship[[1]]))
}
}
# multiply kinship matrix by 2
# see Almasy & Blangero (1998) https://doi.org/10.1086/301844
#
# This can also handle the case of "loco", and of having eigen decomposition pre-computed
double_kinship <-
function(kinship)
{
if(is.null(kinship)) return(NULL)
if(is_kinship_decomposed(kinship)) { # already did decomposition
if(is_kinship_list(kinship)) { # list of decomposed kinship matrices
kinship <- lapply(kinship, function(a) { a$values <- 2*a$values; a })
}
else { # one decomposed kinship matrix
kinship$values <- 2*kinship$values
}
}
else {
if(is.list(kinship))
kinship <- lapply(kinship, function(a) a*2)
else kinship <- 2*kinship
}
kinship
}
# check if alread decomposed
is_kinship_decomposed <-
function(kinship)
{
decomp <- attr(kinship, "eigen_decomp")
(!is.null(decomp) && decomp) || # should have attribute
(length(kinship)==2 && all(names(kinship) == c("values", "vectors"))) || # single-chr case missing attribute
(is.list(kinship) && all(vapply(kinship, length, 1)==2) && all(vapply(kinship, function(a) all(names(a)==c("values", "vectors")), TRUE))) # multi-chr case
}
# is kinship a list with (potentially) multiple chromosomes
is_kinship_list <-
function(kinship)
{
if(is_kinship_decomposed(kinship)) {
if(length(kinship)==2 && all(names(kinship) == c("values", "vectors"))) { # just one chromosome
return(FALSE)
}
else return(TRUE)
}
else {
return(is.list(kinship))
}
}
# check that kinship concerns one chromosome
# and remove outer list if necessary
check_kinship_onechr <-
function(kinship)
{
if(is_kinship_list(kinship)) {
if(length(kinship) > 1)
stop("kinship should concern just one chromosome")
decomp <- attr(kinship, "eigen_decomp") ## preserve attribute
kinship <- kinship[[1]]
attr(kinship, "eigen_decomp") <- decomp
}
kinship
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.