search

gaynorr / AlphaSimR

R package for breeding program simulations
https://gaynorr.github.io/AlphaSimR/
Other
43 stars 18 forks source link

Bug: reduceGenome(), doubleGenome(), and mergeGenome() don't pass info to global pedigree due to incomplete newPop() call #38

Closed gregorgorjanc closed 2 years ago

gregorgorjanc commented 2 years ago

We are using reduceGenome() and mergeGenome() to simulate a 'bee' cross between a diploid female and haploid males - we "reduce" genome of the female to simulate meiosis and then merge it with genomes of males. We noticed that the global pedigree stored in SimParam is not updated when we use these functions - we have all zeroes in the global pedigree. A hunt through the code indicates, that this is the way newPop() is called in the polyploids.R for reduceGenome(), doubleGenome(), and () compared to other crossing functions in crossing.R. For example, makeCross() uses

...
return(newPop(rawPop=rPop,
                mother=pop@id[crossPlan[,1]],
                father=pop@id[crossPlan[,2]],
                simParam=simParam,
                iMother=pop@iid[crossPlan[,1]],
                iFather=pop@iid[crossPlan[,2]],
                femaleParentPop=pop,
                maleParentPop=pop,
                hist=hist))
...

while mergeGenome() uses:

...
  return(newPop(rawPop=rPop,
                mother=mother,
                father=father,
                isDH=FALSE,
                simParam=simParam))
...

in addition mergeGenome() calls simParam$addToRec(newHist) but global record of recombinations is (also/usually) stored via newPop().

We are happy to work on a solution, but guidance on the design would be good so we make effective changes @gaynorr ;)

gaynorr commented 2 years ago

This should be fixed in version 1.1.2

gregorgorjanc commented 2 years ago

@gaynorr I am still experiencing issues with "polyploid" functions and pedigree - for example see this session:

Some progress with AlphaSimR 1.1.2 that is addressing https://github.com/gaynorr/AlphaSimR/issues/37 and https://github.com/gaynorr/AlphaSimR/issues/38, but we are not fully there yet :(

> library(SIMplyBee)
> founderGenomes <- quickHaplo(nInd = 10,
+                              nChr = 16,
+                              segSites = 1000)
> SP <- SimParamBee$new(founderGenomes)
> SP$setTrackPed(TRUE)
> base <- newPop(founderGenomes)
> base@id
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"
> base@mother
 [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
> base@father
 [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
> # all zeroes as expected
> 
> drones <- createFounderDrones(pop = base[2], nDronesPerQueen = 5)
> drones@id
[1] "11" "12" "13" "14" "15"
> drones@mother
[1] "2" "2" "2" "2" "2"
> drones@father
[1] "2" "2" "2" "2" "2"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
11      2      2    0
12      2      2    0
13      2      2    0
14      2      2    0
15      2      2    0
> 
> pop <- beeCross(queen = base[1],
+                 drones = drones)
> pop@id
[1] "17"
> pop@mother
[1] "16"
> pop@father
[1] "15"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
11      2      2    0
12      2      2    0
13      2      2    0
14      2      2    0
15      2      2    0
16      1      1    0
17      0      0    0
> 
> library(SIMplyBee)
> founderGenomes <- quickHaplo(nInd = 10,
+                              nChr = 16,
+                              segSites = 1000)
> SP <- SimParamBee$new(founderGenomes)
> SP$setTrackPed(TRUE)
> base <- newPop(founderGenomes)
> base@id
 [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10"
> base@mother
 [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
> base@father
 [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
> # all zeroes as expected
> 
> drones <- createFounderDrones(pop = base[2], nDronesPerQueen = 5)
> drones@id
[1] "11" "12" "13" "14" "15"
> drones@mother
[1] "2" "2" "2" "2" "2"
> drones@father
[1] "2" "2" "2" "2" "2"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
11      2      2    0
12      2      2    0
13      2      2    0
14      2      2    0
15      2      2    0
> # all good!
> 
> queen <- base[1]
> queen@id
[1] "1"
> # Recombination of queen's genomes to generate gametes from the queen
> gametesFromTheQueen <- reduceGenome(pop = queen, nProgeny = 1,
+                                     keepParents = FALSE, simRecomb = TRUE,
+                                     simParam = SP)
> gametesFromTheQueen@id
[1] "16"
> gametesFromTheQueen@mother
[1] "1"
> gametesFromTheQueen@father
[1] "1"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
11      2      2    0
12      2      2    0
13      2      2    0
14      2      2    0
15      2      2    0
16      1      1    0
> # ok pedigree tracking works, but here we actually don't want it (16 here is a gamete/egg not individual) and asked for keepParents = FALSE
> 
> # Drones are already haploid so we just merge both sets of gametes
> pairs <- cbind(gametesFromTheQueen@id,
+                sample(x = drones@id, size = 1, replace = TRUE))
> ret <- mergeGenome(females = gametesFromTheQueen, males = drones,
+                    crossPlan = pairs, simParam = SP)
> 
> ret@id
[1] "17"
> ret@mother
[1] "16"
> ret@father
[1] "15"
> SP$pedigree
   mother father isDH
1       0      0    0
2       0      0    0
3       0      0    0
4       0      0    0
5       0      0    0
6       0      0    0
7       0      0    0
8       0      0    0
9       0      0    0
10      0      0    0
11      2      2    0
12      2      2    0
13      2      2    0
14      2      2    0
15      2      2    0
16      1      1    0
17      0      0    0
> # Hmm, pedigree is not updated when using mergeGenome()
gregorgorjanc commented 2 years ago

I think that this pull requests solves this issue https://github.com/gaynorr/AlphaSimR/pull/47, but we need an additional function to give us a honeybee pedigree we need.