HighlanderLab / SIMplyBee

SIMplyBee R package extends AlphaSimR for simulating honeybee populations and breeding programmes
http://www.simplybee.info/
Other
0 stars 5 forks source link

Check what caste we assign to diploid drones!? #428

Open gregorgorjanc opened 1 year ago

gregorgorjanc commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 https://github.com/HighlanderLab/SIMplyBee/issues/428#issuecomment-1305188762, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAFRGSSYNB4V75TLLJ4575LWHCVDNANCNFSM6AAAAAARWHN5JA . You are receiving this because you were assigned.Message ID: @.***>

janaobsteter commented 1 year ago

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?