Repository navigation
Extending getSegSiteGeno function to pull genotype instead of a list in multicolony setting #590
Description
In multicolony setting "getSegSiteGeno" function generates a list instead of raw genotype data
Example
Create a founder genome
founderGenomes <- quickHaplo(nInd = 20, nChr = 3, segSites = 100)
Initialise SP
SP <- SimParamBee$new(founderGenomes,
nWorkers = 50,
nDrones = 10,
csdChr = 3,
nCsdAlleles = 32)
baseQueens <- createVirginQueens(founderGenomes)
apiary <- createMultiColony(x = baseQueens[1:20])
DCA <- createDrones(baseQueens[5:10],nInd = 50)
fun <- nFathersPoisson
droneGroups <- pullDroneGroupsFromDCA(DCA, n = 20,nDrones = fun)
apiary <- cross(apiary, drones = droneGroups)
apiary <- buildUp(apiary)
Below I have temporary addressed the issue and it works
start.time <- Sys.time()
library(parallel)
nColony <- nColonies(apiary)
process_colony <- function(colony_data) {
segSiteGeno <- getSegSiteGeno(colony_data, caste = "workers")
pooledGeno <- getPooledGeno(segSiteGeno,
type = "mean",
sex = getCasteSex(colony_data, caste = "workers"))
return(pooledGeno)
}
genoW <- lapply(seq_len(nColony), function(colony) {
process_colony(apiary@colonies[[colony]])
})
pooledGenoW <- do.call(rbind, genoW)
Another option
genoW <- vector("list", length = nColony)
nColony <- colonies(age_0)
for (colony in 1:nColony) {
pooledGenoW <- getPooledGeno(getSnpGeno(age_0, snpChip = 1, caste = "workers",collapse = T),type = "mean",
sex = getCasteSex(age_0@colonies[[colony]], caste = "workers"))
genoW[[colony]] <- pooledGenoW
}
pooledGenoW <- do.call(rbind, genoW)
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
Activity
[-]Extending getWorkersSegSiteGeno function to pull genotype instead of a list in multicolony setting[/-][+]Extending getSegSiteGeno function to pull genotype instead of a list in multicolony setting[/+]gregorgorjanc commentedon Dec 14, 2024
@DrMulusewFikere I am looking at this. So, the
getSegSiteGeno()is designed such that we can work on aPopobject (say a group of individuals), aColony(where we have multiplePopslots), or aMultiColony(where we have multipleColony"slots"). Here is example:giving
then
and finally
But, we did anticipate that people will want one block of genotypes, so we added the
collapseargument togetSegSiteGeno(). So, one could do:But I appreciate that you want to get pooled workers genotype per colony and the collapsed into a matrix across colonies!
So, how would we go about doing this in a clean way? Let's see what can already be done:
We get a list of genotype matrices. So, we should be able to just run
getPooledGeno()on each of these:Right, that's not what I expected! But, looking at the
getPooledGeno()I see this will work (will improve the error message)So, we now have a list of one row matrices (with mean allele dosage of workers per colony) which we can
rbind()using:giving us a matrix with colonies in rows and markers in columns. Alternative would be:
Does this help @DrMulusewFikere?
I appreciate this looks a bit involved, but it's because we work across colonies and we try to have a consistent way of working across
Pop,Colony, andMultiColony, so it gets a bit involved, but not too much.This will be a very common request, I think. So, we should probably add this to examples.
getPooledGeno(): Handle type arg, better default, and pooled example …
gregorgorjanc commentedon Dec 15, 2024
@DrMulusewFikere @janaobsteter I have changed the default to
type = "meanand handle that argument better, and have also included the above into the examples section ofgetPooledGeno()as this will be much sought for thing. With this in mind I will change this into a discussion point.