Check what caste we assign to diploid drones!?

[gregorgorjanc]

We create and kill diploid drones, which means they appear in the SP$pedigree. This is fine and we can't go around it. But, what SP$caste value do we store for them? Using SP$pedigree[SP$caste != "diploidDrone", ] would a neat way of getting rid of them from SP$pedigree for downstream analyses.

[gregorgorjanc]

I am looking into this. The key function is createCastePop(). The key part is

    if (caste == "workers") {
      ret <- vector(mode = "list", length = 2)
      names(ret) <- c("workers", "nHomBrood")
      workers <- combineBeeGametes(
        queen = getQueen(x), drones = getFathers(x),
        nProgeny = nInd, simParamBee = simParamBee
      )
      if (isCsdActive(simParamBee = simParamBee)) {
        sel <- isCsdHeterozygous(pop = workers, simParamBee = simParamBee)
        ret$workers <- workers[sel]
        ret$nHomBrood <- nInd - sum(sel)
        if (exact) {
          if (nInd(ret$workers) < nInd) {
            nMiss <- nInd - nInd(ret$workers)
            while (0 < nMiss) {
              workers <- combineBeeGametes(
                queen = getQueen(x),
                drones = getFathers(x),
                nProgeny = nMiss,
                simParamBee = simParamBee
              )
              sel <- isCsdHeterozygous(pop = workers, simParamBee = simParamBee)
              ret$workers <- c(ret$workers, workers[sel])
              ret$nHomBrood <- ret$nHomBrood + sum(!sel)
              nMiss <- nInd - nInd(ret$workers)
            }
          }
        }
      } else {
        ret$workers <- workers
        ret$nHomBrood <- NA
      }
      ret$workers@sex[] <- "F"
      simParamBee$addToCaste(id = ret$workers@id, caste = "W")

So, at the moment, we could be having SP$caste shorter than SP$pedigree because we only add non-csd-homozygous individuals to it (ret$workers@id), while the homozygous drones are skipped.

Now, we have two options:

1. We keep it as it is, we do save individual id and caste in SP$caste - to see this I have run an example from the createCastePop() and I get this:

> length(SP$caste)
[1] 4247
> SP$caste[4247]
   4324 
workers 
Levels: queen fathers workers drones virginQueens

As we can see above, we have 4247 entries in the SP$caste, but the last individual saved in there is individual with id 4324, which means there were 4324-4247=77 homozygous drones.

So, to get a pedigree with homozygous drones removed, we could do

> str(SP$pedigree[as.integer(names(SP$caste)), ]) # pedigree without homo drones
 int [1:4247, 1:3] 0 0 0 0 0 0 0 0 1 1 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:4247] "1" "2" "3" "4" ...
  ..$ : chr [1:3] "mother" "father" "isDH"

> str(SP$pedigree) # pedigree with homo drones
 int [1:4324, 1:3] 0 0 0 0 0 0 0 0 1 1 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:4324] "1" "2" "3" "4" ...
  ..$ : chr [1:3] "mother" "father" "isDH"

2. Add homo drones in the SP$caste for completeness

So, technically all is fine - no alarm. But will someone get confused? Thoughts?

[janaobsteter]

Hm, this is actually good news - I didn't know we had the functionality to extract the non-homozygous individuals from the pedigree. For us, it makes no difference since we're familiar with the simulator and data structures - but I don't know what would confuse an average user more. @LStrachan , @JerBub , thoughts?

[gregorgorjanc]

We should probably remove them at least in the relationship paper to compare apples with apples? With regards, Gregor

…

On Mon, Nov 7 2022 at 7:24 am, janaobsteter ***@***.***> wrote: Hm, this is actually good news - I didn't know we had the functionality to extract the non-homozygous individuals from the pedigree. For us, it makes no difference since we're familiar with the simulator and data structures - but I don't know what would confuse an average user more. @LStrachan <https://github.com/LStrachan> , @JerBub <https://github.com/JerBub> , thoughts? — Reply to this email directly, view it on GitHub <#428 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAFRGSSYNB4V75TLLJ4575LWHCVDNANCNFSM6AAAAAARWHN5JA> . You are receiving this because you were assigned.Message ID: ***@***.***>

