Imprinting new functions

R/Class-LociMap.R

@@ -201,6 +201,37 @@ isTraitAI = function(x) {
   return(ret)
 }
 
+#TraitADI----
+#' @title Additive dominance and imprinting trait
+#'
+#' @description Extends \code{\link{TraitAD-class}}
+#' to add imprinting
+#'
+#' @slot impEff imprinting effects
+#'
+#' @export
+setClass("TraitADI",
+         slots=c(impEff="numeric"),
+         contains="TraitAD")
+
+setValidity("TraitADI",function(object){
+  errors = character()
+  if(object@nLoci!=length(object@impEff)){
+    errors = c(errors,"nLoci!=length(impEff)")
+  }
+  if(length(errors)==0){
+    return(TRUE)
+  }else{
+    return(errors)
+  }
+})
+
+# Test if object is of a TraitAI class
+isTraitADI = function(x) {
+  ret = is(x, class2 = "TraitADI")
+  return(ret)
+}
+
 #TraitA2D----
 #' @title Sex specific additive and dominance trait
 #'

R/Class-SimParam.R

@@ -671,10 +671,13 @@ SimParam = R6Class(
         name = paste0("Trait",1:nTraits+self$nTraits)
       }
       stopifnot(length(mean)==length(var),
+                length(meanID)==length(mean),
                 isSymmetric(corA),
                 isSymmetric(corID),
-                length(mean)==nrow(corA),
-                length(mean)==length(name))
+                nrow(corA)==nTraits,
+                nrow(corID)==nTraits,
+                length(varID)==nTraits,
+                length(name)==nTraits)
       qtlLoci = private$.pickLoci(nQtlPerChr)
       addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits,
                           corr=corA,gamma=gamma,shape=shape)
@@ -689,7 +692,6 @@ SimParam = R6Class(
                     name=name[i])
         tmp = calcGenParam(trait, self$founderPop,
                            self$nThreads)
