Extending getSegSiteGeno function to pull genotype instead of a list in multicolony setting

[DrMulusewFikere]

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

[gregorgorjanc]

@DrMulusewFikere I am looking at this. So, the getSegSiteGeno() is designed such that we can work on a Pop object (say a group of individuals), a Colony (where we have multiple Pop slots), or a MultiColony (where we have multiple Colony "slots"). Here is example:

founderGenomes <- quickHaplo(nInd = 20, nChr = 3, segSites = 10)
SP <- SimParamBee$new(founderGenomes,
                      nWorkers = 50,
                      nDrones = 10,
                      csdChr = 1,
                      nCsdAlleles = 3)
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)
str(getSegSiteGeno(apiary[[1]]@queen))

giving

 int [1, 1:30] 2 0 1 0 1 1 1 1 0 1 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr "1"
  ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...

then

> str(getSegSiteGeno(apiary[[1]]))
List of 5
 $ queen       : int [1, 1:30] 2 0 1 0 1 1 1 1 0 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr "1"
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ fathers     : num [1:17, 1:30] 1 1 1 0 1 0 0 1 1 0 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:17] "30" "139" "24" "165" ...
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ workers     : int [1:38, 1:30] 2 2 2 2 2 2 1 2 2 2 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:38] "328" "336" "340" "324" ...
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ drones      : num [1:10, 1:30] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:10] "380" "378" "374" "373" ...
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ virginQueens: NULL

and finally

> str(getSegSiteGeno(apiary))
List of 20
 $ 1 :List of 5
  ..$ queen       : int [1, 1:30] 2 0 1 0 1 1 1 1 0 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr "1"
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ fathers     : num [1:17, 1:30] 1 1 1 1 1 0 1 0 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:17] "70" "244" "36" "56" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ workers     : int [1:38, 1:30] 2 2 1 2 2 2 2 2 2 2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:38] "362" "370" "367" "324" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ drones      : num [1:10, 1:30] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:10] "376" "373" "380" "379" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ virginQueens: NULL
 $ 2 :List of 5
  ..$ queen       : int [1, 1:30] 1 0 1 0 0 0 1 0 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr "2"
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ fathers     : num [1:13, 1:30] 0 1 1 0 1 1 1 0 1 0 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:13] "241" "62" "111" "137" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ workers     : int [1:31, 1:30] 1 2 1 0 1 0 0 1 1 2 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:31] "418" "387" "404" "402" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ drones      : num [1:10, 1:30] 1 0 0 1 0 0 1 1 1 1 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:10] "433" "436" "435" "440" ...
  .. .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
  ..$ virginQueens: NULL
...

But, we did anticipate that people will want one block of genotypes, so we added the collapse argument to getSegSiteGeno(). So, one could do:

> str(getSegSiteGeno(apiary, caste = "queen", collapse = TRUE))
 int [1:20, 1:30] 2 1 0 1 2 2 1 1 1 1 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "1" "2" "3" "4" ...
  ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...

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:

> tmp <- getSegSiteGeno(x = apiary, caste = "workers")
> str(tmp)
List of 20
 $ 1 : int [1:38, 1:30] 2 2 1 2 1 1 2 2 2 2 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:38] "336" "362" "334" "364" ...
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ 2 : int [1:31, 1:30] 2 2 0 1 1 1 2 2 1 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:31] "389" "413" "407" "418" ...
...

We get a list of genotype matrices. So, we should be able to just run getPooledGeno() on each of these:

> lapply(X = tmp, FUN = getPooledGeno)
Error in if (type == "mean") { : argument is of length zero
In addition: Warning message:
In FUN(X[[i]], ...) :
  Argument sex is NULL. Assuming that all individuals are female/diploid!

Right, that's not what I expected! But, looking at the getPooledGeno() I see this will work (will improve the error message)

> str(lapply(X = tmp, FUN = getPooledGeno, type = "mean"))
List of 20
 $ 1 : num [1, 1:30] 1.7632 0.0789 0.8947 0.2105 1.2368 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...
 $ 2 : num [1, 1:30] 1.161 0.452 1.258 0.226 0.387 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...

So, we now have a list of one row matrices (with mean allele dosage of workers per colony) which we can rbind() using:

> tmpPooledGeno = lapply(X = tmp, FUN = getPooledGeno, type = "mean")
There were 20 warnings (use warnings() to see them)
> str(do.call(rbind, tmpPooledGeno))
 num [1:20, 1:30] 1.763 1.161 0.618 0.941 1.833 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:30] "1_1" "1_2" "1_3" "1_4" ...

giving us a matrix with colonies in rows and markers in columns. Alternative would be:

> str(t(sapply(X = tmp, FUN = getPooledGeno, type = "mean")))
 num [1:20, 1:30] 1.763 1.161 0.618 0.941 1.833 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "1" "2" "3" "4" ...
  ..$ : NULL
There were 20 warnings (use warnings() to see them)

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, and MultiColony, 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.

[gregorgorjanc]

@DrMulusewFikere @janaobsteter I have changed the default to type = "mean and handle that argument better, and have also included the above into the examples section of getPooledGeno() as this will be much sought for thing. With this in mind I will change this into a discussion point.