[janaobsteter]

I was thinking about that the other day - and I'm not so sure about it anymore. Because, in reality, that is the expected relationship - and you don't actually know based on the expected relationship which of them will die due to csd homozygosity.
On the other hand, we are currently not comparing the same set of individuals - but again, is that actually fine?

[janaobsteter]

Create a test with an inbred colony - create n workers and drones and check whether the pedigree / caste / recHist / sex length is the same.

[janaobsteter]

Test is added in #625. Also, I'm pasting the example code here. SP$pedigree, SP$caste , and SP$recHist all match in length - that is, they all include the homozygous (killed) individuals. The actual number of bees alive differs from that.

founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100)
SP <- SimParamBee$new(founderGenomes)
SP$nThreads = 1 #8
SP$setTrackPed(T)
SP$setTrackRec(T)
#registerDoParallel(cores = 8) # Do this if selecting > 1 cores

nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
basePop <- createVirginQueens(founderGenomes, simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)

drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); (basePop@nInd + drones@nInd)
fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP)

# Create a Colony class object
colony <- createColony(x = basePop[1], simParamBee = SP)
colony <- cross(colony, drones = fatherGroups[[1]], simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)


colony <- buildUp(x = colony, nWorkers = 200, nDrones = 50, simParamBee = SP)
# This is the length of the pedigree and caste - they match
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
# They also match with the expected/specified number of created individuals
(basePop@nInd + drones@nInd + 200 + 50)
# But this is the actual number of created individuals
(basePop@nInd + drones@nInd + nWorkers(colony) + nDrones(colony))

[gregorgorjanc]

@janaobsteter ok, we clearly give them a caste. Out of curiosity, which caste is it (for the diploid drones)?

[gregorgorjanc]

I was curious about this and run the code against the PR #625 and got this

founderGenomes <- quickHaplo(nInd = 8, nChr = 1, segSites = 100)
SP <- SimParamBee$new(founderGenomes)
SP$nThreads = 1 #8
SP$setTrackPed(T)
SP$setTrackRec(T)
#registerDoParallel(cores = 8) # Do this if selecting > 1 cores

nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
# 0 x3
basePop <- createVirginQueens(founderGenomes, simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
# 8 x3
SP$caste
# all virging queens

drones <- createDrones(x = basePop[1], nInd = 1000, simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); (basePop@nInd + drones@nInd)
# 1008 x3
fatherGroups <- pullDroneGroupsFromDCA(drones, n = 10, nDrones = 10, simParamBee = SP)
table(SP$caste)
# drones virginQueens
#   1000            8

# Create a Colony class object
colony <- createColony(x = basePop[1], simParamBee = SP)
colony <- cross(colony, drones = fatherGroups[[1]], simParamBee = SP)
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
# 1008 x3
table(SP$caste)
# drones      fathers        queen virginQueens
#    990           10            1            7

colony <- buildUp(x = colony, nWorkers = 200, nDrones = 250, simParamBee = SP)
# This is the length of the pedigree and caste - they match
nrow(SP$pedigree); length(SP$caste); length(SP$recHist)
# 1458 x3
# They also match with the expected/specified number of created individuals
(basePop@nInd + drones@nInd + 200 + 250)
# 1458
# But this is the actual number of created individuals
(basePop@nInd + drones@nInd + nWorkers(colony) + nDrones(colony))
# 1358
table(SP$caste)
# drones      fathers        queen virginQueens workers 
#   1240           10            1            7     200

colony@queen@misc$nHomBrood # 100
colony@queen@misc$pHomBrood # 0.5

nDrones(colony) # 250
nWorkers(colony) # 100

So, the answer to the original question: "Check what caste we assign to diploid drones!?" is that we give them the caste of worker, which sounds odd, but is all correct biologically/codewise - we set out co create 200 workers above, but we only got out 100 because 100 were homozygous at the csd locus (which matches with the pHomBrood = 0.5). Thence, the code all works very well.

[gregorgorjanc]

@janaobsteter the above example could be a unit test case of generating >0 nHomBrood?