-
         if(useVarA){
           scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1])
         }else{
@@ -707,6 +709,111 @@ SimParam = R6Class(
       invisible(self)
     },
 
+    #' @description
+    #' Randomly assigns eligible QTLs for one or more traits with dominance and imprinting (silencing)
+    #' If simulating more than one trait, all traits will be pleiotropic
+    #' with correlated effects.
+    #'
+    #' @param nQtlPerChr number of QTLs per chromosome. Can be a single value or nChr values.
+    #' @param mean a vector of desired mean genetic values for one or more traits
+    #' @param var a vector of desired genetic variances for one or more traits
+    #' @param meanDD mean dominance degree
+    #' @param varDD variance of dominance degree
+    #' @param meanID mean imprinting degree
+    #' @param varID variance of imprinting degree
+    #' @param corA a matrix of correlations between additive effects
+    #' @param corDD a matrix of correlations between dominance degrees
+    #' @param corID a matrix of correlations between imprinting degrees
+    #' @param useVarA tune according to additive genetic variance if true. If
+    #' FALSE, tuning is performed according to total genetic variance.
+    #' @param gamma should a gamma distribution be used instead of normal
+    #' @param shape the shape parameter for the gamma distribution
+    #'   (the rate/scale parameter of the gamma distribution is accounted
+    #'   for via the desired level of genetic variance, the var argument)
+    #' @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(s)
+    #'
+    #' @examples
+    #' #Create founder haplotypes
+    #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=1)
+    #'
+    #' #Set simulation parameters
+    #' SP = SimParam$new(founderPop)
+    #' \dontshow{SP$nThreads = 1L}
+    #' SP$addTraitADI(1, meanDD=0.5, meanID=0.5)
+    addTraitADI = function(nQtlPerChr,mean=0,var=1,meanDD=0,varDD=0,meanID=0,
+                          varID=0,corA=NULL,corDD=NULL,corID=NULL,useVarA=TRUE,
+                          gamma=FALSE,shape=1,force=FALSE,name=NULL){
+      if(self$founderPop@ploidy > 2){
+        stop("ERROR: Imprinting is not yet implemented for polyploids!")
+      }
+      if(!force){
+        private$.isRunning()
+      }
+      if(length(nQtlPerChr)==1){
+        nQtlPerChr = rep(nQtlPerChr,self$nChr)
+      }
+      nTraits = length(mean)
+      if(length(meanDD)==1) meanDD = rep(meanDD,nTraits)
+      if(length(meanID)==1) meanID = rep(meanID,nTraits)
+      if(length(varDD)==1) varDD = rep(varDD,nTraits)
+      if(length(varID)==1) varID = rep(varID,nTraits)
+      if(length(gamma)==1) gamma = rep(gamma,nTraits)
+      if(length(shape)==1) shape = rep(shape,nTraits)
+      if(is.null(corA)) corA=diag(nTraits)
+      if(is.null(corDD)) corDD=diag(nTraits)
+      if(is.null(corID)) corID=diag(nTraits)
+      if(is.null(name)){
+        name = paste0("Trait",1:nTraits+self$nTraits)
+      }
+      stopifnot(length(mean)==length(var),
+                length(meanDD)==length(mean),
+                length(meanID)==length(mean),
+                isSymmetric(corA),
+                isSymmetric(corDD),
+                isSymmetric(corID),
+                nrow(corA)==nTraits,
+                nrow(corDD)==nTraits,
+                nrow(corID)==nTraits,
+                length(varDD)==nTraits,
+                length(varID)==nTraits,
+                length(name)==nTraits)
+      qtlLoci = private$.pickLoci(nQtlPerChr)
+      addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits,
+                          corr=corA,gamma=gamma,shape=shape)
+      domEff = sampDomEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff,
+                          corDD=corDD,meanDD=meanDD,varDD=varDD)
+      impEff = sampImpEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff,
+                          corID=corID,meanID=meanID,varID=varID)
+      for(i in 1:nTraits){
+        trait = new("TraitADI",
+                    qtlLoci,
+                    addEff=addEff[,i],
+                    domEff=domEff[,i],
+                    impEff=impEff[,i],
+                    intercept=0,
+                    name=name[i])
+        tmp = calcGenParam(trait, self$founderPop,
+                           self$nThreads)
+        if(useVarA){
+          scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1])
+        }else{
+          scale = sqrt(var[i])/sqrt(popVar(tmp$gv)[1])
+        }
+        trait@addEff = trait@addEff*scale
+        trait@domEff = trait@domEff*scale
+        trait@impEff = trait@impEff*scale
+        trait@intercept = mean[i]-mean(tmp$gv*scale)
+        if(useVarA){
+          private$.addTrait(trait,var[i],popVar(tmp$gv*scale)[1])
+        }else{
+          private$.addTrait(trait,popVar(tmp$bv*scale)[1],var[i])
+        }
+      }
+      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
@@ -1527,13 +1634,14 @@ SimParam = R6Class(
     #' formatting the trait as a \code{\link{LociMap-class}}.
     #' The formatting is performed automatically for the user,
     #' with more user friendly data.frames or matrices taken as
-    #' inputs. This function only works for A and AD trait types.
+    #' inputs. This function only works for A, AD, AI and ADI trait types.
     #'
     #' @param markerNames a vector of names for the QTL
     #' @param addEff a matrix of additive effects (nLoci x nTraits).
     #' Alternatively, a vector of length nLoci can be supplied for
     #' a single trait.
     #' @param domEff optional dominance effects for each locus
+    #' @param impEff optional imprinting effects for each locus
     #' @param intercept optional intercepts for each trait
     #' @param name optional name(s) for the trait(s)
     #' @param varE default error variance for phenotype, optional
@@ -1542,6 +1650,7 @@ SimParam = R6Class(
     importTrait = function(markerNames,
                            addEff,
                            domEff=NULL,
+                           impEff=NULL,
                            intercept=NULL,
                            name=NULL,
                            varE=NULL,
@@ -1562,6 +1671,16 @@ SimParam = R6Class(
         stopifnot(nrow(addEff)==nrow(domEff),
                   ncol(addEff)==ncol(domEff))
       }
+      if(is.null(impEff)){
+        useImp = FALSE
+      }else{
+        useImp = TRUE
+        impEff = as.matrix(impEff)
+        stopifnot(nrow(addEff)==nrow(impEff),
+                  ncol(addEff)==ncol(impEff),
+                  nrow(domEff)==nrow(impEff),
+                  ncol(domEff)==ncol(impEff))
+      }
 
       # Prepare the intercept
       if(is.null(intercept)){
@@ -1593,16 +1712,16 @@ SimParam = R6Class(
       # Create trait variables
       lociPerChr = integer(self$nChr)
       lociLoc = vector("list", self$nChr)
-      addEffList = domEffList = vector("list", nTraits)
+      addEffList = domEffList = impEffList = vector("list", nTraits)
       for(i in 1:nTraits){
-        addEffList[[i]] = domEffList[[i]] = vector("list", self$nChr)
+        addEffList[[i]] = domEffList[[i]] = impEffList[[i]] = vector("list", self$nChr)
       }
 
       # Loop through chromosomes
       for(i in 1:self$nChr){
         # Working on trait 1
         # Initialize variables
-        addEffList[[1]][[i]] = domEffList[[1]][[i]] = numeric()
+        addEffList[[1]][[i]] = domEffList[[1]][[i]] = impEffList[[1]][[i]] = numeric()
         lociLoc[[i]] = integer()
 
         # Find matches if they exist
@@ -1615,17 +1734,23 @@ SimParam = R6Class(
           if(useDom){
             domEffList[[1]][[i]] = domEff[na.omit(take),1]
           }
+          if(useImp){
+            impEffList[[1]][[i]] = impEff[na.omit(take),1]
+          }
         }
 
         # Work on additional traits?
         if(nTraits>1){
           for(j in 2:nTraits){
-            addEffList[[j]][[i]] = domEffList[[j]][[i]] = numeric()
+            addEffList[[j]][[i]] = domEffList[[j]][[i]] = impEffList[[j]][[i]] = numeric()
             if(lociPerChr[i]>0L){
               addEffList[[j]][[i]] = addEff[na.omit(take),j]
               if(useDom){
                 domEffList[[j]][[i]] = domEff[na.omit(take),j]
               }
+              if(useImp){
+                impEffList[[j]][[i]] = impEff[na.omit(take),j]
+              }
             }
           }
         }
@@ -1639,6 +1764,18 @@ SimParam = R6Class(
         addEff = unlist(addEffList[[i]])
         if(useDom){
           domEff = unlist(domEffList[[i]])
+          if(useImp){
+            impEff = unlist(impEffList[[i]])
+            trait = new("TraitADI",
+                        addEff=addEff,
+                        domEff=domEff,
+                        impEff=impEff,
+                        intercept=intercept[i],
+                        nLoci=nLoci,
+                        lociPerChr=lociPerChr,
+                        lociLoc=lociLoc,
+                        name=name[i])
+          }else{
           trait = new("TraitAD",
                       addEff=addEff,
                       domEff=domEff,
@@ -1647,14 +1784,27 @@ SimParam = R6Class(
                       lociPerChr=lociPerChr,
                       lociLoc=lociLoc,
                       name=name[i])
+          }
         }else{
+          if(useImp){
+            impEff = unlist(impEffList[[i]])
+            trait = new("TraitAI",
+                        addEff=addEff,
+                        impEff=impEff,
+                        intercept=intercept[i],
+                        nLoci=nLoci,
+                        lociPerChr=lociPerChr,
+                        lociLoc=lociLoc,
+                        name=name[i])
+          }else{
           trait = new("TraitA",
                       addEff=addEff,
                       intercept=intercept[i],
                       nLoci=nLoci,
                       lociPerChr=lociPerChr,
                       lociLoc=lociLoc,
                       name=name[i])
+          }
         }
 
         # Add trait to simParam

R/importData.R

@@ -276,7 +276,7 @@ importHaplo = function(haplo, genMap, ploidy=2L, ped=NULL){
       genMap[[i]] = genMap[[i]] - genMap[[i]]-genMap[[i]][1]
       take = na.omit(take)
     }
-    haplotypes[[i]] = haplo[,take]
+    haplotypes[[i]] = as.matrix(haplo[,take])
   }
 
   founderPop = newMapPop(genMap=genMap,