Skip to content

Instantly share code, notes, and snippets.

@kbroman
Created February 2, 2021 16:04
Show Gist options
  • Save kbroman/e09f8116401180b80005d47d4358764f to your computer and use it in GitHub Desktop.
Save kbroman/e09f8116401180b80005d47d4358764f to your computer and use it in GitHub Desktop.
Expand RIL data for individual-level analysis
# functions for expanding probs and kinship matrices
# for individual-level analysis of RILs
# create mapping, individuals -> lines
#
# n = sample size for each line
n2mapping <-
function(n)
{
mapping <- rep("", sum(n))
cur <- 0
for(i in seq_along(n)) {
mapping[cur + 1:n[i]] <- names(n)[i]
names(mapping)[cur + 1:n[i]] <- paste0(names(n)[i], "_", seq_len(n[i]))
cur <- cur + n[i]
}
mapping
}
# create mapping using individual names, assuming their like "strain_individual"
names2mapping <-
function(ind_names)
{
setNames( sapply(strsplit(ind_names, "_"), "[", 1), ind_names )
}
# expand genoprobs
#
# probs = output of calc_genoprob()
# mapping = vector of character strings corresponding to individual lines in probs
# names are the new names
expand_probs <-
function(probs, mapping)
{
for(i in seq_along(probs)) {
rn <- rownames(probs[[i]])
if(!all(!is.na(mapping) & mapping %in% rn)) {
missing <- mapping[is.na(mapping) | !(mapping %in% rn)]
stop("Some strains not in probs: ", paste(missing, collapse=" "))
}
probs[[i]] <- probs[[i]][match(mapping, rownames(probs[[i]])),,]
rownames(probs[[i]]) <- names(mapping)
}
probs
}
# expand kinship
expand_kinship <-
function(kinship, mapping)
{
if(is.list(kinship)) {
for(i in seq_along(kinship)) {
kinship[[i]] <- expand_kinship(kinship[[i]], mapping)
}
return(kinship)
}
rn <- rownames(kinship)
if(!all(!is.na(mapping) & mapping %in% rn)) {
missing <- mapping[is.na(mapping) | !(mapping %in% rn)]
stop("Some strains not in kinship: ", paste(missing, collapse=" "))
}
m <- match(mapping, rownames(kinship))
newk <- matrix(nrow=length(mapping), ncol=length(mapping))
dimnames(newk) <- list(names(mapping), names(mapping))
for(i in seq_along(m)) {
newk[i,] <- kinship[m[i], m]
}
newk
}
# create_plain_kinship
create_plain_kinship <-
function(mapping)
{
# create identity matrix for the strains
u <- unique(mapping)
k <- diag(rep(1, length(u)))
dimnames(k) <- list(u,u)
# expand
expand_kinship(k, mapping)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment