Simulation_Crossbreeding/CreateFounderPops.R at main · SeyiV/Simulation_Crossbreeding · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

# Founder populations

# This part of simulation with runMacs() and repeat {} is quite slow!
# There are two options one can take regarding runMacs():
# 1) Save the simulated founder genomes and load them in the next run
#    (only for testing purposes - we want variation between replicates)
# 2) Use the quickHaplo() function to speed up the simulation
#    (only for testing purposes - it does not generate properly structured genomes)

# ---- Simulate Founders genome ----

# We can speed this up by saving founder genomes and loading them instead
# Beware, we must however resimulate founder genomes for each replicate!
FounderGenomes <- runMacs(
  nInd = NFounders,
  nChr = nChr,
  segSites = nSNP + nQTL,
  split = nSplit,
  species = "CATTLE")
if (FALSE) {
  FounderGenomes <- quickHaplo(
    nInd = NFounders,
    nChr = nChr,
    segSites = nSNP + nQTL)
}
# save(FounderGenomes, file = "FounderGenomes.RData")
# load(file = "FounderGenomes.RData")

# ---- Traits with additive and dominance genetic effects ----

# Using repeat to get a desired setting with additive and dominance effects
# (ongoing work by AlphaSimR developers will remove the need for such an approach;
#  it will enable specifying the desired level of dominance variance and inbreeding
#  depression)
repeat {
  SP <- SimParam$new(FounderGenomes)
  SP$addSnpChip(nSnpPerChr=nSNP)
  SP$addTraitAD(nQTL,
                mean = PhenoMean,
                var = AdditVar,
                corA = TraitCor,
                meanDD = DomMeanDD,
                varDD = DomVarDD,
                name = c("BodyWeight_local", "TickCount_local",
                         "BodyWeight_exotic", "TickCount_exotic"))

  # Systematically assign sexes to individuals
  SP$setSexes("yes_sys")

  # Generate initial founder population
  Founders <- newPop(FounderGenomes)
  Founders@misc <- list(gen = rep(0, times = nInd(Founders)))

  # Split the initial founder into local and exotic founders
  LocalFounders <- Founders[1:2500]
  LocalFounders@sex <- sample(rep(c("F", "M"), c(2000, 500)), 2500, replace = FALSE)
  LocalFounders <- setPheno(pop = LocalFounders, h2 = h2)

  ExoticFounders <- Founders[2501:5000]
  ExoticFounders@sex <- sample(rep(c("F", "M"), c(2000, 500)), 2500, replace = FALSE)
  ExoticFounders <- setPheno(pop = ExoticFounders, h2 = h2)

  # Ensure that the expected heritabilities (±10%) are obtained in the founder populations
  h2_local <- diag(varA(LocalFounders)/varP(LocalFounders))
  h2_exotic <- diag(varA(ExoticFounders)/varP(ExoticFounders))
  print(h2_local)
  print(h2_exotic)
  if (h2_local[1]  > 0.27 & h2_local[1]  < 0.33 &
      h2_local[2]  > 0.09 & h2_local[2]  < 0.11 &
      h2_exotic[1] > 0.27 & h2_exotic[1] < 0.33 &
      h2_exotic[2] > 0.09 & h2_exotic[2] < 0.11 &
      h2_exotic[3] > 0.27 & h2_exotic[3] < 0.33 &
      h2_exotic[4] > 0.09 & h2_exotic[4] < 0.11) {
    break
  }
}
cat("Ratio of dominance variance to phenotypic variance\n")
print(diag(varD(LocalFounders)/varP(LocalFounders)))
print(diag(varD(LocalFounders)/varP(LocalFounders)))

# ---- Evaluate genetic distance between local and exotic founders ----

Fst <- calcFst(LocalFounders, ExoticFounders, Founders)
HetFounders <- calcHet(Founders)
HetLocalFounders <- calcHet(LocalFounders)
HetExoticFounders <- calcHet(ExoticFounders)