-fixed bug in writeRecords #154

-added `altAddTraitAD` #140

-add miscPop slot to pop #120

-expanded functionality of candidates in selection #138

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: AlphaSimR
 Type: Package
 Title: Breeding Program Simulations
-Version: 1.4.2.9990
-Date: 2023-6-29
+Version: 1.5.0
+Date: 2023-10-30
 Authors@R: c(person("Chris", "Gaynor", email = "gaynor.robert@hotmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")),
   person("Gregor", "Gorjanc", role = "ctb", 

NEWS.md

@@ -4,6 +4,12 @@
 
 *fixed bug in `writePlink` to correctly export map positions in cM
 
+*fixed bug in `writeRecords` due to removed reps slot in pops
+
+*added `altAddTraitAD` for specifying traits with dominance effects using dominance variance and inbreeding depression
+
+*add miscPop slot to class `Pop`
+
 # AlphaSimR 1.4.2
 
 *updated MaCS citation to https site
@@ -34,8 +40,6 @@
 
 # AlphaSimR 1.3.2
 
-*changed name of `MegaPop` to `MultiPop`
-
 *fixed column name bug with multiple traits in `setEBV`
 
 *fixed CTD caused by `runMacs` when too many segSites are requested

R/Class-Pop.R

@@ -379,6 +379,8 @@ isNamedMapPop = function(x) {
 #' population. This list is normally empty and exists solely as an
 #' open slot available for uses to store extra information about
 #' individuals.
+#' @slot miscPop a list of any length containing optional meta data for the 
+#' populations. This list is empty unless information is supplied by the user.
 #'
 #' @export
 setClass("Pop",
@@ -393,7 +395,8 @@ setClass("Pop",
                  ebv="matrix",
                  gxe="list",
                  fixEff="integer",
-                 misc="list"),
+                 misc="list",
+                 miscPop="list"),
          contains="RawPop")
 
 setValidity("Pop",function(object){
@@ -500,6 +503,7 @@ setMethod("[",
             for(chr in 1:x@nChr){
               x@geno[[chr]] = x@geno[[chr]][,,i,drop=FALSE]
             }
+            x@miscPop = list()
             return(x)
           }
 )
@@ -692,7 +696,8 @@ newPop = function(rawPop,simParam=NULL,...){
                ebv=matrix(NA_real_,
                           nrow=rawPop@nInd,
                           ncol=0),
-               misc=vector("list",rawPop@nInd))
+               misc=vector("list",rawPop@nInd),
+               miscPop=list())
   if(simParam$nTraits>=1){
     output = setPheno(output, varE=NULL, reps=1,
                       fixEff=1L, p=NULL, onlyPheno=FALSE,
@@ -879,7 +884,8 @@ newEmptyPop = function(ploidy=2L, simParam=NULL){
                ebv = matrix(NA_real_,
                             nrow=0L,
                             ncol=0L),
-               misc = list())
+               misc = list(),
+               miscPop = list())
   return(output)
 }
 

R/Class-SimParam.R

@@ -609,6 +609,124 @@ SimParam = R6Class(
       }
       invisible(self)
     },
