Calculate Genomic Relatedness Matrix (GRM) for honeybees from Identical By State genomic data

Source: R/Functions_L0_auxilary.R

Level 0 function that returns Genomic Relatedness Matrix (GRM) for honeybees from Identical By State genomic data (bi-allelic SNP represented as allele dosages) following the method for the sex X chromosome (Druet and Legarra, 2020)

calcBeeGRMIbs(x, sex, alleleFreq = NULL)

calcBeeAlleleFreq(x, sex)

Arguments

x

matrix of genotypes represented as allele dosage coded as 0, 1, or 2 in females (queens or workers) and as 0 or 1 in males (fathers or drones); individuals are in rows and sites are in columns; no missing values are allowed (this is not checked - you will get NAs!)

sex

character vector denoting sex for individuals with genotypes in x - "F" for female and "M" for male

alleleFreq

numeric, vector of allele frequencies for the sites in x; if NULL, then calcBeeAlleleFreq is used

Value

matrix of genomic relatedness coefficients

Functions

  • calcBeeAlleleFreq(): Calculate allele frequencies from honeybee genotypes

References

Druet and Legarra (2020) Theoretical and empirical comparisons of expected and realized relationships for the X-chromosome. Genetics Selection Evolution, 52:50 doi:/10.1186/s12711-020-00570-6

Examples

founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100)
SP <- SimParamBee$new(founderGenomes)
SP$nThreads = 1L
SP$setTrackRec(TRUE)
SP$setTrackPed(isTrackPed = TRUE)

basePop <- createVirginQueens(founderGenomes)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
drones <- createDrones(x = basePop[1], nInd = 1000)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
droneGroups <- pullDroneGroupsFromDCA(drones, n = 1, nDrones = nFathersPoisson)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
colony <- createColony(basePop[2])
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
colony <- cross(x = colony, drones = droneGroups[[1]])
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found

geno <- getSegSiteGeno(colony, collapse = TRUE)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
sex <- getCasteSex(x = colony, collapse = TRUE)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found

GRM <- calcBeeGRMIbs(x = geno, sex = sex)
#> Error in eval(expr, envir, enclos): object 'geno' not found
# You can visualise this matrix with the function image() from the package 'Matrix'

#Look at the diagonal at the relationship matrix
x <- diag(GRM)
#> Error in eval(expr, envir, enclos): object 'GRM' not found
hist(x)
#> Error in eval(expr, envir, enclos): object 'x' not found
summary(x)
#> Error in eval(expr, envir, enclos): object 'x' not found

#Look at the off-diagonal at the relationship matrix
x <- GRM[lower.tri(x = GRM, diag = FALSE)]
#> Error in eval(expr, envir, enclos): object 'GRM' not found
hist(x)
#> Error in eval(expr, envir, enclos): object 'x' not found
summary(x)
#> Error in eval(expr, envir, enclos): object 'x' not found

# Compare relationship between castes
ids <- getCasteId(colony)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
idQueen <- ids$queen
#> Error in eval(expr, envir, enclos): object 'ids' not found
idWorkers <- ids$workers
#> Error in eval(expr, envir, enclos): object 'ids' not found
idDrones <- ids$drones
#> Error in eval(expr, envir, enclos): object 'ids' not found

# Queen vs others
GRM[idQueen, idWorkers]
#> Error in eval(expr, envir, enclos): object 'GRM' not found
GRM[idQueen, idDrones]
#> Error in eval(expr, envir, enclos): object 'GRM' not found

# Workers vs worker
GRM[idWorkers, idWorkers]
#> Error in eval(expr, envir, enclos): object 'GRM' not found

# Workers vs drones
GRM[idWorkers, idDrones]
#> Error in eval(expr, envir, enclos): object 'GRM' not found

# Calculating allele frequencies ourselves (say, to "shift" base population)
aF <- calcBeeAlleleFreq(x = geno, sex = sex)
#> Error in eval(expr, envir, enclos): object 'geno' not found
hist(aF)
#> Error in eval(expr, envir, enclos): object 'aF' not found
GRM2 <- calcBeeGRMIbs(x = geno, sex = sex, alleleFreq = aF)
#> Error in eval(expr, envir, enclos): object 'geno' not found
stopifnot(identical(GRM2, GRM))
#> Error in eval(expr, envir, enclos): object 'GRM2' not found

# You can also create relationships with pooled genomes
pooledGenoW <- getPooledGeno(getWorkersSegSiteGeno(colony),
                             type = "mean",
                             sex = getCasteSex(colony, caste="workers"))
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
queenGeno <- getQueenSegSiteGeno(colony)
#> Error in get(x = "SP", envir = .GlobalEnv): object 'SP' not found
# Compute relationship between pooled workers genotype and the queen
calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"))
#> Error in eval(expr, envir, enclos): object 'queenGeno' not found
# You can now compare how this compare to relationships between the queen
# individual workers!