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)
drones <- createDrones(x = basePop[1], nInd = 1000)
droneGroups <- pullDroneGroupsFromDCA(drones, n = 1, nDrones = nFathersPoisson)
colony <- createColony(basePop[2])
colony <- cross(x = colony, drones = droneGroups[[1]])
colony <- buildUp(x = colony, nWorkers = 6, nDrones = 3)
geno <- getSegSiteGeno(colony, collapse = TRUE)
sex <- getCasteSex(x = colony, collapse = TRUE)
GRM <- calcBeeGRMIbs(x = geno, sex = sex)
# You can visualise this matrix with the function image() from the package 'Matrix'
#Look at the diagonal at the relationship matrix
x <- diag(GRM)
hist(x)
summary(x)
#Look at the off-diagonal at the relationship matrix
x <- GRM[lower.tri(x = GRM, diag = FALSE)]
hist(x)
summary(x)
# Compare relationship between castes
ids <- getCasteId(colony)
idQueen <- ids$queen
idWorkers <- ids$workers
idDrones <- ids$drones
# Queen vs others
GRM[idQueen, idWorkers]
GRM[idQueen, idDrones]
# Workers vs worker
GRM[idWorkers, idWorkers]
# Workers vs drones
GRM[idWorkers, idDrones]
# Calculating allele frequencies ourselves (say, to "shift" base population)
aF <- calcBeeAlleleFreq(x = geno, sex = sex)
hist(aF)
GRM2 <- calcBeeGRMIbs(x = geno, sex = sex, alleleFreq = aF)
stopifnot(identical(GRM2, GRM))
# You can also create relationships with pooled genomes
pooledGenoW <- getPooledGeno(getWorkersSegSiteGeno(colony),
type = "mean",
sex = getCasteSex(colony, caste="workers"))
queenGeno <- getQueenSegSiteGeno(colony)
# Compute relationship between pooled workers genotype and the queen
calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"))
# You can now compare how this compare to relationships between the queen
# individual workers!
Run the code above in your browser using DataLab