+
+    #' @description 
+    #' An alternative method for adding a trait with additive  and dominance effects 
+    #' to an AlphaSimR simulation. The function attempts to create a trait matching 
+    #' user defined values for number of QTL, inbreeding depression, additive genetic 
+    #' variance and dominance genetic variance.
+    #' 
+    #' @param nQtlPerChr number of QTLs per chromosome. 
+    #' Can be a single value or nChr values.
+    #' @param mean desired mean of the trait
+    #' @param varA desired additive variance
+    #' @param varD desired dominance variance
+    #' @param inbrDepr desired inbreeding depression, see details
+    #' @param limMeanDD limits for meanDD, see details
+    #' @param limVarDD limits for varDD, see details
+    #' @param silent should summary details be printed to the console
+    #' @param force should the check for a running simulation be
+    #' ignored. Only set to TRUE if you know what you are doing.
+    #' @param name optional name for trait
+    #' 
+    #' @details 
+    #' This function will always add a trait to 'SimParam', unless an error occurs 
+    #' with picking QTLs. The resulting trait will always have the desired mean and 
+    #' additive genetic variance. However, it may not have the desired values for 
+    #' inbreeding depression and dominance variance. Thus, it is strongly recommended 
+    #' to check the output printed to the console to determine how close the trait's 
+    #' parameters came to these desired values.
+    #' 
+    #' The mean and additive genetic variance will always be achieved exactly. The 
+    #' function attempts to achieve the desired dominance variance and inbreeding 
+    #' depression while staying within the user supplied constraints for the 
+    #' acceptable range of dominance degree mean and variance. If the desired values
+    #' are not being achieved, the acceptable range need to be increased and/or the 
+    #' number of QTL may need to be increased. There are not limits to setting the 
+    #' range for dominance degree mean and variance, but care should be taken to 
+    #' with regards to the biological feasibility of the limits that are supplied. 
+    #' The default limits were somewhat arbitrarily set, so I make not claim to 
+    #' how reasonable these limits are for routine use.
+    #' 
+    #' Inbreeding depression in this function is defined as the difference in mean 
+    #' genetic value between a population with the same allele frequency as the 
+    #' reference population (population used to initialize SimParam) in 
+    #' Hardy-Weinberg equilibrium compared to a population with the same allele 
+    #' frequency that is fully inbred. This is equivalent to the amount the mean of 
+    #' a population increases when going from an inbreeding coefficient of 1 (fully 
+    #' inbred) to a population with an inbreeding coefficient of 0 (Hardy-Weinberg 
+    #' equilibrium). Note that the sign of the value should (usually) be positive. 
+    #' This corresponds to a detrimental effect of inbreeding when higher values of 
+    #' the trait are considered biologically beneficial.
+    #' 
+    #' Summary information on this trait is printed to the console when silent=FALSE. 
+    #' The summary information reports the inbreeding depression and dominance 
+    #' variance for the population as well as the dominance degree mean and variance 
+    #' applied to the trait.
+    #' 
+    #' @examples
+    #' #Create founder haplotypes
+    #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10)
+    #'
+    #' #Set simulation parameters
+    #' SP = SimParam$new(founderPop)
+    #' SP$altAddTraitAD(nQtlPerChr=10, mean=0, varA=1, varD=0.05, inbrDepr=0.2)
+    altAddTraitAD = function(nQtlPerChr,mean=0,varA=1,varD=0,inbrDepr=0, 
+                             limMeanDD=c(0,1.5),limVarDD=c(0,0.5),
+                             silent=FALSE,force=FALSE,name=NULL){
+      if(!force){
+        private$.isRunning()
+      }
+      if(length(nQtlPerChr)==1){
+        nQtlPerChr = rep(nQtlPerChr,self$nChr)
+      }
+      if(is.null(name)){
+        name = paste0("Trait",self$nTraits+1)
+      }
+
+      # Pick QTL
+      qtlLoci = private$.pickLoci(nQtlPerChr)
+
+      # Create list of arguments for optimization
+      argsList = argAltAD(LociMap = qtlLoci,
+                          Pop = self$founderPop,
+                          mean = mean,
+                          varA = varA,
+                          varD = varD,
+                          inbrDepr = inbrDepr,
+                          nThreads = self$nThreads)
+
+      # Run optim to optimize meanDD and varDD
+      optOut = optim(par = c(mean(limMeanDD), mean(sqrt(limVarDD))),
+                     fn = objAltAD, 
+                     gr = NULL,
+                     method = "L-BFGS-B",
+                     lower = c(limMeanDD[1], sqrt(limVarDD[1])),
+                     upper = c(limMeanDD[2], sqrt(limVarDD[2])),
+                     args = argsList)
+
+      # Finalize creation of trait
+      output = finAltAD(input = optOut$par, args = argsList)
+      trait = new("TraitAD",
+                  qtlLoci,
+                  addEff=c(output$a),
+                  domEff=c(output$d),
+                  intercept=c(output$intercept),
+                  name=name)
+      private$.addTrait(trait,varA,output$varG)
+
+      # Report trait details
+      if(!silent){
+        cat("A new trait called", name, "was added. \n")
+        cat("   varD =", output$varD, "\n")
+        cat("   inbrDepr =", output$inbrDepr, "\n")
+        cat("   meanDD =", output$meanDD, "\n")
+        cat("   varDD =", output$varDD, "\n")
+      }
+
+      invisible(self)
+    },
+
 
     #' @description
     #' Randomly assigns eligible QTLs for one or more additive GxE traits.

R/RcppExports.R

@@ -212,6 +212,18 @@ writeASHaplotypes <- function(g, locations, allLocations, snpchips, names, missi
     invisible(.Call(`_AlphaSimR_writeASHaplotypes`, g, locations, allLocations, snpchips, names, missing, fname))
 }
 
+argAltAD <- function(LociMap, Pop, mean, varA, varD, inbrDepr, nThreads) {
+    .Call(`_AlphaSimR_argAltAD`, LociMap, Pop, mean, varA, varD, inbrDepr, nThreads)
+}
+
+objAltAD <- function(input, args) {
+    .Call(`_AlphaSimR_objAltAD`, input, args)
+}
+
+finAltAD <- function(input, args) {
+    .Call(`_AlphaSimR_finAltAD`, input, args)
+}
+
 calcGenParam <- function(trait, pop, nThreads) {
     .Call(`_AlphaSimR_calcGenParam`, trait, pop, nThreads)
 }

R/mergePops.R

@@ -144,5 +144,6 @@ mergePops = function(popList){
              gxe=gxe,
              pheno=pheno,
              ebv=ebv,
-             misc=misc))
+             misc=misc,
+             miscPop=list()))
 }

R/selection.R

@@ -56,6 +56,27 @@ getResponse = function(pop,trait,use,simParam=NULL,...){
   return(response)
 }
 
+# Converts candidates to a vector of positive numbers
+# for the individuals that are candidates. This function
+# handle indexing by id and negative value indexing
+getCandidates = function(pop, candidates){
+  if(is.character(candidates)){
+    candidates = match(candidates, pop@id)
+    if(any(is.na(candidates))){
+      stop("Trying to select invalid individuals")
+    }
+    if(any(is.null(candidates))){
+      stop("Not valid ids")
+    }
+  }else{
+    if(any(abs(candidates)>pop@nInd)){
+      stop("Trying to select invalid individuals")
+    }
+    candidates = (1:pop@nInd)[candidates]
+  }
+  return(candidates)
+}
+
 # Returns a vector of individuals in a population with the required sex
 checkSexes = function(pop,sex,simParam,...){
   sex = toupper(sex)
@@ -167,6 +188,7 @@ selectInd = function(pop,nInd,trait=1,use="pheno",sex="B",
   }
   eligible = checkSexes(pop=pop,sex=sex,simParam=simParam,...)
   if(!is.null(candidates)){
+    candidates = getCandidates(pop=pop,candidates=candidates)
     eligible = eligible[eligible%in%candidates]
   }
   if(length(eligible)<nInd){
@@ -257,6 +279,7 @@ selectFam = function(pop,nFam,trait=1,use="pheno",sex="B",
   }
   eligible = checkSexes(pop=pop,sex=sex,simParam=simParam,...)
   if(!is.null(candidates)){
+    candidates = getCandidates(pop=pop,candidates=candidates)
     eligible = eligible[eligible%in%candidates]
   }
   allFam = getFam(pop=pop,famType=famType)
@@ -357,6 +380,7 @@ selectWithinFam = function(pop,nInd,trait=1,use="pheno",sex="B",
   }
   eligible = checkSexes(pop=pop,sex=sex,simParam=simParam,...)
   if(!is.null(candidates)){
+    candidates = getCandidates(pop=pop,candidates=candidates)
     eligible = eligible[eligible%in%candidates]
   }
   families = getFam(pop=pop,famType=famType)