diff --git a/DESCRIPTION b/DESCRIPTION index e2cd7099..eb3f7bf2 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: AlphaSimR Type: Package Title: Breeding Program Simulations -Version: 1.5.3 -Date: 2023-11-30 +Version: 1.6.0 +Date: 2024-08-15 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", @@ -39,7 +39,7 @@ Depends: R (>= 4.0.0), methods, R6 Imports: Rcpp (>= 0.12.7), Rdpack RdMacros: Rdpack LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.2 Suggests: knitr, rmarkdown, testthat VignetteBuilder: knitr NeedsCompilation: true diff --git a/NAMESPACE b/NAMESPACE index 6c76c2dc..144729ec 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,7 +28,6 @@ export(genicVarAA) export(genicVarD) export(genicVarG) export(getGenMap) -export(getMisc) export(getNumThreads) export(getPed) export(getQtlMap) @@ -49,6 +48,7 @@ export(isRawPop) export(makeCross) export(makeCross2) export(makeDH) +export(meanEBV) export(meanG) export(meanP) export(mergeGenome) @@ -89,7 +89,6 @@ export(selectWithinFam) export(self) export(setEBV) export(setMarkerHaplo) -export(setMisc) export(setPheno) export(setPhenoGCA) export(setPhenoProgTest) @@ -108,6 +107,7 @@ export(usefulness) export(varA) export(varAA) export(varD) +export(varEBV) export(varG) export(varP) export(writePlink) diff --git a/NEWS.md b/NEWS.md index 9a73d75a..00f593c0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +# AlphaSimR 1.6.0 + +*exported `meanEBV` and added `varEBV` to complement `meanP`/`varP` and `meanG`/`varG` + +*Changed all parameters of the CATTLE demographic model to exactly match Macleod et al. (2013) - specifically reducing the mutation rate from 2.5e-8 (from human literature) to 1.2e-8 (used in Macleod et al., 2013) and recombination rate from 1e-8 (generic) to 9.26e-9 (used in Macleod et al., 2013). These changes will reduce number of segregating sites to ~240K per chromosome for 100 samples and will run faster. + +*changed misc slot in Pop class from a list organised as ind x nodes to to a list organised as nodes x ind (this simplified code and increased speed) + +*removed `setMisc` and `getMisc` because the new misc slot structure makes it easy to set and get misc components with base R code + +*added `length` method for Pop class that returns number of individuals (like `nInd`) + +*added `length` method for MultiPop class that returns number of populations + +*fixed bug in quadrivalent pairing resulting in distribution of double reductions not respecting the centromere + # AlphaSimR 1.5.3 *fixed bug in `SimParam$restrSegSites` with excluding sites at end of chromosome diff --git a/R/Class-Pop.R b/R/Class-Pop.R index 84ac2105..c7862046 100644 --- a/R/Class-Pop.R +++ b/R/Class-Pop.R @@ -148,7 +148,7 @@ setValidity("MapPop",function(object){ if(object@nChr!=length(object@genMap)){ errors = c(errors,"nInd!=length(id)") } - for(i in 1:object@nChr){ + for(i in seq_len(object@nChr)){ if(object@nLoci[i]!=length(object@genMap[[i]])){ errors = c(errors, paste0("nLoci[",i,"]!=length(genMap[[",i,"]]")) @@ -168,7 +168,7 @@ setMethod("[", if(any(abs(i)>x@nInd)){ stop("Trying to select invalid individuals") } - for(chr in 1:x@nChr){ + for(chr in seq_len(x@nChr)){ x@geno[[chr]] = x@geno[[chr]][,,i,drop=FALSE] } x@nInd = dim(x@geno[[1]])[3] @@ -375,13 +375,19 @@ isNamedMapPop = function(x) { #' @slot gxe list containing GxE slopes for GxE traits #' @slot fixEff a fixed effect relating to the phenotype. #' Used by genomic selection models but otherwise ignored. -#' @slot misc a list whose elements correspond to individuals in the -#' population. This list is normally empty and exists solely as an +#' @slot misc a list whose elements correspond to additional miscellaneous +#' nodes with the items for individuals in the population (see example in +#' \code{\link{newPop}}). +#' 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 -#' population. This list is empty unless information is supplied by the user. -#' Note that the list is emptied every time the population is subsetted. +#' @slot miscPop a list of any length containing optional meta data for the +#' population (see example in \code{\link{newPop}}). +#' This list is empty unless information is supplied by the user. +#' Note that the list is emptied every time the population is subsetted or +#' combined because the meta data for old population might not be valid anymore. +#' +#' @seealso \code{\link{newPop}}, \code{\link{newEmptyPop}}, \code{\link{resetPop}} #' #' @export setClass("Pop", @@ -456,8 +462,8 @@ setValidity("Pop",function(object){ if(object@nInd!=length(object@fixEff)){ errors = c(errors,"nInd!=length(fixEff)") } - if(object@nInd!=length(object@misc)){ - errors = c(errors,"nInd!=length(misc)") + if(any(object@nInd!=sapply(object@misc, length))){ + errors = c(errors,"any(nInd!=sapply(misc, length))") } if(length(errors)==0){ return(TRUE) @@ -488,7 +494,8 @@ setMethod("[", x@mother = x@mother[i] x@father = x@father[i] x@fixEff = x@fixEff[i] - x@misc = x@misc[i] + x@misc = lapply(x@misc, FUN = function(z) z[i]) + x@miscPop = list() x@gv = x@gv[i,,drop=FALSE] x@pheno = x@pheno[i,,drop=FALSE] x@ebv = x@ebv[i,,drop=FALSE] @@ -504,7 +511,6 @@ setMethod("[", for(chr in 1:x@nChr){ x@geno[[chr]] = x@geno[[chr]][,,i,drop=FALSE] } - x@miscPop = list() return(x) } ) @@ -534,12 +540,20 @@ setMethod("show", } ) +#' @describeIn Pop Number of individuals in Pop (the same as nInd()) +setMethod("length", + signature(x = "Pop"), + function (x){ + return(x@nInd) + } +) + #' @title Create new population #' #' @description #' Creates an initial \code{\link{Pop-class}} from an object of #' \code{\link{MapPop-class}} or \code{\link{NamedMapPop-class}}. -#' The function is intended for us with output from functions such +#' The function is intended for use with output from functions such #' as \code{\link{runMacs}}, \code{\link{newMapPop}}, or #' \code{\link{quickHaplo}}. #' @@ -550,6 +564,14 @@ setMethod("show", #' #' @return Returns an object of \code{\link{Pop-class}} #' +#' @details Note that \code{newPop} takes genomes from the +#' \code{rawPop} and uses them without recombination! Hence, if you +#' call \code{newPop(rawPop = founderGenomes)} twice, you will get +#' two sets of individuals with different id but the same genomes. +#' To get genetically different sets of individuals you can subset the +#' \code{rawPop} input, say first half for one set and the second half +#' for the other set. +#' #' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) @@ -562,6 +584,13 @@ setMethod("show", #' pop = newPop(founderPop, simParam=SP) #' isPop(pop) #' +#' #Misc +#' pop@misc$tmp1 = rnorm(n=2) +#' pop@misc$tmp2 = rnorm(n=2) +#' +#' #MiscPop +#' pop@miscPop$tmp1 = sum(pop@misc$tmp1) +#' pop@miscPop$tmp2 = sum(pop@misc$tmp2) #' @export newPop = function(rawPop,simParam=NULL,...){ if(is.null(simParam)){ @@ -602,7 +631,7 @@ newPop = function(rawPop,simParam=NULL,...){ stopifnot(sapply(simParam$genMap,length)==rawPop@nLoci) lastId = simParam$lastId - iid = (1:rawPop@nInd) + lastId + iid = seq_len(rawPop@nInd) + lastId lastId = max(iid) if(is.null(id)){ @@ -666,7 +695,7 @@ newPop = function(rawPop,simParam=NULL,...){ pheno = gv if(simParam$nTraits>=1){ - for(i in 1:simParam$nTraits){ + for(i in seq_len(simParam$nTraits)){ tmp = getGv(simParam$traits[[i]], rawPop, simParam$nThreads) gv[,i] = tmp[[1]] @@ -690,15 +719,15 @@ newPop = function(rawPop,simParam=NULL,...){ mother=mother, father=father, fixEff=rep(1L,rawPop@nInd), + misc=list(), + miscPop=list(), nTraits=simParam$nTraits, gv=gv, gxe=gxe, pheno=pheno, ebv=matrix(NA_real_, nrow=rawPop@nInd, - ncol=0), - misc=vector("list",rawPop@nInd), - miscPop=list()) + ncol=0)) if(simParam$nTraits>=1){ output = setPheno(output, varE=NULL, reps=1, fixEff=1L, p=NULL, onlyPheno=FALSE, @@ -752,10 +781,10 @@ resetPop = function(pop,simParam=NULL){ simParam = get("SP",envir=.GlobalEnv) } pop@nTraits = simParam$nTraits - + # Extract names to add back at the end traitNames = colnames(pop@gv) - + # Create empty slots for traits pop@pheno = matrix(NA_real_, nrow=pop@nInd, @@ -769,19 +798,17 @@ resetPop = function(pop,simParam=NULL){ pop@fixEff = rep(1L,pop@nInd) # Calculate genetic values - if(simParam$nTraits>=1){ - for(i in 1:simParam$nTraits){ - tmp = getGv(simParam$traits[[i]],pop,simParam$nThreads) - pop@gv[,i] = tmp[[1]] - if(length(tmp)>1){ - pop@gxe[[i]] = tmp[[2]] - } + for(i in seq_len(simParam$nTraits)){ + tmp = getGv(simParam$traits[[i]],pop,simParam$nThreads) + pop@gv[,i] = tmp[[1]] + if(length(tmp)>1){ + pop@gxe[[i]] = tmp[[2]] } } # Add back trait names colnames(pop@pheno) = colnames(pop@gv) = traitNames - + return(pop) } @@ -848,20 +875,18 @@ newEmptyPop = function(ploidy=2L, simParam=NULL){ ncol = simParam$nTraits) traitNames = character(simParam$nTraits) - - if(simParam$nTraits > 0L){ - # Get trait names - for(i in 1:simParam$nTraits){ - traitNames[i] = simParam$traits[[i]]@name - } + + # Get trait names + for(i in seq_len(simParam$nTraits)){ + traitNames[i] = simParam$traits[[i]]@name } - + colnames(traitMat) = traitNames # Create empty geno list nLoci = unname(sapply(simParam$genMap, length)) geno = vector("list", simParam$nChr) - for(i in 1:simParam$nChr){ + for(i in seq_len(simParam$nChr)){ DIM1 = nLoci[i]%/%8L + (nLoci[i]%%8L > 0L) geno[[i]] = array(as.raw(0), dim=c(DIM1, ploidy, 0)) } @@ -878,15 +903,15 @@ newEmptyPop = function(ploidy=2L, simParam=NULL){ mother = character(), father = character(), fixEff = integer(), + misc = list(), + miscPop = list(), nTraits = simParam$nTraits, gv = traitMat, gxe = vector("list", simParam$nTraits), pheno = traitMat, ebv = matrix(NA_real_, nrow=0L, - ncol=0L), - misc = list(), - miscPop = list()) + ncol=0L)) return(output) } @@ -913,7 +938,7 @@ setClass("MultiPop", setValidity("MultiPop",function(object){ errors = character() # Check that all populations are valid - for(i in 1:length(object@pops)){ + for(i in seq_len(length(object@pops))){ if(!validObject(object@pops[[i]]) & (is(object@pops[[i]], "Pop") | is(object@pops[[i]],"MultiPop"))){ @@ -964,6 +989,15 @@ setMethod("c", } ) +#' @describeIn MultiPop Number of pops in MultiPop +setMethod("length", + signature(x = "MultiPop"), + function (x){ + n = length(x@pops) + return(n) + } +) + #' @title Create new Multi Population #' #' @description diff --git a/R/Class-SimParam.R b/R/Class-SimParam.R index 33549cad..8fd0031b 100644 --- a/R/Class-SimParam.R +++ b/R/Class-SimParam.R @@ -219,7 +219,7 @@ SimParam = R6Class( # Make exclusions restr = self$invalidQtl - for(i in 1:self$nChr){ + for(i in seq_len(self$nChr)){ restr[[i]] = sort(union(restr[[i]], matchList[[i]])) } self$invalidQtl = restr @@ -242,7 +242,7 @@ SimParam = R6Class( # Make exclusions restr = self$invalidSnp - for(i in 1:self$nChr){ + for(i in seq_len(self$nChr)){ restr[[i]] = sort(union(restr[[i]], matchList[[i]])) } self$invalidSnp = restr @@ -396,7 +396,7 @@ SimParam = R6Class( lociLoc = vector("list", length(genMap)) # Loop through chromosomes - for(i in 1:length(genMap)){ + for(i in seq_len(length(genMap))){ # Initialize lociLoc lociLoc[[i]] = integer() @@ -450,11 +450,11 @@ SimParam = R6Class( }) lociLoc = do.call("c",lociLoc) - for (i in 1:nrow(structure)){ + for (i in seq_len(nrow(structure))){ snps = lociLoc[structure[i,]] start = 1 numChr = numeric(length(nSnpPerChr)) - for (j in 1:length(nSnpPerChr)){ + for (j in seq_len(length(nSnpPerChr))){ end = start + nSnpPerChr[j] - 1 numChr[j] = sum(structure[i,start:end]) start = end + 1 @@ -517,7 +517,7 @@ SimParam = R6Class( qtlLoci = private$.pickLoci(nQtlPerChr) addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corA,gamma=gamma,shape=shape) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitA", qtlLoci, addEff=addEff[,i], @@ -592,7 +592,7 @@ SimParam = R6Class( corr=corA,gamma=gamma,shape=shape) domEff = sampDomEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, corDD=corDD,meanDD=meanDD,varDD=varDD) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitAD", qtlLoci, addEff=addEff[,i], @@ -796,7 +796,7 @@ SimParam = R6Class( corr=corA,gamma=gamma,shape=shape) gxeEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corGxE,gamma=FALSE,shape=NULL) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitA", qtlLoci, addEff=addEff[,i], @@ -906,7 +906,7 @@ SimParam = R6Class( corDD=corDD,meanDD=meanDD,varDD=varDD) gxeEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corGxE,gamma=FALSE,shape=NULL) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitAD", qtlLoci, addEff=addEff[,i], @@ -1019,7 +1019,7 @@ SimParam = R6Class( corr=corA,gamma=gamma,shape=shape, relVar=relAA) E = matrix(sample.int(sum(nQtlPerChr),sum(nQtlPerChr)),ncol=2) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitAE", qtlLoci, addEff=addEff[,i], @@ -1119,7 +1119,7 @@ SimParam = R6Class( corr=corA,gamma=gamma,shape=shape, relVar=relAA) E = matrix(sample.int(sum(nQtlPerChr),sum(nQtlPerChr)),ncol=2) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitADE", qtlLoci, addEff=addEff[,i], @@ -1222,7 +1222,7 @@ SimParam = R6Class( E = matrix(sample.int(sum(nQtlPerChr),sum(nQtlPerChr)),ncol=2) gxeEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corGxE,gamma=FALSE,shape=NULL) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitAE", qtlLoci, addEff=addEff[,i], @@ -1357,7 +1357,7 @@ SimParam = R6Class( E = matrix(sample.int(sum(nQtlPerChr),sum(nQtlPerChr)),ncol=2) gxeEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corGxE,gamma=FALSE,shape=NULL) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = new("TraitADE", qtlLoci, addEff=addEff[,i], @@ -1505,12 +1505,12 @@ SimParam = R6Class( lociPerChr = integer(self$nChr) lociLoc = vector("list", self$nChr) addEffList = domEffList = vector("list", nTraits) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ addEffList[[i]] = domEffList[[i]] = vector("list", self$nChr) } # Loop through chromosomes - for(i in 1:self$nChr){ + for(i in seq_len(self$nChr)){ # Working on trait 1 # Initialize variables addEffList[[1]][[i]] = domEffList[[1]][[i]] = numeric() @@ -1546,7 +1546,7 @@ SimParam = R6Class( nLoci = sum(lociPerChr) # Create Trait(s) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ addEff = unlist(addEffList[[i]]) if(useDom){ domEff = unlist(domEffList[[i]]) @@ -1597,12 +1597,10 @@ SimParam = R6Class( private$.varA[traitPos] = popVar(tmp$bv)[1] private$.varG[traitPos] = popVar(tmp$gv)[1] if(is.matrix(private$.varE)){ - private$.varE[traitPos,] = 0 - private$.varE[,traitPos] = 0 - private$.varE[traitPos,traitPos] = varE - }else{ - private$.varE[traitPos] = varE + warning("Error correlations have been removed, use setVarE to reinstate") + private$.varE = diag(private$.varE) } + private$.varE[traitPos] = varE invisible(self) }, @@ -1662,7 +1660,7 @@ SimParam = R6Class( all(private$.varG>0), all(private$.varA>0)) varE = numeric(self$nTraits) - for(i in 1:length(h2)){ + for(i in seq_len(length(h2))){ tmp = private$.varA[i]/h2[i]-private$.varG[i] if(tmp<0){ stop(paste0("h2=",h2[i]," is not possible for trait ",i)) @@ -1673,7 +1671,7 @@ SimParam = R6Class( }else if(!is.null(H2)){ stopifnot(length(H2)==self$nTraits) varE = numeric(self$nTraits) - for(i in 1:length(H2)){ + for(i in seq_len(length(H2))){ tmp = private$.varG[i]/H2[i]-private$.varG[i] varE[i] = tmp } @@ -1728,7 +1726,8 @@ SimParam = R6Class( warning("This function has been deprecated. Use simParam$setVarE instead.") stopifnot(isSymmetric(corE), nrow(corE)==self$nTraits, - length(private$.varE)==self$nTraits) + length(private$.varE)==self$nTraits, + !any(is.na(private$.varE))) if(is.matrix(private$.varE)){ varE = diag(private$.varE) }else{ @@ -1783,7 +1782,7 @@ SimParam = R6Class( length(var)==self$nTraits, length(varEnv)==self$nTraits, length(varGxE)==self$nTraits) - for(i in 1:self$nTraits){ + for(i in seq_len(self$nTraits)){ trait = private$.traits[[i]] trait@intercept = 0 tmp = calcGenParam(trait, @@ -1994,12 +1993,12 @@ SimParam = R6Class( newRecHist = vector("list",nNewInd) tmpLastHaplo = private$.lastHaplo if(all(isDH==1L)){ - for(i in 1:nNewInd){ + for(i in seq_len(nNewInd)){ tmpLastHaplo = tmpLastHaplo + 1L newRecHist[[i]] = rep(tmpLastHaplo, ploidy) } }else{ - for(i in 1:nNewInd){ + for(i in seq_len(nNewInd)){ newRecHist[[i]] = (tmpLastHaplo+1L):(tmpLastHaplo+ploidy) tmpLastHaplo = tmpLastHaplo + ploidy } @@ -2169,7 +2168,7 @@ SimParam = R6Class( # Identify potential sites pot = vector('list', self$nChr) - for(i in 1:self$nChr){ + for(i in seq_len(self$nChr)){ pot[[i]] = setdiff(1:private$.segSites[i], restr[[i]]) } @@ -2178,7 +2177,7 @@ SimParam = R6Class( if(is.null(refPop)){ refPop = self$founderPop } - for(chr in 1:self$nChr){ + for(chr in seq_len(self$nChr)){ q = calcChrFreq(refPop@geno[[chr]]) q = 0.5-abs(q-0.5) #Convert to minor allele frequency tmp = which(q>=minFreq) @@ -2267,7 +2266,7 @@ SimParam = R6Class( if(length(value)!=self$nTraits){ stop("length of traitNames vector must equal ",self$nTraits) } - for(i in 1:self$nTraits){ + for(i in seq_len(self$nTraits)){ private$.traits[[i]]@name = value[i] } } @@ -2285,7 +2284,7 @@ SimParam = R6Class( if(length(value)!=self$nSnpChips){ stop("length of snpChipNames vector must equal ",self$nSnpChips) } - for(i in 1:self$nSnpChips){ + for(i in seq_len(self$nSnpChips)){ self$snpChips[[i]]@name = value[i] } } @@ -2360,7 +2359,7 @@ SimParam = R6Class( if(missing(value)){ if(private$.sepMap){ genMap = vector("list",self$nChr) - for(i in 1:self$nChr){ + for(i in seq_len(self$nChr)){ genMap[[i]] = (private$.femaleMap[[i]]+ private$.maleMap[[i]])/2 } diff --git a/R/GS.R b/R/GS.R index 6ace27fb..0f65ac2c 100644 --- a/R/GS.R +++ b/R/GS.R @@ -17,68 +17,68 @@ convertTraitsToNames = function(traits, simParam){ #' @title Fast RR-BLUP #' #' @description -#' Solves an RR-BLUP model for genomic predictions given known variance -#' components. This implementation is meant as a fast and low memory -#' alternative to \code{\link{RRBLUP}} or \code{\link{RRBLUP2}}. Unlike -#' the those functions, the fastRRBLUP does not fit fixed effects (other -#' than the intercept) or account for unequal replication. +#' Solves an RR-BLUP model for genomic predictions given known variance +#' components. This implementation is meant as a fast and low memory +#' alternative to \code{\link{RRBLUP}} or \code{\link{RRBLUP2}}. Unlike +#' the those functions, the fastRRBLUP does not fit fixed effects (other +#' than the intercept) or account for unequal replication. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population -#' @param traits an integer indicating the trait to model, a trait name, -#' or a function of the traits returning a single value. Only univariate models +#' @param traits an integer indicating the trait to model, a trait name, +#' or a function of the traits returning a single value. Only univariate models #' are supported. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. -#' @param maxIter maximum number of iterations. -#' @param Vu marker effect variance. If value is NULL, a +#' @param maxIter maximum number of iterations. +#' @param Vu marker effect variance. If value is NULL, a #' reasonable value is chosen automatically. -#' @param Ve error variance. If value is NULL, a +#' @param Ve error variance. If value is NULL, a #' reasonable value is chosen automatically. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = fastRRBLUP(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=1000, Vu=NULL, Ve=NULL, +fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=1000, Vu=NULL, Ve=NULL, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + #fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -88,7 +88,7 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + # Sort out Vu and Ve if(is.function(traits)){ if(is.null(Vu)){ @@ -112,12 +112,12 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, } } } - + #Fit model ans = callFastRRBLUP(y,pop@geno,lociPerChr, lociLoc,Vu,Ve,maxIter, simParam$nThreads) - + bv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -125,7 +125,7 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$beta), name=paste0("est_BV_",traits)) - + gv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -133,7 +133,7 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + output = new("RRsol", bv = list(bv), gv = list(gv), @@ -141,7 +141,7 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, male = as.list(NULL), Vu = as.matrix(Vu), Ve = as.matrix(Ve)) - + return(output) } @@ -153,57 +153,57 @@ fastRRBLUP = function(pop, traits=1, use="pheno", snpChip=1, #' Fits an RR-BLUP model for genomic predictions. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population -#' @param traits an integer indicating the trait or traits to model, a vector of trait names, +#' @param traits an integer indicating the trait or traits to model, a vector of trait names, #' or a function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. -#' @param maxIter maximum number of iterations. Only used +#' @param maxIter maximum number of iterations. Only used #' when number of traits is greater than 1. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=1000L, +RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=1000L, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -213,7 +213,7 @@ RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + #Fit model if(ncol(y)>1){ ans = callRRBLUP_MV(y, fixEff, pop@geno, lociPerChr, @@ -222,12 +222,12 @@ RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, ans = callRRBLUP(y, fixEff, pop@geno, lociPerChr, lociLoc, simParam$nThreads) } - + markerEff=ans$u - + bv = gv = vector("list",ncol(y)) - for(i in 1:ncol(y)){ + for(i in seq_len(ncol(y))){ bv[[i]] = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -235,7 +235,7 @@ RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, addEff=ans$alpha[,i], intercept=ans$beta[i], name=paste0("est_BV_",traits[i])) - + gv[[i]] = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -244,7 +244,7 @@ RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, intercept=ans$mu[i], name=paste0("est_GV_",traits[i])) } - + output = new("RRsol", bv = bv, gv = gv, @@ -252,107 +252,107 @@ RRBLUP = function(pop, traits=1, use="pheno", snpChip=1, male = as.list(NULL), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } #' @title RR-BLUP Model 2 #' #' @description -#' Fits an RR-BLUP model for genomic predictions. This implementation is -#' meant for situations where \code{\link{RRBLUP}} is too slow. Note that -#' RRBLUP2 is only faster in certain situations, see details below. Most +#' Fits an RR-BLUP model for genomic predictions. This implementation is +#' meant for situations where \code{\link{RRBLUP}} is too slow. Note that +#' RRBLUP2 is only faster in certain situations, see details below. Most #' users should use \code{\link{RRBLUP}}. -#' +#' #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a -#' function of the traits returning a single value. Unlike \code{\link{RRBLUP}}, +#' function of the traits returning a single value. Unlike \code{\link{RRBLUP}}, #' only univariate models are supported. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. -#' @param maxIter maximum number of iterations. -#' @param Vu marker effect variance. If value is NULL, a +#' @param maxIter maximum number of iterations. +#' @param Vu marker effect variance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param Ve error variance. If value is NULL, a +#' @param Ve error variance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param useEM use EM to solve variance components. If false, +#' @param useEM use EM to solve variance components. If false, #' the initial values are considered true. #' @param tol tolerance for EM algorithm convergence #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits -#' -#' @details -#' The RRBLUP2 function works best when the number of markers is not -#' too large. This is because it solves the RR-BLUP problem by setting -#' up and solving Henderson's mixed model equations. Solving these equations -#' involves a square matrix with dimensions equal to the number of fixed -#' effects plus the number of random effects (markers). Whereas the \code{\link{RRBLUP}} -#' function solves the RR-BLUP problem using the EMMA approach. This approach involves -#' a square matrix with dimensions equal to the number of phenotypic records. This means -#' that the RRBLUP2 function uses less memory than RRBLUP when the number of markers -#' is approximately equal to or smaller than the number of phenotypic records. -#' -#' The RRBLUP2 function is not recommend for cases where the variance components are -#' unknown. This is uses the EM algorithm to solve for unknown variance components, -#' which is generally considerably slower than the EMMA approach of \code{\link{RRBLUP}}. -#' The number of iterations for the EM algorithm is set by maxIter. The default value -#' is typically too small for convergence. When the algorithm fails to converge a -#' warning is displayed, but results are given for the last iteration. These results may -#' be "good enough". However we make no claim to this effect, because we can not generalize +#' +#' @details +#' The RRBLUP2 function works best when the number of markers is not +#' too large. This is because it solves the RR-BLUP problem by setting +#' up and solving Henderson's mixed model equations. Solving these equations +#' involves a square matrix with dimensions equal to the number of fixed +#' effects plus the number of random effects (markers). Whereas the \code{\link{RRBLUP}} +#' function solves the RR-BLUP problem using the EMMA approach. This approach involves +#' a square matrix with dimensions equal to the number of phenotypic records. This means +#' that the RRBLUP2 function uses less memory than RRBLUP when the number of markers +#' is approximately equal to or smaller than the number of phenotypic records. +#' +#' The RRBLUP2 function is not recommend for cases where the variance components are +#' unknown. This is uses the EM algorithm to solve for unknown variance components, +#' which is generally considerably slower than the EMMA approach of \code{\link{RRBLUP}}. +#' The number of iterations for the EM algorithm is set by maxIter. The default value +#' is typically too small for convergence. When the algorithm fails to converge a +#' warning is displayed, but results are given for the last iteration. These results may +#' be "good enough". However we make no claim to this effect, because we can not generalize #' to all possible use cases. -#' -#' The RRBLUP2 function can quickly solve the mixed model equations without estimating variance -#' components. The variance components are set by defining Vu and Ve. Estimation of components -#' is suppressed by setting useEM to false. This may be useful if the model is being retrained -#' multiple times during the simulation. You could run \code{\link{RRBLUP}} function the first -#' time the model is trained, and then use the variance components from this output for all -#' future runs with the RRBLUP2 functions. Again, we can make no claim to the general robustness +#' +#' The RRBLUP2 function can quickly solve the mixed model equations without estimating variance +#' components. The variance components are set by defining Vu and Ve. Estimation of components +#' is suppressed by setting useEM to false. This may be useful if the model is being retrained +#' multiple times during the simulation. You could run \code{\link{RRBLUP}} function the first +#' time the model is trained, and then use the variance components from this output for all +#' future runs with the RRBLUP2 functions. Again, we can make no claim to the general robustness #' of this approach. -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP2(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=10, Vu=NULL, Ve=NULL, - useEM=TRUE, tol=1e-6, simParam=NULL, +RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=10, Vu=NULL, Ve=NULL, + useEM=TRUE, tol=1e-6, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -362,7 +362,7 @@ RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + # Sort out Vu and Ve if(is.function(traits)){ if(is.null(Vu)){ @@ -386,12 +386,12 @@ RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, } } } - + #Fit model ans = callRRBLUP2(y, fixEff, pop@geno, lociPerChr, lociLoc, Vu, Ve, tol, maxIter, useEM, simParam$nThreads) - + bv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -399,7 +399,7 @@ RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$beta), name=paste0("est_BV_",traits)) - + gv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -407,7 +407,7 @@ RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + output = new("RRsol", bv = list(bv), gv = list(gv), @@ -415,68 +415,68 @@ RRBLUP2 = function(pop, traits=1, use="pheno", snpChip=1, male = as.list(NULL), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } #' @title RR-BLUP Model with Dominance #' #' @description -#' Fits an RR-BLUP model for genomic predictions that includes +#' Fits an RR-BLUP model for genomic predictions that includes #' dominance effects. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. -#' @param maxIter maximum number of iterations. Only used +#' @param maxIter maximum number of iterations. Only used #' when number of traits is greater than 1. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitAD(10, meanDD=0.5) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_D(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=40L, +RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=40L, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -486,12 +486,12 @@ RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + #Fit model stopifnot(ncol(y)==1) ans = callRRBLUP_D(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, simParam$nThreads) - + bv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -499,7 +499,7 @@ RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$beta), name=paste0("est_BV_",traits)) - + gv = new("TraitAD", nLoci=nLoci, lociPerChr=lociPerChr, @@ -508,7 +508,7 @@ RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, domEff=c(ans$d), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + output = new("RRsol", bv = list(bv), gv = list(gv), @@ -516,7 +516,7 @@ RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, male = as.list(NULL), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } @@ -524,74 +524,74 @@ RRBLUP_D = function(pop, traits=1, use="pheno", snpChip=1, #' @title RR-BLUP with Dominance Model 2 #' #' @description -#' Fits an RR-BLUP model for genomic predictions that includes -#' dominance effects. This implementation is meant for situations where -#' \code{\link{RRBLUP_D}} is too slow. Note that RRBLUP_D2 -#' is only faster in certain situations. Most users should use +#' Fits an RR-BLUP model for genomic predictions that includes +#' dominance effects. This implementation is meant for situations where +#' \code{\link{RRBLUP_D}} is too slow. Note that RRBLUP_D2 +#' is only faster in certain situations. Most users should use #' \code{\link{RRBLUP_D}}. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. -#' @param maxIter maximum number of iterations. Only used +#' @param maxIter maximum number of iterations. Only used #' when number of traits is greater than 1. -#' @param Va marker effect variance for additive effects. If value is NULL, +#' @param Va marker effect variance for additive effects. If value is NULL, #' a reasonable starting point is chosen automatically. -#' @param Vd marker effect variance for dominance effects. If value is NULL, +#' @param Vd marker effect variance for dominance effects. If value is NULL, #' a reasonable starting point is chosen automatically. -#' @param Ve error variance. If value is NULL, a +#' @param Ve error variance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param useEM use EM to solve variance components. If false, +#' @param useEM use EM to solve variance components. If false, #' the initial values are considered true. #' @param tol tolerance for EM algorithm convergence #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitAD(10, meanDD=0.5) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_D2(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=10, Va=NULL, Vd=NULL, - Ve=NULL, useEM=TRUE, tol=1e-6, +RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=10, Va=NULL, Vd=NULL, + Ve=NULL, useEM=TRUE, tol=1e-6, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -601,7 +601,7 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + # Sort out Va, Vd and Ve if(is.function(traits)){ if(is.null(Va)){ @@ -634,13 +634,13 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, } } } - + #Fit model stopifnot(ncol(y)==1) ans = callRRBLUP_D2(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, Va, Vd, Ve, tol, useEM, simParam$nThreads) - + bv = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -648,7 +648,7 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha), intercept=c(ans$beta), name=paste0("est_BV_",traits)) - + gv = new("TraitAD", nLoci=nLoci, lociPerChr=lociPerChr, @@ -657,7 +657,7 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, domEff=c(ans$d), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + output = new("RRsol", bv = list(bv), gv = list(gv), @@ -665,7 +665,7 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, male = as.list(NULL), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } @@ -674,60 +674,60 @@ RRBLUP_D2 = function(pop, traits=1, use="pheno", snpChip=1, #' @description #' Fits an RR-BLUP model that estimates seperate marker effects for #' females and males. Useful for predicting GCA of parents -#' in single cross hybrids. Can also predict performance of specific +#' in single cross hybrids. Can also predict performance of specific #' single cross hybrids. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. #' @param maxIter maximum number of iterations for convergence. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_GCA(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=40L, +RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=40L, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -737,13 +737,13 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + #Fit model stopifnot(ncol(y)==1) ans = callRRBLUP_GCA(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, simParam$nThreads) - + gv = new("TraitA2", nLoci=nLoci, lociPerChr=lociPerChr, @@ -752,7 +752,7 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, addEffMale=c(ans$alpha2), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + female = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -760,7 +760,7 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha1), intercept=c(ans$beta1), name=paste0("est_female_",traits)) - + male = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -768,7 +768,7 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha2), intercept=c(ans$beta2), name=paste0("est_male_",traits)) - + output = new("RRsol", gv = list(gv), bv = as.list(NULL), @@ -776,7 +776,7 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, male = list(male), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } @@ -784,72 +784,72 @@ RRBLUP_GCA = function(pop, traits=1, use="pheno", snpChip=1, #' #' @description #' Fits an RR-BLUP model that estimates seperate marker effects for -#' females and males. This implementation is meant for situations where -#' \code{\link{RRBLUP_GCA}} is too slow. Note that RRBLUP_GCA2 -#' is only faster in certain situations. Most users should use +#' females and males. This implementation is meant for situations where +#' \code{\link{RRBLUP_GCA}} is too slow. Note that RRBLUP_GCA2 +#' is only faster in certain situations. Most users should use #' \code{\link{RRBLUP_GCA}}. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. #' @param maxIter maximum number of iterations for convergence. -#' @param VuF marker effect variance for females. If value is NULL, a +#' @param VuF marker effect variance for females. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param VuM marker effect variance for males. If value is NULL, a +#' @param VuM marker effect variance for males. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param Ve error variance. If value is NULL, a +#' @param Ve error variance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param useEM use EM to solve variance components. If false, +#' @param useEM use EM to solve variance components. If false, #' the initial values are considered true. #' @param tol tolerance for EM algorithm convergence #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_GCA2(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=10, VuF=NULL, VuM=NULL, - Ve=NULL, useEM=TRUE, tol=1e-6, +RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=10, VuF=NULL, VuM=NULL, + Ve=NULL, useEM=TRUE, tol=1e-6, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -859,7 +859,7 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + # Sort out VuF, VuM and Ve if(is.function(traits)){ if(is.null(VuF)){ @@ -892,15 +892,15 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, } } } - + #Fit model stopifnot(ncol(y)==1) - + ans = callRRBLUP_GCA2(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, VuF, VuM, Ve, tol, useEM, simParam$nThreads) - + gv = new("TraitA2", nLoci=nLoci, lociPerChr=lociPerChr, @@ -909,7 +909,7 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, addEffMale=c(ans$alpha2), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + female = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -917,7 +917,7 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha1), intercept=c(ans$beta1), name=paste0("est_female_",traits)) - + male = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -925,7 +925,7 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha2), intercept=c(ans$beta2), name=paste0("est_male_",traits)) - + output = new("RRsol", gv = list(gv), bv = as.list(NULL), @@ -933,68 +933,68 @@ RRBLUP_GCA2 = function(pop, traits=1, use="pheno", snpChip=1, male = list(male), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } #' @title RR-BLUP SCA Model #' #' @description -#' An extention of \code{\link{RRBLUP_GCA}} that adds dominance effects. -#' Note that we have not seen any consistent benefit of this model over +#' An extention of \code{\link{RRBLUP_GCA}} that adds dominance effects. +#' Note that we have not seen any consistent benefit of this model over #' \code{\link{RRBLUP_GCA}}. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. #' @param maxIter maximum number of iterations for convergence. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_SCA(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=40L, +RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=40L, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -1004,13 +1004,13 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + #Fit model stopifnot(ncol(y)==1) ans = callRRBLUP_SCA(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, simParam$nThreads) - + gv = new("TraitA2D", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1020,7 +1020,7 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, domEff=c(ans$d), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + female = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1028,7 +1028,7 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha1), intercept=c(ans$beta1), name=paste0("est_female_",traits)) - + male = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1036,7 +1036,7 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha2), intercept=c(ans$beta2), name=paste0("est_male_",traits)) - + output = new("RRsol", gv = list(gv), bv = as.list(NULL), @@ -1044,7 +1044,7 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, male = list(male), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } @@ -1052,74 +1052,74 @@ RRBLUP_SCA = function(pop, traits=1, use="pheno", snpChip=1, #' #' @description #' Fits an RR-BLUP model that estimates seperate additive effects for -#' females and males and a dominance effect. This implementation is meant -#' for situations where \code{\link{RRBLUP_SCA}} is too slow. Note that -#' RRBLUP_SCA2 is only faster in certain situations. Most users should use +#' females and males and a dominance effect. This implementation is meant +#' for situations where \code{\link{RRBLUP_SCA}} is too slow. Note that +#' RRBLUP_SCA2 is only faster in certain situations. Most users should use #' \code{\link{RRBLUP_SCA}}. #' #' @param pop a \code{\link{Pop-class}} to serve as the training population #' @param traits an integer indicating the trait to model, a trait name, or a #' function of the traits returning a single value. -#' @param use train model using phenotypes "pheno", genetic values "gv", +#' @param use train model using phenotypes "pheno", genetic values "gv", #' estimated breeding values "ebv", breeding values "bv", or randomly "rand" -#' @param snpChip an integer indicating which SNP chip genotype +#' @param snpChip an integer indicating which SNP chip genotype #' to use -#' @param useQtl should QTL genotypes be used instead of a SNP chip. -#' If TRUE, snpChip specifies which trait's QTL to use, and thus these +#' @param useQtl should QTL genotypes be used instead of a SNP chip. +#' If TRUE, snpChip specifies which trait's QTL to use, and thus these #' QTL may not match the QTL underlying the phenotype supplied in traits. #' @param maxIter maximum number of iterations for convergence. -#' @param VuF marker effect variance for females. If value is NULL, a +#' @param VuF marker effect variance for females. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param VuM marker effect variance for males. If value is NULL, a +#' @param VuM marker effect variance for males. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param VuD marker effect variance for dominance. If value is NULL, a +#' @param VuD marker effect variance for dominance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param Ve error variance. If value is NULL, a +#' @param Ve error variance. If value is NULL, a #' reasonable starting point is chosen automatically. -#' @param useEM use EM to solve variance components. If false, +#' @param useEM use EM to solve variance components. If false, #' the initial values are considered true. #' @param tol tolerance for EM algorithm convergence #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' traits #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP_SCA2(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, - useQtl=FALSE, maxIter=10, VuF=NULL, VuM=NULL, - VuD=NULL, Ve=NULL, useEM=TRUE, tol=1e-6, +RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, + useQtl=FALSE, maxIter=10, VuF=NULL, VuM=NULL, + VuD=NULL, Ve=NULL, useEM=TRUE, tol=1e-6, simParam=NULL, ...){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + y = getResponse(pop=pop,trait=traits,use=use, simParam=simParam,...) - + traits = convertTraitsToNames(traits, simParam) - + fixEff = as.integer(factor(pop@fixEff)) - + if(useQtl){ nLoci = simParam$traits[[snpChip]]@nLoci lociPerChr = simParam$traits[[snpChip]]@lociPerChr @@ -1129,7 +1129,7 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, lociPerChr = simParam$snpChips[[snpChip]]@lociPerChr lociLoc = simParam$snpChips[[snpChip]]@lociLoc } - + # Sort out VuF, VuM, VuD and Ve if(is.function(traits)){ if(is.null(VuF)){ @@ -1171,14 +1171,14 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, } } } - + #Fit model stopifnot(ncol(y)==1) ans = callRRBLUP_SCA2(y, fixEff, pop@geno, lociPerChr, lociLoc, maxIter, VuF, VuM, VuD, Ve, tol, useEM, simParam$nThreads) - + gv = new("TraitA2D", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1188,7 +1188,7 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, domEff=c(ans$d), intercept=c(ans$mu), name=paste0("est_GV_",traits)) - + female = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1196,7 +1196,7 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha1), intercept=c(ans$beta1), name=paste0("est_female_",traits)) - + male = new("TraitA", nLoci=nLoci, lociPerChr=lociPerChr, @@ -1204,7 +1204,7 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, addEff=c(ans$alpha2), intercept=c(ans$beta2), name=paste0("est_male_",traits)) - + output = new("RRsol", gv = list(gv), bv = as.list(NULL), @@ -1212,119 +1212,118 @@ RRBLUP_SCA2 = function(pop, traits=1, use="pheno", snpChip=1, male = list(male), Vu = as.matrix(ans$Vu), Ve = as.matrix(ans$Ve)) - + return(output) } -#' @title Set EBV +#' @title Set estimated breeding values (EBV) #' #' @description -#' Adds genomic estimated values to a populations's EBV -#' slot using output from a genomic selection functions. -#' The genomic estimated values can be either estimated -#' breeding values, estimated genetic values, or +#' Adds genomic estimated values to a populations's EBV +#' slot using output from a genomic selection functions. +#' The genomic estimated values can be either estimated +#' breeding values, estimated genetic values, or #' estimated general combining values. #' #' @param pop an object of \code{\link{Pop-class}} #' @param solution an object of \code{\link{RRsol-class}} -#' @param value the genomic value to be estimated. Can be +#' @param value the genomic value to be estimated. Can be #' either "gv", "bv", "female", or "male". -#' @param targetPop an optional target population that can -#' be used when value is "bv", "female", or "male". When -#' supplied, the allele frequency in the targetPop is used +#' @param targetPop an optional target population that can +#' be used when value is "bv", "female", or "male". When +#' supplied, the allele frequency in the targetPop is used #' to set these values. -#' @param append should estimated values be appended to -#' existing data in the EBV slot. If TRUE, a new column is -#' added. If FALSE, existing data is replaced with the +#' @param append should estimated values be appended to +#' existing data in the EBV slot. If TRUE, a new column is +#' added. If FALSE, existing data is replaced with the #' new estimates. #' @param simParam an object of \code{\link{SimParam}} #' -#' +#' #' @return Returns an object of \code{\link{Pop-class}} #' -#' @examples +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=20) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) #' SP$setVarE(h2=0.5) #' SP$addSnpChip(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Run GS model and set EBV #' ans = RRBLUP(pop, simParam=SP) #' pop = setEBV(pop, ans, simParam=SP) -#' +#' #' #Evaluate accuracy #' cor(gv(pop), ebv(pop)) -#' +#' #' @export -setEBV = function(pop, solution, value="gv", targetPop=NULL, +setEBV = function(pop, solution, value="gv", targetPop=NULL, append=FALSE, simParam=NULL){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - + nTraits = length(solution@gv) - + ebv = matrix(NA_real_, nrow=pop@nInd, ncol=nTraits) - + # Placeholder names colnames(ebv) = as.character(1:nTraits) - + value = tolower(value) - + if(value=="gv"){ - - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ tmp = getGv(solution@gv[[i]],pop,simParam$nThreads) ebv[,i] = tmp[[1]] colnames(ebv)[i] = solution@gv[[i]]@name } - + }else if(value=="bv"){ - + if(is.null(targetPop)){ - + if(length(solution@bv)==0){ stop("This genomic selection model does not produce breeding value estimates.") } - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ tmp = getGv(solution@bv[[i]],pop,simParam$nThreads) ebv[,i] = tmp[[1]] colnames(ebv)[i] = solution@bv[[i]]@name } - + }else{ - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = solution@gv[[i]] if(.hasSlot(trait,"addEffMale")){ stop("This genomic selection model does not produce breeding value estimates. Try value='male' or value='female' instead.") } - - p = calcGenoFreq(targetPop@geno, - trait@lociPerChr, - trait@lociLoc, + + p = calcGenoFreq(targetPop@geno, + trait@lociPerChr, + trait@lociLoc, simParam$nThreads) p = c(p) q = 1-p - + a = trait@addEff if(.hasSlot(trait,"domEff")){ d = trait@domEff }else{ d = rep(0, length(a)) } - + alpha = a+d*(q-p) intercept = -sum((p-q)*alpha) trait = new("TraitA", @@ -1333,51 +1332,51 @@ setEBV = function(pop, solution, value="gv", targetPop=NULL, lociLoc=trait@lociLoc, addEff=alpha, intercept=intercept) - + tmp = getGv(trait, pop, simParam$nThreads) ebv[,i] = tmp[[1]] - + # changing original name from "est_GV_..." to "est_BV_..." tmp = solution@gv[[i]]@name tmp = strsplit(tmp, "_")[[1]] tmp[2] = "BV" colnames(ebv)[i] = paste(tmp,collapse="_") } - + } - + }else if(value=="female"){ - + if(is.null(targetPop)){ - + if(length(solution@female)==0){ stop("This genomic selection model does not produce GCA estimates for females.") } - - for(i in 1:nTraits){ + + for(i in seq_len(nTraits)){ tmp = getGv(solution@female[[i]],pop,simParam$nThreads) ebv[,i] = tmp[[1]] colnames(ebv)[i] = solution@female[[i]]@name } - + }else{ - - for(i in 1:nTraits){ + + for(i in seq_len(nTraits)){ trait = solution@gv[[i]] - p = calcGenoFreq(targetPop@geno, - trait@lociPerChr, - trait@lociLoc, + p = calcGenoFreq(targetPop@geno, + trait@lociPerChr, + trait@lociLoc, simParam$nThreads) p = c(p) q = 1-p - + a = trait@addEff if(.hasSlot(trait,"domEff")){ d = trait@domEff }else{ d = rep(0, length(a)) } - + alpha = (a+d*(q-p))/2 intercept = -sum((p-q)*alpha) trait = new("TraitA", @@ -1386,44 +1385,44 @@ setEBV = function(pop, solution, value="gv", targetPop=NULL, lociLoc=trait@lociLoc, addEff=alpha, intercept=intercept) - + tmp = getGv(trait, pop, simParam$nThreads) ebv[,i] = tmp[[1]] - + # changing original name from "est_GV_..." to "est_female_..." tmp = solution@gv[[i]]@name tmp = strsplit(tmp, "_")[[1]] tmp[2] = "female" colnames(ebv)[i] = paste(tmp,collapse="_") } - + } - + }else if(value=="male"){ - + if(is.null(targetPop)){ - + if(length(solution@male)==0){ stop("This genomic selection model does not produce GCA estimates for males.") } - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ tmp = getGv(solution@male[[i]],pop,simParam$nThreads) ebv[,i] = tmp[[1]] colnames(ebv)[i] = solution@male[[i]]@name } - + }else{ - - for(i in 1:nTraits){ + + for(i in seq_len(nTraits)){ trait = solution@gv[[i]] - p = calcGenoFreq(targetPop@geno, - trait@lociPerChr, - trait@lociLoc, + p = calcGenoFreq(targetPop@geno, + trait@lociPerChr, + trait@lociLoc, simParam$nThreads) p = c(p) q = 1-p - + if(.hasSlot(trait,"addEffMale")){ a = trait@addEffMale }else{ @@ -1434,7 +1433,7 @@ setEBV = function(pop, solution, value="gv", targetPop=NULL, }else{ d = rep(0, length(a)) } - + alpha = (a+d*(q-p))/2 intercept = -sum((p-q)*alpha) trait = new("TraitA", @@ -1443,50 +1442,49 @@ setEBV = function(pop, solution, value="gv", targetPop=NULL, lociLoc=trait@lociLoc, addEff=alpha, intercept=intercept) - + tmp = getGv(trait, pop, simParam$nThreads) ebv[,i] = tmp[[1]] - + # changing original name from "est_BV_..." to "est_male_..." tmp = solution@gv[[i]]@name tmp = strsplit(tmp, "_")[[1]] tmp[2] = "male" colnames(ebv)[i] = paste(tmp,collapse="_") } - + } - + }else{ stop(paste0("value=",value," is not a valid option")) } - + if(append){ pop@ebv = cbind(pop@ebv,ebv) }else{ pop@ebv = ebv } - + return(pop) } - #' @title RRBLUP Memory Usage #' #' @description #' Estimates the amount of RAM needed to run the \code{\link{RRBLUP}} -#' and its related functions for a given training population size. +#' and its related functions for a given training population size. #' Note that this function may underestimate total usage. #' #' @param nInd the number of individuals in the training population #' @param nMarker the number of markers per individual -#' @param model either "REG", "GCA", or "SCA" for \code{\link{RRBLUP}} +#' @param model either "REG", "GCA", or "SCA" for \code{\link{RRBLUP}} #' \code{\link{RRBLUP_GCA}} and \code{\link{RRBLUP_SCA}} respectively. #' #' @return Returns an estimate for the required gigabytes of RAM #' -#' @examples +#' @examples #' RRBLUPMemUse(nInd=1000, nMarker=5000) -#' +#' #' @export RRBLUPMemUse = function(nInd,nMarker,model="REG"){ y = nInd @@ -1518,10 +1516,10 @@ RRBLUPMemUse = function(nInd,nMarker,model="REG"){ objects = objects[objects!="model"] bytes = sapply(objects,function(x) get(x)) bytes = 8*sum(bytes) - + # Using base 2 calculation #return(bytes/((2^10)^3)) #GB - + # Using base 10 calculation, more conservative # Incorrect, but accounts for sources of data usage return(bytes/10^9) #GB diff --git a/R/crossing.R b/R/crossing.R index 81e2ad9e..bf6642ff 100644 --- a/R/crossing.R +++ b/R/crossing.R @@ -159,7 +159,7 @@ randCross = function(pop,nCrosses,nProgeny=1, tmp = male[sample.int(nMale, nMale)] n = nCrosses%/%nMale + 1 male = NULL - for(i in 1:n){ + for(i in seq_len(n)){ take = nMale - (i:(nMale+i-1))%%nMale male = c(male, tmp[take]) } @@ -435,7 +435,7 @@ randCross2 = function(females,males,nCrosses,nProgeny=1, tmp = male[sample.int(nMale, nMale)] n = nCrosses%/%nMale + 1 male = NULL - for(i in 1:n){ + for(i in seq_len(n)){ take = nMale - (i:(nMale+i-1))%%nMale male = c(male, tmp[take]) } @@ -666,7 +666,7 @@ sortPed = function(id, mother, father, maxCycle=100){ motherID=as.character(mother), fatherID=as.character(father)) unsorted = rep(TRUE, nInd) - for(gen in 1:maxCycle){ + for(gen in seq_len(maxCycle)){ for(i in which(unsorted)){ if(is.na(output$mother[i])&is.na(output$father[i])){ # Is a founder @@ -751,7 +751,6 @@ sortPed = function(id, mother, father, maxCycle=100){ #' father = c(0,0,2,3:9) #' pop2 = pedigreeCross(pop, id, mother, father, simParam=SP) #' -#' #' @export pedigreeCross = function(founderPop, id, mother, father, matchID=FALSE, maxCycle=100, DH=NULL, nSelf=NULL, useFemale=TRUE, @@ -844,7 +843,7 @@ pedigreeCross = function(founderPop, id, mother, father, matchID=FALSE, # Create individuals crossPlan = matrix(c(1,1),ncol=2) - for(gen in 1:max(ped$gen)){ + for(gen in seq_len(max(ped$gen))){ for(i in which(ped$gen==gen)){ if(isFounder[i]){ # Copy over founder individual @@ -870,13 +869,11 @@ pedigreeCross = function(founderPop, id, mother, father, matchID=FALSE, simParam=simParam) } } - + # Self? - if(nSelf[i]>0){ - for(j in 1:nSelf[i]){ - output[[i]] = self(output[[i]], - simParam=simParam) - } + for(j in seq_len(nSelf[i])){ + output[[i]] = self(output[[i]], + simParam=simParam) } # Make the individual a DH? diff --git a/R/founderPop.R b/R/founderPop.R index ca887676..5cfacd43 100644 --- a/R/founderPop.R +++ b/R/founderPop.R @@ -1,44 +1,44 @@ #' @title New MapPop #' -#' @description -#' Creates a new \code{\link{MapPop-class}} from user supplied +#' @description +#' Creates a new \code{\link{MapPop-class}} from user supplied #' genetic maps and haplotypes. -#' +#' #' @param genMap a list of genetic maps -#' @param haplotypes a list of matrices or data.frames that +#' @param haplotypes a list of matrices or data.frames that #' can be coerced to matrices. See details. #' @param inbred are individuals fully inbred #' @param ploidy ploidy level of the organism -#' +#' #' @details -#' Each item of genMap must be a vector of ordered genetic lengths in -#' Morgans. The first value must be zero. The length of the vector +#' Each item of genMap must be a vector of ordered genetic lengths in +#' Morgans. The first value must be zero. The length of the vector #' determines the number of segregating sites on the chromosome. -#' -#' Each item of haplotypes must be coercible to a matrix. The columns -#' of this matrix correspond to segregating sites. The number of rows -#' must match the number of individuals times the ploidy if using -#' inbred=FALSE. If using inbred=TRUE, the number of rows must equal -#' the number of individuals. The haplotypes can be stored as numeric, +#' +#' Each item of haplotypes must be coercible to a matrix. The columns +#' of this matrix correspond to segregating sites. The number of rows +#' must match the number of individuals times the ploidy if using +#' inbred=FALSE. If using inbred=TRUE, the number of rows must equal +#' the number of individuals. The haplotypes can be stored as numeric, #' integer or raw. The underlying C++ function will use raw. -#' +#' #' @return an object of \code{\link{MapPop-class}} -#' -#' @examples +#' +#' @examples #' # Create genetic map for two chromosomes, each 1 Morgan long #' # Each chromosome contains 11 equally spaced segregating sites #' genMap = list(seq(0,1,length.out=11), #' seq(0,1,length.out=11)) -#' +#' #' # Create haplotypes for 10 outbred individuals #' chr1 = sample(x=0:1,size=20*11,replace=TRUE) #' chr1 = matrix(chr1,nrow=20,ncol=11) #' chr2 = sample(x=0:1,size=20*11,replace=TRUE) #' chr2 = matrix(chr2,nrow=20,ncol=11) #' haplotypes = list(chr1,chr2) -#' +#' #' founderPop = newMapPop(genMap=genMap, haplotypes=haplotypes) -#' +#' #' @export newMapPop = function(genMap, haplotypes, inbred=FALSE, ploidy=2L){ @@ -68,7 +68,7 @@ newMapPop = function(genMap, haplotypes, inbred=FALSE, stop("Number of segregating sites in haplotypes and genMap don't match") } output = vector("list",length(genMap)) - for(chr in 1:length(genMap)){ + for(chr in seq_len(length(genMap))){ geno = packHaplo(as.matrix(haplotypes[[chr]]), ploidy=ploidy,inbred=inbred) output[[chr]] = new("MapPop", @@ -93,49 +93,49 @@ newMapPop = function(genMap, haplotypes, inbred=FALSE, #' @title Create founder haplotypes using MaCS #' -#' @description Uses the MaCS software to produce founder haplotypes +#' @description Uses the MaCS software to produce founder haplotypes #' \insertCite{MaCS}{AlphaSimR}. -#' +#' #' @param nInd number of individuals to simulate #' @param nChr number of chromosomes to simulate -#' @param segSites number of segregating sites to keep per chromosome. A +#' @param segSites number of segregating sites to keep per chromosome. A #' value of NULL results in all sites being retained. #' @param inbred should founder individuals be inbred #' @param species species history to simulate. See details. #' @param split an optional historic population split in terms of generations ago. #' @param ploidy ploidy level of organism #' @param manualCommand user provided MaCS options. For advanced users only. -#' @param manualGenLen user provided genetic length. This must be supplied if using -#' manualCommand. If not using manualCommand, this value will replace the predefined -#' genetic length for the species. However, this the genetic length is only used by -#' AlphaSimR and is not passed to MaCS, so MaCS still uses the predefined genetic length. +#' @param manualGenLen user provided genetic length. This must be supplied if using +#' manualCommand. If not using manualCommand, this value will replace the predefined +#' genetic length for the species. However, this the genetic length is only used by +#' AlphaSimR and is not passed to MaCS, so MaCS still uses the predefined genetic length. #' For advanced users only. -#' @param nThreads if OpenMP is available, this will allow for simulating chromosomes in parallel. +#' @param nThreads if OpenMP is available, this will allow for simulating chromosomes in parallel. #' If the value is NULL, the number of threads is automatically detected. -#' +#' #' @details #' There are currently three species histories available: GENERIC, CATTLE, WHEAT, and MAIZE. -#' -#' The GENERIC history is meant to be a reasonable all-purpose choice. It runs quickly and -#' models a population with an effective populations size that has gone through several historic -#' bottlenecks. This species history is used as the default arguments in the \code{\link{runMacs2}} +#' +#' The GENERIC history is meant to be a reasonable all-purpose choice. It runs quickly and +#' models a population with an effective populations size that has gone through several historic +#' bottlenecks. This species history is used as the default arguments in the \code{\link{runMacs2}} #' function, so the user should examine this function for the details of how the species is modeled. -#' +#' #' The CATTLE history is based off of real genome sequence data \insertCite{cattle}{AlphaSimR}. -#' -#' The WHEAT \insertCite{gaynor_2017}{AlphaSimR} and MAIZE \insertCite{hickey_2014}{AlphaSimR} -#' histories have been included due to their use in previous simulations. However, it should -#' be noted that neither faithfully simulates its respective species. This is apparent by -#' the low number of segregating sites simulated by each history relative to their real-world -#' analogs. Adjusting these histories to better represent their real-world analogs would result +#' +#' The WHEAT \insertCite{gaynor_2017}{AlphaSimR} and MAIZE \insertCite{hickey_2014}{AlphaSimR} +#' histories have been included due to their use in previous simulations. However, it should +#' be noted that neither faithfully simulates its respective species. This is apparent by +#' the low number of segregating sites simulated by each history relative to their real-world +#' analogs. Adjusting these histories to better represent their real-world analogs would result #' in a drastic increase to runtime. #' #' @return an object of \code{\link{MapPop-class}} -#' -#' @references +#' +#' @references #' \insertAllCited{} -#' -#' @examples +#' +#' @examples #' # Creates a populations of 10 outbred individuals #' # Their genome consists of 1 chromosome and 100 segregating sites #' \dontrun{ @@ -154,61 +154,74 @@ runMacs = function(nInd,nChr=1, segSites=NULL, inbred=FALSE, species="GENERIC", if(nChr0 if(any(nLoci[isLimited] != segSites[isLimited])){ stop("MaCS did not return enough segSites, use segSites=NULL to return all sites generated by MaCS") } - + genMap = vector("list",nChr) - for(i in 1:nChr){ + for(i in seq_len(nChr)){ genMap[[i]] = genLen[i]*c(macsOut$genMap[[i]]-macsOut$genMap[[i]][1]) names(genMap[[i]]) = paste(i,1:length(genMap[[i]]),sep="_") } names(genMap) = as.character(1:nChr) - + output = new("MapPop", nInd=nInd, nChr=nChr, @@ -261,15 +274,15 @@ runMacs = function(nInd,nChr=1, segSites=NULL, inbred=FALSE, species="GENERIC", #' @title Alternative wrapper for MaCS #' -#' @description -#' A wrapper function for \code{\link{runMacs}}. This wrapper is designed +#' @description +#' A wrapper function for \code{\link{runMacs}}. This wrapper is designed #' to provide a more intuitive interface for writing custom commands -#' in MaCS \insertCite{MaCS}{AlphaSimR}. It effectively automates the creation -#' of an appropriate line for the manualCommand argument in \code{\link{runMacs}} -#' using user supplied variables, but only allows for a subset of the functionality -#' offered by this argument. The default arguments of this function were chosen to match +#' in MaCS \insertCite{MaCS}{AlphaSimR}. It effectively automates the creation +#' of an appropriate line for the manualCommand argument in \code{\link{runMacs}} +#' using user supplied variables, but only allows for a subset of the functionality +#' offered by this argument. The default arguments of this function were chosen to match #' species="GENERIC" in \code{\link{runMacs}}. -#' +#' #' @param nInd number of individuals to simulate #' @param nChr number of chromosomes to simulate #' @param segSites number of segregating sites to keep per chromosome @@ -277,32 +290,45 @@ runMacs = function(nInd,nChr=1, segSites=NULL, inbred=FALSE, species="GENERIC", #' @param bp base pair length of chromosome #' @param genLen genetic length of chromosome in Morgans #' @param mutRate per base pair mutation rate -#' @param histNe effective population size in previous +#' @param histNe effective population size in previous #' generations -#' @param histGen number of generations ago for effective +#' @param histGen number of generations ago for effective #' population sizes given in histNe #' @param inbred should founder individuals be inbred #' @param split an optional historic population split in terms of generations ago #' @param ploidy ploidy level of organism -#' @param returnCommand should the command passed to manualCommand in -#' \code{\link{runMacs}} be returned. If TRUE, MaCS will not be called and +#' @param returnCommand should the command passed to manualCommand in +#' \code{\link{runMacs}} be returned. If TRUE, MaCS will not be called and #' the command is returned instead. -#' @param nThreads if OpenMP is available, this will allow for simulating chromosomes in parallel. +#' @param nThreads if OpenMP is available, this will allow for simulating chromosomes in parallel. #' If the value is NULL, the number of threads is automatically detected. #' -#' @return an object of \code{\link{MapPop-class}} or if -#' returnCommand is true a string giving the MaCS command passed to +#' @return an object of \code{\link{MapPop-class}} or if +#' returnCommand is true a string giving the MaCS command passed to #' the manualCommand argument of \code{\link{runMacs}}. -#' -#' @references +#' +#' @references #' \insertAllCited{} -#' -#' @examples +#' +#' @examples #' # Creates a populations of 10 outbred individuals #' # Their genome consists of 1 chromosome and 100 segregating sites #' # The command is equivalent to using species="GENERIC" in runMacs #' \dontrun{ #' founderPop = runMacs2(nInd=10,nChr=1,segSites=100) +#' +#' # runMacs() Implementation of the cattle demography following +#' # Macleod et al. (2013) https://doi.org/10.1093/molbev/mst125 +#' cattleChrSum = 2.8e9 # https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_002263795.3/ +#' (cattleChrBp = cattleChrSum / 30) +#' recRate = 9.26e-09 +#' (cattleGenLen = recRate * cattleChrBp) +#' mutRate = 1.20e-08 +#' runMacs2(nInd = 10, nChr = 1, Ne = 90, bp = cattleChrBp, +#' genLen = cattleGenLen, mutRate = 1.20e-08, +#' histNe = c(120, 250, 350, 1000, 1500, 2000, 2500, 3500, 7000, 10000, 17000, 62000), +#' histGen = c( 3, 6, 12, 18, 24, 154, 454, 654, 1754, 2354, 3354, 33154), +#' returnCommand = TRUE) #' } #' @export runMacs2 = function(nInd,nChr=1,segSites=NULL,Ne=100, @@ -323,7 +349,7 @@ runMacs2 = function(nInd,nChr=1,segSites=NULL,Ne=100, if(length(histNe)>0){ histNe = histNe/Ne histGen = histGen/(4*Ne) - for(i in 1:length(histNe)){ + for(i in seq_len(length(histNe))){ speciesHist = paste(speciesHist,"-eN", histGen[i],histNe[i]) } @@ -348,24 +374,24 @@ runMacs2 = function(nInd,nChr=1,segSites=NULL,Ne=100, #' @title Sample haplotypes from a MapPop #' -#' @description -#' Creates a new \code{\link{MapPop-class}} from an existing +#' @description +#' Creates a new \code{\link{MapPop-class}} from an existing #' \code{\link{MapPop-class}} by randomly sampling haplotypes. -#' -#' @param mapPop the \code{\link{MapPop-class}} used to +#' +#' @param mapPop the \code{\link{MapPop-class}} used to #' sample haplotypes #' @param nInd the number of individuals to create #' @param inbred should new individuals be fully inbred -#' @param ploidy new ploidy level for organism. If NULL, +#' @param ploidy new ploidy level for organism. If NULL, #' the ploidy level of the mapPop is used. #' @param replace should haplotypes be sampled with replacement -#' +#' #' @return an object of \code{\link{MapPop-class}} -#' -#' @examples +#' +#' @examples #' founderPop = quickHaplo(nInd=2,nChr=1,segSites=11,inbred=TRUE) #' founderPop = sampleHaplo(mapPop=founderPop,nInd=20) -#' +#' #' @export sampleHaplo = function(mapPop,nInd,inbred=FALSE,ploidy=NULL,replace=TRUE){ if(is.null(ploidy)) ploidy = mapPop@ploidy @@ -378,24 +404,24 @@ sampleHaplo = function(mapPop,nInd,inbred=FALSE,ploidy=NULL,replace=TRUE){ nBin = mapPop@nLoci%/%8L + (mapPop@nLoci%%8L > 0L) if(!replace) stopifnot(nHaplo>=nSamp) output = vector("list",mapPop@nChr) - for(chr in 1:mapPop@nChr){ + for(chr in seq_len(mapPop@nChr)){ haplo = sample.int(nHaplo,nSamp,replace=replace) geno = array(data=as.raw(0), dim=c(nBin[chr], ploidy,nInd)) - outHap = 1L + outHap = 1L outInd = 1L - for(i in 1:length(haplo)){ + for(i in seq_len(length(haplo))){ inHap = (haplo[i]-1L)%%mapPop@ploidy + 1L inInd = (haplo[i]-1L)%/%mapPop@ploidy + 1L if(inbred){ - for(outHap in 1:ploidy){ - geno[,outHap,outInd] = + for(outHap in seq_len(ploidy)){ + geno[,outHap,outInd] = mapPop@geno[[chr]][,inHap,inInd] } outInd = outInd+1L }else{ - geno[,outHap,outInd] = + geno[,outHap,outInd] = mapPop@geno[[chr]][,inHap,inInd] outHap = outHap%%ploidy+1L if(outHap==1L){ @@ -419,10 +445,10 @@ sampleHaplo = function(mapPop,nInd,inbred=FALSE,ploidy=NULL,replace=TRUE){ #' @title Quick founder haplotype simulation #' -#' @description Rapidly simulates founder haplotypes by randomly -#' sampling 0s and 1s. This is equivalent to having all loci with +#' @description Rapidly simulates founder haplotypes by randomly +#' sampling 0s and 1s. This is equivalent to having all loci with #' allele frequency 0.5 and being in linkage equilibrium. -#' +#' #' @param nInd number of individuals to simulate #' @param nChr number of chromosomes to simulate #' @param segSites number of segregating sites per chromosome @@ -431,12 +457,12 @@ sampleHaplo = function(mapPop,nInd,inbred=FALSE,ploidy=NULL,replace=TRUE){ #' @param inbred should founder individuals be inbred #' #' @return an object of \code{\link{MapPop-class}} -#' -#' @examples +#' +#' @examples #' # Creates a populations of 10 outbred individuals #' # Their genome consists of 1 chromosome and 100 segregating sites #' founderPop = quickHaplo(nInd=10,nChr=1,segSites=100) -#' +#' #' @export quickHaplo = function(nInd,nChr,segSites,genLen=1,ploidy=2L,inbred=FALSE){ ploidy = as.integer(ploidy) @@ -447,17 +473,17 @@ quickHaplo = function(nInd,nChr,segSites,genLen=1,ploidy=2L,inbred=FALSE){ if(length(genLen)==1) genLen = rep(genLen,nChr) nBins = segSites%/%8L + (segSites%%8L > 0L) centromere = genLen/2 - + genMap = vector("list",nChr) geno = vector("list",nChr) - for(i in 1:nChr){ + for(i in seq_len(nChr)){ genMap[[i]] = seq(0,genLen[i],length.out=segSites[i]) names(genMap[[i]]) = paste(i, 1:segSites[i], sep="_") geno[[i]] = array(sample(as.raw(0:255), nInd*ploidy*nBins[i], replace=TRUE), dim = c(nBins[i],ploidy,nInd)) - + if(inbred){ if(ploidy>1){ for(j in 2:ploidy){ @@ -480,10 +506,10 @@ quickHaplo = function(nInd,nChr,segSites,genLen=1,ploidy=2L,inbred=FALSE){ #' @title Add segregating site to MapPop #' -#' @description This function allows for adding a new -#' segregating site with user supplied genotypes to a MapPop. +#' @description This function allows for adding a new +#' segregating site with user supplied genotypes to a MapPop. #' The position of the site is set using a genetic map position. -#' +#' #' @param mapPop an object of \code{\link{MapPop-class}} #' @param siteName name to give the segregating site #' @param chr which chromosome to add the site @@ -491,53 +517,53 @@ quickHaplo = function(nInd,nChr,segSites,genLen=1,ploidy=2L,inbred=FALSE){ #' @param haplo haplotypes for the site #' #' @return an object of \code{\link{MapPop-class}} -#' -#' @examples +#' +#' @examples #' # Creates a populations of 10 outbred individuals #' # Their genome consists of 1 chromosome and 2 segregating sites #' founderPop = quickHaplo(nInd=10,nChr=1,segSites=2) -#' +#' #' # Add a locus a the 0.5 Morgan map position #' haplo = matrix(sample(x=0:1, size=20, replace=TRUE), ncol=1) #' #' founderPop2 = addSegSite(founderPop, siteName="x", chr=1, mapPos=0.5, haplo=haplo) #' #' pullSegSiteHaplo(founderPop2) -#' +#' #' @export addSegSite = function(mapPop, siteName, chr, mapPos, haplo){ # Check validity of input data stopifnot(is(mapPop, "MapPop")) stopifnot(length(haplo)==(mapPop@nInd*mapPop@ploidy)) - + # Check that name isn't already present allSiteNames = unlist(unname(lapply(mapPop@genMap, names))) if(siteName%in%allSiteNames){ stop("The siteName '",siteName, "' is already in use.") } - + # Coerce haplo to a raw matrix haplo = matrix(as.raw(haplo), ncol=1) - + # Convert chr to a number, if needed if(is.character(chr)){ chr = match(chr, names(mapPop@genMap)) } - + # Extract temporary map and geno for target chromosome chrMap = mapPop@genMap[[chr]] markerNames = names(chrMap) chrGeno = mapPop@geno[chr] # Leaving in list format for getHaplo - + # Extract haplotype matrix from the bit array - haploMat = getHaplo(chrGeno, - mapPop@nLoci[chr], - 1:mapPop@nLoci[chr], + haploMat = getHaplo(chrGeno, + mapPop@nLoci[chr], + 1:mapPop@nLoci[chr], getNumThreads()) - + # Find position of insertion pos = findInterval(x=mapPos, vec=chrMap) - + # Insert site if(pos==length(chrMap)){ # Inserting site at the end @@ -552,24 +578,24 @@ addSegSite = function(mapPop, siteName, chr, mapPos, haplo){ haploMat = cbind(haplo, haploMat) }else{ # Inserting site in the middle - chrMap = c(chrMap[1:pos], mapPos, + chrMap = c(chrMap[1:pos], mapPos, chrMap[(pos+1L):length(chrMap)]) - names(chrMap) = c(markerNames[1:pos], siteName, + names(chrMap) = c(markerNames[1:pos], siteName, markerNames[(pos+1L):length(markerNames)]) - haploMat = cbind(haploMat[,1:pos,drop=FALSE], haplo, + haploMat = cbind(haploMat[,1:pos,drop=FALSE], haplo, haploMat[,(pos+1L):ncol(haploMat),drop=FALSE]) } - + # Convert haplotype matrix back to a bit array chrGeno = packHaplo(haploMat, ploidy=mapPop@ploidy, inbred=FALSE) - + # Return temporary map and geno mapPop@genMap[[chr]] = chrMap mapPop@geno[[chr]] = chrGeno - + # Add one to locus count mapPop@nLoci[chr] = mapPop@nLoci[chr] + 1L - + # Check if an inbred population has become outbred if(mapPop@inbred){ M = pullMarkerGeno(mapPop, markers=siteName) @@ -580,6 +606,6 @@ addSegSite = function(mapPop, siteName, chr, mapPos, haplo){ mapPop@inbred = FALSE } } - + return(mapPop) } diff --git a/R/hybrids.R b/R/hybrids.R index 239e0264..24c74b66 100644 --- a/R/hybrids.R +++ b/R/hybrids.R @@ -1,7 +1,7 @@ #' @title Hybrid crossing #' #' @description -#' A convience function for hybrid plant breeding simulations. Allows for +#' A convenient function for hybrid plant breeding simulations. Allows for #' easy specification of a test cross scheme and/or creation of an object #' of \code{\link{HybridPop-class}}. Note that the \code{\link{HybridPop-class}} #' should only be used if the parents were created using the \code{\link{makeDH}} @@ -10,7 +10,7 @@ #' #' @param females female population, an object of \code{\link{Pop-class}} #' @param males male population, an object of \code{\link{Pop-class}} -#' @param crossPlan either "testcross" for all possible combinantions +#' @param crossPlan either "testcross" for all possible combinations #' or a matrix with two columns for designed crosses #' @param returnHybridPop should results be returned as #' \code{\link{HybridPop-class}}. If false returns results as diff --git a/R/importData.R b/R/importData.R index 34056927..8393bbed 100644 --- a/R/importData.R +++ b/R/importData.R @@ -34,7 +34,7 @@ importGenMap = function(genMap){ names(genMap) = uniqueChr # Iterate through chromosomes - for(i in 1:length(uniqueChr)){ + for(i in seq_len(length(uniqueChr))){ take = (chromosome==uniqueChr[i]) tmpPos = position[take] @@ -155,7 +155,7 @@ importInbredGeno = function(geno, genMap, ped=NULL){ haplotypes = vector("list", length=length(genMap)) # Order haplotypes by chromosome - for(i in 1:length(genMap)){ + for(i in seq_len(length(genMap))){ mapMarkers = names(genMap[[i]]) take = match(mapMarkers, markerName) if(any(is.na(take))){ @@ -267,7 +267,7 @@ importHaplo = function(haplo, genMap, ploidy=2L, ped=NULL){ haplotypes = vector("list", length=length(genMap)) # Order haplotypes by chromosome - for(i in 1:length(genMap)){ + for(i in seq_len(length(genMap))){ mapMarkers = names(genMap[[i]]) take = match(mapMarkers, markerName) if(any(is.na(take))){ @@ -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]] = haplo[,take,drop=FALSE] } founderPop = newMapPop(genMap=genMap, diff --git a/R/mergePops.R b/R/mergePops.R index 6ba27777..9ea60fa2 100644 --- a/R/mergePops.R +++ b/R/mergePops.R @@ -18,13 +18,16 @@ #' #' #Create a list of populations and merge list #' pop = newPop(founderPop, simParam=SP) +#' pop@misc$tmp = rnorm(n=10) +#' pop@misc$tmp2 = rnorm(n=10) +#' #' popList = list(pop, pop) #' pop2 = mergePops(popList) #' #' @export mergePops = function(popList){ if(is(popList,"MultiPop")){ - for(i in 1:length(popList@pops)){ + for(i in seq_len(length(popList@pops))){ if(is(popList@pops[i],"MultiPop")){ popList@pops[i] = mergePops(popList@pops[i]) } @@ -40,16 +43,19 @@ mergePops = function(popList){ classes = classes[-remove] } stopifnot(all(classes=="Pop")) + #nChr nChr = do.call("c",lapply(popList, function(x) x@nChr)) stopifnot(all(nChr==nChr[1])) nChr = nChr[1] + #ploidy ploidy = do.call("c",lapply(popList, function(x) x@ploidy)) stopifnot(all(ploidy==ploidy[1])) ploidy = ploidy[1] + #nLoci nLoci = do.call("c",lapply(popList, function(x){ @@ -57,50 +63,82 @@ mergePops = function(popList){ })) stopifnot(all(nLoci)) nLoci = popList[[1]]@nLoci + #id id = do.call("c", lapply(popList, function(x) x@id)) + #iid iid = do.call("c", lapply(popList, function(x) x@iid)) + #mother mother = do.call("c", lapply(popList, function(x) x@mother)) + #father father= do.call("c", lapply(popList, function(x) x@father)) + #fixEff fixEff= do.call("c", lapply(popList, function(x) x@fixEff)) - + #misc - misc = do.call("c", - lapply(popList, - function(x) x@misc)) - + tmp = sapply(popList, function(x) length(x@misc)) + if(all(tmp == tmp[1]) & tmp[1]>0) { + tmp = lapply(popList, function(x) names(x@misc)) + allMatch = TRUE + if(length(tmp)>1){ + for(i in 2:length(tmp)){ + if(!all(tmp[[1]]==tmp[[i]])){ + allMatch = FALSE + break + } + } + } + if(allMatch){ + misc = vector("list", length=length(tmp[[1]])) + for(i in seq_len(length(tmp[[1]]))){ + miscTmp = lapply(popList, function(x) x@misc[[i]]) + misc[[i]] = do.call("c", miscTmp) + } + names(misc) = tmp[[1]] + }else{ + misc = list() + } + } else { + misc = list() + } + #sex sex = do.call("c", lapply(popList, function(x) x@sex)) + #nTraits nTraits = do.call("c",lapply(popList, function(x) x@nTraits)) stopifnot(all(nTraits==nTraits[1])) nTraits = nTraits[1] + #nInd nInd = do.call("c",lapply(popList, function(x) x@nInd)) + #gv gv = do.call("rbind",lapply(popList, function(x) x@gv)) + #pheno pheno = do.call("rbind",lapply(popList, function(x) x@pheno)) + #ebv ebv = do.call("c",lapply(popList, function(x) ncol(x@ebv))) @@ -110,10 +148,11 @@ mergePops = function(popList){ }else{ ebv = matrix(NA_real_,nrow=sum(nInd),ncol=0) } + #gxe if(nTraits>=1){ gxe = vector("list",length=nTraits) - for(trait in 1:nTraits){ + for(trait in seq_len(nTraits)){ if(!is.null(popList[[1]]@gxe[[trait]])){ tmp = lapply(popList,function(x) x@gxe[[trait]]) tmp = do.call("c",tmp) @@ -123,10 +162,13 @@ mergePops = function(popList){ }else{ gxe = list() } + #geno nBin = as.integer(nLoci%/%8L + (nLoci%%8L > 0L)) geno = mergeMultGeno(popList,nInd=nInd,nBin=nBin,ploidy=ploidy) dim(geno) = NULL # Account for matrix bug in RcppArmadillo + + #wrap it all up into a Pop nInd = sum(nInd) return(new("Pop", nInd=nInd, @@ -140,11 +182,11 @@ mergePops = function(popList){ mother=mother, father=father, fixEff=fixEff, + misc=misc, + miscPop=list(), nTraits=nTraits, gv=gv, gxe=gxe, pheno=pheno, - ebv=ebv, - misc=misc, - miscPop=list())) + ebv=ebv)) } diff --git a/R/misc.R b/R/misc.R index 01960c55..6f6d0ef0 100644 --- a/R/misc.R +++ b/R/misc.R @@ -18,13 +18,13 @@ convToImat = function(X){ #' SP <- SimParam$new(founderGenomes) #' SP$setSexes(sexes = "yes_sys") #' pop <- newPop(founderGenomes) -#' +#' #' isFemale(pop) #' isMale(pop) -#' +#' #' pop[isFemale(pop)] #' pop[isFemale(pop)]@sex -#' +#' #' @export isFemale <- function(x) { if (isPop(x)) { @@ -48,145 +48,35 @@ isMale <- function(x) { return(ret) } - -#' @rdname setMisc -#' @title Set miscelaneous information in a population -#' -#' @description Set miscelaneous information in a population -#' -#' @param x \code{\link{Pop-class}} -#' @param node character, name of the node to set within the \code{x@misc} slot -#' @param value, value to be saved into \code{x@misc[[*]][[node]]}; length of -#' \code{value} should be equal to \code{nInd(x)}; if its length is 1, then -#' it is repeated using \code{rep} (see examples) -#' -#' @details A \code{NULL} in \code{value} is ignored -#' -#' @return \code{\link{Pop-class}} -#' -#' @export -setMisc <- function(x, node = NULL, value = NULL) { - if (isPop(x)) { - if (is.null(node)) { - stop("Argument node must be provided!") - } - if (is.null(value)) { - stop("Argument value must be provided!") - } - n <- nInd(x) - if (length(value) == 1 && n > 1) { - value <- rep(x = value, times = n) - } - if (length(value) != n) { - stop("Argument value must be of length 1 or nInd(x)!") - } - for (ind in seq_len(n)) { - if (!is.null(value[ind])) { - x@misc[[ind]][node] <- value[ind] - } - } - } else { - stop("Argument x must be a Pop class object!") - } - return(x) -} - -#' @rdname getMisc -#' @title Get miscelaneous information in a population -#' -#' @description Get miscelaneous information in a population +#' @title Get pedigree #' -#' @param x \code{\link{Pop-class}} -#' @param node character, name of the node to get from the \code{x@misc} slot; -#' if \code{NULL} the whole \code{x@misc} slot is returned +#' @description +#' Returns the population's pedigree as stored in the +#' id, mother and father slots. NULL is returned if the +#' input population lacks the required. #' -#' @return The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]} +#' @param pop a population #' #' @examples -#' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) -#' SP <- SimParam$new(founderGenomes) -#' \dontshow{SP$nThreads = 1L} -#' basePop <- newPop(founderGenomes) -#' -#' basePop <- setMisc(basePop, node = "info", value = 1) -#' basePop@misc -#' getMisc(x = basePop, node = "info") -#' -#' basePop <- setMisc(basePop, node = "info2", value = c("A", "B", "C")) -#' basePop@misc -#' getMisc(x = basePop, node = "info2") -#' -#' n <- nInd(basePop) -#' location <- vector(mode = "list", length = n) -#' for (ind in seq_len(n)) { -#' location[[ind]] <- runif(n = 2, min = 0, max = 100) -#' } -#' location -#' basePop <- setMisc(basePop, node = "location", value = location) -#' basePop@misc -#' getMisc(x = basePop, node = "location") -#' -#' n <- nInd(basePop) -#' location <- vector(mode = "list", length = n) -#' for (ind in c(1, 3)) { -#' location[[ind]] <- runif(n = 2, min = 0, max = 100) -#' } -#' location -#' basePop <- setMisc(basePop, node = "location", value = location) -#' basePop@misc -#' getMisc(x = basePop, node = "location") -#' -#' getMisc(x = basePop) -#' -#' @export -getMisc <- function(x, node = NULL) { - if (isPop(x)) { - if (is.null(node)) { - ret <- x@misc - } else { - nInd <- nInd(x) - ret <- vector(mode = "list", length = nInd) - for (ind in seq_len(nInd)) { - if (!is.null(x@misc[[ind]][[node]])) { - ret[ind] <- x@misc[[ind]][node] - } - } - } - } else { - stop("Argument x must be a Pop class object!") - } - return(ret) -} - -#' @title Get pedigree -#' -#' @description -#' Returns the population's pedigree as stored in the -#' id, mother and father slots. NULL is returned if the -#' input population lacks the required. -#' -#' @param pop a population -#' -#' @examples #' # Create a founder population #' founderPop = quickHaplo(2,1,2) -#' +#' #' # Set simulation parameters #' SP = SimParam$new(founderPop) -#' +#' #' # Create a population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' # Get the pedigree #' getPed(pop) -#' +#' #' # Returns NULL when a population lacks a pedigree #' getPed(founderPop) -#' +#' #' @export getPed = function(pop){ - if(.hasSlot(pop, "id") & - .hasSlot(pop, "mother") & + if(.hasSlot(pop, "id") & + .hasSlot(pop, "mother") & .hasSlot(pop, "father")){ df = data.frame(id = pop@id, mother = pop@mother, @@ -199,60 +89,60 @@ getPed = function(pop){ #' @title Selection intensity -#' -#' @description +#' +#' @description #' Calculates the standardized selection intensity -#' +#' #' @param p the proportion of individuals selected -#' -#' @examples +#' +#' @examples #' selInt(0.1) -#' +#' #' @export selInt = function(p){ return(dnorm(qnorm(1-p))/p) } #' @title Calculate Smith-Hazel weights -#' +#' #' @description -#' Calculates weights for Smith-Hazel index given economice weights +#' Calculates weights for Smith-Hazel index given economice weights #' and phenotypic and genotypic variance-covariance matrices. -#' +#' #' @param econWt vector of economic weights #' @param varG the genetic variance-covariance matrix #' @param varP the phenotypic variance-covariance matrix -#' +#' #' @return a vector of weight for calculating index values -#' +#' #' @examples #' G = 1.5*diag(2)-0.5 #' E = diag(2) #' P = G+E #' wt = c(1,1) #' smithHazel(wt, G, P) -#' +#' #' @export smithHazel = function(econWt,varG,varP){ return(solve(varP)%*%varG%*%econWt) } #' @title Selection index -#' +#' #' @description -#' Calculates values of a selection index given trait values and -#' weights. This function is intended to be used in combination with -#' selection functions working on populations such as +#' Calculates values of a selection index given trait values and +#' weights. This function is intended to be used in combination with +#' selection functions working on populations such as #' \code{\link{selectInd}}. -#' +#' #' @param Y a matrix of trait values #' @param b a vector of weights #' @param scale should Y be scaled and centered -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} @@ -260,19 +150,19 @@ smithHazel = function(econWt,varG,varP){ #' G = 1.5*diag(2)-0.5 #Genetic correlation matrix #' SP$addTraitA(10, mean=c(0,0), var=c(1,1), corA=G) #' SP$setVarE(h2=c(0.5,0.5)) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Calculate Smith-Hazel weights #' econWt = c(1, 1) #' b = smithHazel(econWt, varG(pop), varP(pop)) -#' +#' #' #Selection 2 best individuals using Smith-Hazel index #' #selIndex is used as a trait -#' pop2 = selectInd(pop, nInd=2, trait=selIndex, +#' pop2 = selectInd(pop, nInd=2, trait=selIndex, #' simParam=SP, b=b) -#' +#' #' @export selIndex = function(Y,b,scale=FALSE){ if(scale){ @@ -282,40 +172,40 @@ selIndex = function(Y,b,scale=FALSE){ } #' @title Edit genome -#' +#' #' @description -#' Edits selected loci of selected individuals to a homozygous -#' state for either the 1 or 0 allele. The gv slot is recalculated to +#' Edits selected loci of selected individuals to a homozygous +#' state for either the 1 or 0 allele. The gv slot is recalculated to #' reflect the any changes due to editing, but other slots remain the same. -#' +#' #' @param pop an object of \code{\link{Pop-class}} #' @param ind a vector of individuals to edit -#' @param chr a vector of chromosomes to edit. +#' @param chr a vector of chromosomes to edit. #' Length must match length of segSites. -#' @param segSites a vector of segregating sites to edit. Length must +#' @param segSites a vector of segregating sites to edit. Length must #' match length of chr. #' @param allele either 0 or 1 for desired allele #' @param simParam an object of \code{\link{SimParam}} -#' +#' #' @return Returns an object of \code{\link{Pop-class}} -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' -#' #Change individual 1 to homozygous for the 1 allele +#' +#' #Change individual 1 to homozygous for the 1 allele #' #at locus 1, chromosome 1 -#' pop2 = editGenome(pop, ind=1, chr=1, segSites=1, +#' pop2 = editGenome(pop, ind=1, chr=1, segSites=1, #' allele=1, simParam=SP) -#' +#' #' @export editGenome = function (pop, ind, chr, segSites, allele, simParam = NULL) { if (is.null(simParam)) { @@ -339,7 +229,7 @@ editGenome = function (pop, ind, chr, segSites, allele, simParam = NULL) { BYTE = (segSites[i] - 1L)%/%8L + 1L BIT = (segSites[i] - 1L)%%8L + 1L for (selInd in ind) { - for (j in 1:pop@ploidy) { + for (j in seq_len(pop@ploidy)) { TMP = pop@geno[[selChr]][BYTE, j, selInd] TMP = rawToBits(TMP) TMP[BIT] = allele[i] @@ -358,36 +248,36 @@ editGenome = function (pop, ind, chr, segSites, allele, simParam = NULL) { } #' @title Edit genome - the top QTL -#' +#' #' @description -#' Edits the top QTL (with the largest additive effect) to a homozygous +#' Edits the top QTL (with the largest additive effect) to a homozygous #' state for the allele increasing. Only nonfixed QTL are edited The gv slot is #' recalculated to reflect the any changes due to editing, but other slots remain the same. -#' +#' #' @param pop an object of \code{\link{Pop-class}} #' @param ind a vector of individuals to edit #' @param nQtl number of QTL to edit #' @param trait which trait effects should guide selection of the top QTL #' @param increase should the trait value be increased or decreased #' @param simParam an object of \code{\link{SimParam}} -#' +#' #' @return Returns an object of \code{\link{Pop-class}} -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' -#' #Change up to 10 loci for individual 1 +#' +#' #Change up to 10 loci for individual 1 #' pop2 = editGenomeTopQtl(pop, ind=1, nQtl=10, simParam=SP) -#' +#' #' @export editGenomeTopQtl = function(pop, ind, nQtl, trait = 1, increase = TRUE, simParam = NULL) { if (is.null(simParam)) { @@ -397,7 +287,7 @@ editGenomeTopQtl = function(pop, ind, nQtl, trait = 1, increase = TRUE, simParam stopifnot(all(ind %in% (1:pop@nInd))) nQtl = as.integer(nQtl) stopifnot(nQtl > 0 & nQtl <= simParam$traits[[trait]]@nLoci) - + findTopQtl = function(pop, ind, nQtl, trait, increase, simParam) { # @title Find the top non fixed QTL for use in editGenome() # @param pop an object of \code{\link{Pop-class}} @@ -412,14 +302,14 @@ editGenomeTopQtl = function(pop, ind, nQtl, trait = 1, increase = TRUE, simParam # third indicating chromosome of the QTL # fourth indicates which allele we want to fix (edit to) QtlGeno = pullQtlGeno(pop=pop[ind],trait=trait,simParam=simParam) - + QtlEff = simParam$traits[[trait]]@addEff ret = vector(mode = "list", length = 4) ret[[1]] = ret[[2]] = ret[[3]] = ret[[4]] = rep(NA, times = nQtl) QtlEffRank = order(abs(QtlEff), decreasing = TRUE) nQtlInd = 0 Qtl = 0 - + while (nQtlInd < nQtl) { Qtl = Qtl + 1 if(Qtl>ncol(QtlGeno)){ @@ -455,15 +345,15 @@ editGenomeTopQtl = function(pop, ind, nQtl, trait = 1, increase = TRUE, simParam } } } - + # Locate QTL segsite to chromosomes tmp = cumsum(simParam$traits[[trait]]@lociPerChr) - for (Qtl in 1:nQtl) { + for (Qtl in seq_len(nQtl)) { ret[[3]][Qtl] = which(ret[[1]][Qtl] <= tmp)[1] } ret } - + for (ind2 in ind) { targetQtl = findTopQtl(pop = pop, ind = ind2, nQtl = nQtl, trait = trait, increase = increase, simParam = simParam) @@ -478,43 +368,43 @@ editGenomeTopQtl = function(pop, ind, nQtl, trait = 1, increase = TRUE, simParam } #' @title Usefulness criterion -#' +#' #' @description Calculates the usefulness criterion -#' -#' @param pop and object of \code{\link{Pop-class}} or +#' +#' @param pop and object of \code{\link{Pop-class}} or #' \code{\link{HybridPop-class}} -#' @param trait the trait for selection. Either a number indicating +#' @param trait the trait for selection. Either a number indicating #' a single trait or a function returning a vector of length nInd. #' @param use select on genetic values (\code{gv}, default), estimated -#' breeding values (\code{ebv}), breeding values (\code{bv}), +#' breeding values (\code{ebv}), breeding values (\code{bv}), #' or phenotypes (\code{pheno}) #' @param p the proportion of individuals selected -#' @param selectTop selects highest values if true. +#' @param selectTop selects highest values if true. #' Selects lowest values if false. #' @param simParam an object of \code{\link{SimParam}} -#' @param ... additional arguments if using a function for +#' @param ... additional arguments if using a function for #' trait -#' +#' #' @return Returns a numeric value -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' -#' #Determine usefulness of population +#' +#' #Determine usefulness of population #' usefulness(pop, simParam=SP) -#' +#' #' #Should be equivalent to GV of best individual #' max(gv(pop)) -#' +#' #' @export usefulness = function(pop,trait=1,use="gv",p=0.1, selectTop=TRUE,simParam=NULL,...){ @@ -529,59 +419,59 @@ usefulness = function(pop,trait=1,use="gv",p=0.1, } #' @title Linear transformation matrix -#' -#' @description -#' Creates an m by m linear transformation matrix that -#' can be applied to n by m uncorrelated deviates +#' +#' @description +#' Creates an m by m linear transformation matrix that +#' can be applied to n by m uncorrelated deviates #' sampled from a standard normal distribution to produce -#' correlated deviates with an arbitrary correlation -#' of R. If R is not positive semi-definite, the function +#' correlated deviates with an arbitrary correlation +#' of R. If R is not positive semi-definite, the function #' returns smoothing and returns a warning (see details). -#' +#' #' @param R a correlation matrix -#' -#' @details -#' An eigendecomposition is applied to the correlation -#' matrix and used to test if it is positive semi-definite. -#' If the matrix is not positive semi-definite, it is not a -#' valid correlation matrix. In this case, smoothing is -#' applied to the matrix (as described in the 'cor.smooth' of -#' the 'psych' library) to obtain a valid correlation matrix. -#' The resulting deviates will thus not exactly match the -#' desired correlation, but will hopefully be close if the -#' input matrix wasn't too far removed from a valid +#' +#' @details +#' An eigendecomposition is applied to the correlation +#' matrix and used to test if it is positive semi-definite. +#' If the matrix is not positive semi-definite, it is not a +#' valid correlation matrix. In this case, smoothing is +#' applied to the matrix (as described in the 'cor.smooth' of +#' the 'psych' library) to obtain a valid correlation matrix. +#' The resulting deviates will thus not exactly match the +#' desired correlation, but will hopefully be close if the +#' input matrix wasn't too far removed from a valid #' correlation matrix. -#' -#' @examples +#' +#' @examples #' # Create an 2x2 correlation matrix #' R = 0.5*diag(2) + 0.5 -#' -#' # Sample 1000 uncorrelated deviates from a +#' +#' # Sample 1000 uncorrelated deviates from a #' # bivariate standard normal distribution #' X = matrix(rnorm(2*1000), ncol=2) -#' +#' #' # Compute the transformation matrix #' T = transMat(R) -#' +#' #' # Apply the transformation to the deviates #' Y = X%*%T -#' +#' #' # Measure the sample correlation #' cor(Y) -#' +#' #' @export transMat = function(R){ - # Check if matrix is symmetric + # Check if matrix is symmetric # Stop if it is not nameR = deparse(substitute(R)) if(!isSymmetric(R)){ stop(nameR, " is not a symmetric matrix") } - + # Check if matrix is positive semi-definite # Provide a warning if it is not eig = eigen(R, symmetric=TRUE) - + if(min(eig$values)<.Machine$double.eps){ warning("Matrix is not positive semi-definite, see ?transMat for details") # Performing correlation matrix smoothing @@ -593,73 +483,73 @@ transMat = function(R){ newR = cov2cor(newR) eig = eigen(newR, symmetric=TRUE) } - + return( - t(eig$vectors %*% + t(eig$vectors %*% (t(eig$vectors)*sqrt(pmax(eig$values, 0))) ) ) } #' @title Add Random Mutations -#' +#' #' @description -#' Adds random mutations to individuals in a -#' population. Note that any existing phenotypes -#' or EBVs are kept. Thus, the user will need to run -#' \code{\link{setPheno}} and/or \code{\link{setEBV}} -#' to generate new phenotypes or EBVs that reflect +#' Adds random mutations to individuals in a +#' population. Note that any existing phenotypes +#' or EBVs are kept. Thus, the user will need to run +#' \code{\link{setPheno}} and/or \code{\link{setEBV}} +#' to generate new phenotypes or EBVs that reflect #' changes introduced by the new mutations. -#' +#' #' @param pop an object of \code{\link{Pop-class}} #' @param mutRate rate of new mutations #' @param returnPos should the positions of mutations be returned #' @param simParam an object of \code{\link{SimParam}} #' -#' @return an object of \code{\link{Pop-class}} if -#' returnPos=FALSE or a list containing a -#' \code{\link{Pop-class}} and a data.frame containing the +#' @return an object of \code{\link{Pop-class}} if +#' returnPos=FALSE or a list containing a +#' \code{\link{Pop-class}} and a data.frame containing the #' postions of mutations if returnPos=TRUE -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} #' SP$addTraitA(10) -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Introduce mutations #' pop = mutate(pop, simParam=SP) -#' +#' #' @export mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ # Mutation history variable IND=NULL; CHR=NULL; HAP=NULL; SITE=NULL - + # Number of haplotypes per chromosome nHap = pop@nInd*pop@ploidy - + # Number of total sites s = sum(pop@nLoci) - + # Number of mutations per haplotype nMut = rbinom(nHap, s, mutRate) - + if(any(nMut>0L)){ for(take in which(nMut>0L)){ # Determine haplotype and individual ind = (take-1L)%/%pop@ploidy + 1L hap = (take-1L)%%pop@ploidy + 1L - + # Sample mutation sites sites = sampleInt(nMut[take], s) + 1L - + # Resolve all mutations chr = 1L for(i in sites){ @@ -671,14 +561,14 @@ mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ break } } - + # Find site if(chr>1L){ - site = i - sum(pop@nLoci[1L:(chr-1L)]) + site = i - sum(pop@nLoci[1L:(chr-1L)]) }else{ site = i } - + # Create mutation BYTE = (site-1L)%/%8L + 1L BIT = (site-1L)%%8L + 1L @@ -687,7 +577,7 @@ mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ TMP[BIT] = ifelse(TMP[BIT], as.raw(0L), as.raw(1L)) TMP = packBits(TMP) pop@geno[[chr]][BYTE,hap,ind] = TMP - + # Record results if(returnPos){ IND = c(IND, ind) @@ -697,7 +587,7 @@ mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ } } } - + # Reset population PHENO = pop@pheno EBV = pop@ebv @@ -705,7 +595,7 @@ mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ pop@pheno = PHENO pop@ebv = EBV } - + # Return results if(returnPos){ return(list(pop,data.frame(individual=IND,chromosome=CHR,haplotype=HAP,site=SITE))) @@ -715,32 +605,32 @@ mutate = function(pop, mutRate=2.5e-8, returnPos=FALSE, simParam=NULL){ } #' @title Lose individuals at random -#' +#' #' @description -#' Samples individuals at random to remove from the population. -#' The user supplies a probability for the individuals to be +#' Samples individuals at random to remove from the population. +#' The user supplies a probability for the individuals to be #' removed from the population. -#' +#' #' @param pop an object of \code{\link{Pop-class}} #' @param p the expected proportion of individuals that will #' be lost to attrition. #' #' @return an object of \code{\link{Pop-class}} -#' -#' @examples +#' +#' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=100, nChr=1, segSites=10) -#' +#' #' #Set simulation parameters #' SP = SimParam$new(founderPop) #' \dontshow{SP$nThreads = 1L} -#' +#' #' #Create population #' pop = newPop(founderPop, simParam=SP) -#' +#' #' #Lose an expected 5% of individuals #' pop = attrition(pop, p=0.05) -#' +#' #' @export attrition = function(pop, p){ take = as.logical(rbinom(pop@nInd, size=1, prob=1-p)) @@ -758,7 +648,7 @@ rnormWithSeed = function(n, u){ if(is.null(origSeed)){ rm(list =".Random.seed", envir=glbEnv) }else{ - assign(".Random.seed", value=origSeed, + assign(".Random.seed", value=origSeed, envir=glbEnv) } }) diff --git a/R/phenotypes.R b/R/phenotypes.R index 197e1873..19f73e09 100644 --- a/R/phenotypes.R +++ b/R/phenotypes.R @@ -33,7 +33,7 @@ calcPheno = function(pop, varE, reps, p, traits, simParam){ } gv = pop@gv - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ if(.hasSlot(simParam$traits[[traits[i]]], "envVar")){ stdDev = sqrt(simParam$traits[[traits[i]]]@envVar) gv[,traits[i]] = gv[,traits[i]] + @@ -183,7 +183,7 @@ setPheno = function(pop, h2=NULL, H2=NULL, varE=NULL, corE=NULL, all(varA>0), all(varG>0)) varE = numeric(nTraits) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ tmp = varA[i]/h2[i]-varG[i] if(tmp<0){ stop(paste0("h2=",h2[i]," is not possible for trait ",traits[i])) @@ -198,7 +198,7 @@ setPheno = function(pop, h2=NULL, H2=NULL, varE=NULL, corE=NULL, stopifnot(length(H2)==nTraits) varE = numeric(nTraits) - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ tmp = varG[i]/H2[i]-varG[i] varE[i] = tmp } diff --git a/R/polyploids.R b/R/polyploids.R index e60b7142..64536898 100644 --- a/R/polyploids.R +++ b/R/polyploids.R @@ -51,7 +51,7 @@ reduceGenome = function(pop,nProgeny=1,useFemale=TRUE,keepParents=TRUE, }else{ # Create dummy map with zero genetic distance map = vector("list",pop@nChr) - for(i in 1:pop@nChr){ + for(i in seq_len(pop@nChr)){ map[[i]] = rep(0,pop@nLoci[i]) } map = as.matrix(map) @@ -142,7 +142,7 @@ doubleGenome = function(pop, keepParents=TRUE, } geno = pop@geno - for(i in 1:pop@nChr){ + for(i in seq_len(pop@nChr)){ geno[[i]] = geno[[i]][,rep(1:pop@ploidy,each=2),,drop=FALSE] } @@ -156,7 +156,7 @@ doubleGenome = function(pop, keepParents=TRUE, if(simParam$isTrackRec){ # Match haplotypes according to original ploidy hist = vector("list", pop@ploidy) - for(i in 1:pop@ploidy){ + for(i in seq_len(pop@ploidy)){ hist[[i]] = cbind(i, 1L) } @@ -248,12 +248,12 @@ mergeGenome = function(females,males,crossPlan,simParam=NULL){ # Merge genotype data geno = vector("list", females@nChr) - for(i in 1:females@nChr){ + for(i in seq_len(females@nChr)){ geno[[i]] = array(as.raw(0), dim = c(dim(females@geno[[i]])[1], females@ploidy+males@ploidy, nrow(crossPlan))) - for(j in 1:nrow(crossPlan)){ + for(j in seq_len(nrow(crossPlan))){ # Add female gametes geno[[i]][,1:females@ploidy,j] = females@geno[[i]][,,crossPlan[j,1]] @@ -277,12 +277,12 @@ mergeGenome = function(females,males,crossPlan,simParam=NULL){ hist = vector("list", rPop@ploidy) # Add female contribution - for(i in 1:females@ploidy){ + for(i in seq_len(females@ploidy)){ hist[[i]] = cbind(i, 1L) } # Add male contribution - for(i in 1:males@ploidy){ + for(i in seq_len(males@ploidy)){ hist[[i+females@ploidy]] = cbind(i, 1L) } diff --git a/R/popSummary.R b/R/popSummary.R index e2b75033..79d17fd5 100644 --- a/R/popSummary.R +++ b/R/popSummary.R @@ -1,8 +1,3 @@ -# Internal function for calculating mean EBV of populations -# Used selectPop MultiPop-class -meanEBV = function(pop){ - colMeans(pop@ebv) -} #' @title Mean genetic values #' @@ -54,6 +49,33 @@ meanP = function(pop){ colMeans(pop@pheno) } +#' @title Mean estimated breeding values +#' +#' @description Returns the mean estimated breeding values for all traits +#' +#' @param pop an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}} +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitA(10) +#' trtH2 = 0.5 +#' SP$setVarE(h2=trtH2) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' pop@ebv = trtH2 * (pop@pheno - meanP(pop)) #ind performance based EBV +#' meanEBV(pop) +#' +#' @export +meanEBV = function(pop){ + colMeans(pop@ebv) +} + #' @title Total genetic variance #' #' @description Returns total genetic variance for all traits @@ -108,6 +130,36 @@ varP = function(pop){ return(P) } +#' @title Variance of estimated breeding values +#' +#' @description Returns variance of estimated breeding values for all traits +#' +#' @param pop an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}} +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitA(10) +#' trtH2 = 0.5 +#' SP$setVarE(h2=trtH2) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' pop@ebv = trtH2 * (pop@pheno - meanP(pop)) #ind performance based EBV +#' varA(pop) +#' varEBV(pop) +#' +#' @export +varEBV = function(pop){ + ebv = popVar(pop@ebv) + rownames(ebv) = colnames(ebv) = colnames(pop@ebv) + return(ebv) +} + #' @title Sumarize genetic parameters #' #' @description @@ -148,6 +200,8 @@ varP = function(pop){ #' \item{gv_a}{a matrix of additive genetic values with dimensions nInd by nTraits} #' \item{gv_d}{a matrix of dominance genetic values with dimensions nInd by nTraits} #' \item{gv_aa}{a matrix of additive-by-additive genetic values with dimensions nInd by nTraits} +#' \item{alpha}{a list of average allele subsitution effects with length nTraits} +#' \item{alpha_HW}{a list of average allele subsitution effects at Hardy-Weinberg equilibrium with length nTraits} #' } #' #' @examples @@ -169,7 +223,7 @@ genParam = function(pop,simParam=NULL){ if(is.null(simParam)){ simParam = get("SP",envir=.GlobalEnv) } - stopifnot(class(pop)=="Pop") + nInd = nInd(pop) nTraits = simParam$nTraits traitNames = simParam$traitNames @@ -186,8 +240,13 @@ genParam = function(pop,simParam=NULL){ covG_HW = mu = mu_HW = gv_mu = covAAA_L = covDAA_L = covAD_L = genicVarA + # Average effect of an allele substitution + alpha = vector("list", length=nTraits) + names(alpha) = traitNames + alpha_HW = alpha + #Loop through trait calculations - for(i in 1:nTraits){ + for(i in seq_len(nTraits)){ trait = simParam$traits[[i]] tmp = calcGenParam(trait,pop,simParam$nThreads) genicVarA[i] = tmp$genicVarA2 @@ -229,6 +288,8 @@ genParam = function(pop,simParam=NULL){ covAAA_L[i] = popVar(cbind(bv[,i],aa[,i]))[1,2] covDAA_L[i] = popVar(cbind(dd[,i],aa[,i]))[1,2] } + alpha[[i]] = tmp$alpha + alpha_HW[[i]] = tmp$alpha_HW } varA = popVar(bv) @@ -274,7 +335,9 @@ genParam = function(pop,simParam=NULL){ gv_mu=gv_mu, gv_a=gv_a, gv_d=gv_d, - gv_aa=gv_aa) + gv_aa=gv_aa, + alpha=alpha, + alpha_HW=alpha_HW) return(output) } diff --git a/R/pullGeno.R b/R/pullGeno.R index fafbf8b1..07480821 100644 --- a/R/pullGeno.R +++ b/R/pullGeno.R @@ -10,7 +10,7 @@ selectLoci = function(chr, inLociPerChr, inLociLoc){ outLociPerChr[chr] = inLociPerChr[chr] outLociLoc = numeric(sum(outLociPerChr)) inStart = outStart = inEnd = outEnd = 0L - for(i in 1:nChr){ + for(i in seq_len(nChr)){ inStart = inStart + 1L inEnd = inEnd + inLociPerChr[i] if(outLociPerChr[i]>0){ @@ -32,7 +32,7 @@ selectLoci = function(chr, inLociPerChr, inLociLoc){ getLociNames = function(lociPerChr, lociLoc, genMap){ lociNames = character(length(lociLoc)) start = end = 0L - for(chr in 1:length(lociPerChr)){ + for(chr in seq_len(length(lociPerChr))){ if(lociPerChr[chr]>0){ start = end + 1L end = end + lociPerChr[chr] @@ -54,7 +54,7 @@ mapLoci = function(markers, genMap){ lociLoc = vector("list", length(genMap)) # Loop through chromosomes - for(i in 1:length(genMap)){ + for(i in seq_len(length(genMap))){ # Initialize lociLoc lociLoc[[i]] = integer() diff --git a/R/selection.R b/R/selection.R index 2f2133f6..9d70d09c 100644 --- a/R/selection.R +++ b/R/selection.R @@ -1,53 +1,70 @@ # Returns a vector response from a population # pop is an object of class Pop or HybridPop # trait is a vector of traits or a function -# use is "rand", "gv", "ebv", "pheno", or "bv" -# "bv" doesn't work on class HybridPop +# (must work like trait(pop@use, ...) or trait(use(pop, ...), ...)). +# See examples in \code{selectInd}. +# use is a character ("rand", "gv", "ebv", "pheno", or "bv"; note that +# "bv" doesn't work on class HybridPop) +# or a function (must work like use(pop, trait, ...)). +# See examples in \code{selectInd}. # simParam is an object of class SimParam, it is only called when use="bv" # ... are additional arguments passed to trait when trait is a function getResponse = function(pop,trait,use,simParam=NULL,...){ - use = tolower(use) - if(use=="rand"){ - return(rnorm(pop@nInd)) - } if(is(trait,"function")){ - if(use=="gv"){ - response = trait(pop@gv,...) - }else if(use=="ebv"){ - response = trait(pop@ebv,...) - }else if(use=="pheno"){ - response = trait(pop@pheno,...) - }else if(use=="bv"){ - if(is(pop,"HybridPop")){ - stop("Use='bv' is not a valid option for HybridPop") + if(is.character(use)){ + use = tolower(use) + if(use=="rand"){ + return(rnorm(pop@nInd)) + }else if(use=="gv"){ + response = trait(pop@gv,...) + }else if(use=="ebv"){ + response = trait(pop@ebv,...) + }else if(use=="pheno"){ + response = trait(pop@pheno,...) + }else if(use=="bv"){ + if(is(pop,"HybridPop")){ + stop("Use='bv' is not a valid option for HybridPop") + } + response = genParam(pop,simParam=simParam)$bv + response = trait(response,...) + }else{ + stop(paste0("Use=",use," is not an option")) } - response = genParam(pop,simParam=simParam)$bv - response = trait(response,...) + }else if(is(use,"function")){ + response = trait(use(pop, ...), ...) }else{ - stop(paste0("Use=",use," is not an option")) + stop("use must be a character or a function") } - }else{ - if(is.character(trait)){ - # Suspect trait is a name + }else{ # trait is not a function, so must be numeric or character + if(is.character(trait)){ # Suspect trait is a name take = match(trait, simParam$traitNames) if(is.na(take)){ stop("'",trait,"' did not match any trait names") } trait = take } - if(use == "gv"){ - response = pop@gv[,trait,drop=FALSE] - }else if(use=="ebv"){ - response = pop@ebv[,trait,drop=FALSE] - }else if(use=="pheno"){ - response = pop@pheno[,trait,drop=FALSE] - }else if(use=="bv"){ - if(is(pop,"HybridPop")){ - stop("Use='bv' is not a valid option for HybridPop") + if(is.character(use)){ + use = tolower(use) + if(use=="rand"){ + return(rnorm(pop@nInd)) + }else if(use == "gv"){ + response = pop@gv[,trait,drop=FALSE] + }else if(use=="ebv"){ + response = pop@ebv[,trait,drop=FALSE] + }else if(use=="pheno"){ + response = pop@pheno[,trait,drop=FALSE] + }else if(use=="bv"){ + if(is(pop,"HybridPop")){ + stop("Use='bv' is not a valid option for HybridPop") + } + response = genParam(pop,simParam=simParam)$bv[,trait,drop=FALSE] + }else{ + stop(paste0("Use=",use," is not an option")) } - response = genParam(pop,simParam=simParam)$bv[,trait,drop=FALSE] + }else if(is(use,"function")){ + response = use(pop, trait=trait, ...) }else{ - stop(paste0("Use=",use," is not an option")) + stop("use must be a character or a function") } } if(any(is.na(response))){ @@ -84,15 +101,17 @@ checkSexes = function(pop,sex,simParam,...){ if(simParam$sexes=="no"){ return(eligible) }else{ + # Check in gender is incorrectly being used + args = list(...) + if(any(names(args)=="gender")){ + stop("The discontinued 'gender' argument appears to be in use. This argument was renamed as 'sex' in AlphaSimR version 0.13.0.") + } if(sex=="B"){ - # Check in gender is incorrectly being used - args = list(...) - if(any(names(args)=="gender")){ - stop("The discontinued 'gender' argument appears to be in use. This argument was renamed as 'sex' in AlphaSimR version 0.13.0.") - } return(eligible) - }else{ + }else if(sex=="F" | sex=="M"){ return(eligible[pop@sex%in%sex]) + }else{ + stop("Invalid option supplied for 'sex'") } } } @@ -117,27 +136,32 @@ getFam = function(pop,famType){ #' population. #' #' @param pop and object of \code{\link{Pop-class}}, -#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} +#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} #' @param nInd the number of individuals to select #' @param trait the trait for selection. Either a number indicating -#' a single trait or a function returning a vector of length nInd. -#' The function must work on a vector or matrix of \code{use} values. -#' See the examples and \code{\link{selIndex}}. -#' @param use select on genetic values "gv", estimated -#' breeding values "ebv", breeding values "bv", phenotypes "pheno", -#' or randomly "rand" +#' a single trait or a function returning a vector of length nInd. +#' The function must work on a vector or matrix of \code{use} values as +#' \code{trait(pop@use, ...)} - depending on what \code{use} is. +#' See the examples and \code{\link{selIndex}}. +#' @param use the selection criterion. Either a character +#' (genetic values "gv", estimated breeding values "ebv", breeding values "bv", +#' phenotypes "pheno", or randomly "rand") or +#' a function returning a vector of length nInd. +#' The function must work on \code{pop} as \code{use(pop, trait, ...)} or +#' as \code{trait(pop@use, ...)} depending on what \code{trait} is. +#' See the examples. #' @param sex which sex to select. Use "B" for both, "F" for -#' females and "M" for males. If the simulation is not using sexes, -#' the argument is ignored. +#' females and "M" for males. If the simulation is not using sexes, +#' the argument is ignored. #' @param selectTop selects highest values if true. -#' Selects lowest values if false. +#' Selects lowest values if false. #' @param returnPop should results be returned as a -#' \code{\link{Pop-class}}. If FALSE, only the index of selected -#' individuals is returned. +#' \code{\link{Pop-class}}. If FALSE, only the index of selected +#' individuals is returned. #' @param candidates an optional vector of eligible selection candidates. #' @param simParam an object of \code{\link{SimParam}} #' @param ... additional arguments if using a function for -#' trait +#' \code{trait} or \code{use} #' #' @return Returns an object of \code{\link{Pop-class}}, #' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} @@ -157,19 +181,31 @@ getFam = function(pop,famType){ #' #' #Select top 5 (directional selection) #' pop2 = selectInd(pop, 5, simParam=SP) -#' hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -#' abline(v = pop2@pheno, col = "red", lwd = 2) +#' hist(pop@pheno); abline(v=pop@pheno, lwd=2) +#' abline(v=pop2@pheno, col="red", lwd=2) #' #' #Select 5 most deviating from an optima (disruptive selection) -#' squaredDeviation = function(x, optima = 0) (x - optima)^2 -#' pop3 = selectInd(pop, 5, simParam=SP, trait = squaredDeviation, selectTop = TRUE) -#' hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -#' abline(v = pop3@pheno, col = "red", lwd = 2) +#' squaredDeviation = function(x, optima=0) (x - optima)^2 +#' pop3 = selectInd(pop, 5, trait=squaredDeviation, selectTop=TRUE, simParam=SP) +#' hist(pop@pheno); abline(v=pop@pheno, lwd=2) +#' abline(v=pop3@pheno, col="red", lwd=2) #' #' #Select 5 least deviating from an optima (stabilising selection) -#' pop4 = selectInd(pop, 5, simParam=SP, trait = squaredDeviation, selectTop = FALSE) -#' hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -#' abline(v = pop4@pheno, col = "red", lwd = 2) +#' pop4 = selectInd(pop, 5, trait=squaredDeviation, selectTop=FALSE, simParam=SP) +#' hist(pop@pheno); abline(v=pop@pheno, lwd=2) +#' abline(v=pop4@pheno, col="red", lwd=2) +#' +#' #Select 5 individuals based on miscelaneous information with use function +#' pop@misc = list(smth=rnorm(10), smth2=rnorm(10)) +#' useFunc = function(pop, trait=NULL) pop@misc$smth + pop@misc$smth2 +#' pop5 = selectInd(pop, 5, use=useFunc, simParam=SP) +#' pop5@id +#' +#' #... equivalent result with the use & trait function +#' useFunc2 = function(pop, trait=NULL) cbind(pop@misc$smth, pop@misc$smth2) +#' trtFunc = function(x) rowSums(x) +#' pop6 = selectInd(pop, 5, trait=trtFunc, use=useFunc2, simParam=SP) +#' pop6@id #' #' @export selectInd = function(pop,nInd,trait=1,use="pheno",sex="B", @@ -216,30 +252,35 @@ selectInd = function(pop,nInd,trait=1,use="pheno",sex="B", #' population. #' #' @param pop and object of \code{\link{Pop-class}}, -#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} +#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} #' @param nFam the number of families to select #' @param trait the trait for selection. Either a number indicating -#' a single trait or a function returning a vector of length nInd. -#' The function must work on a vector or matrix of \code{use} values. -#' See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}. -#' @param use select on genetic values "gv", estimated -#' breeding values "ebv", breeding values "bv", phenotypes "pheno", -#' or randomly "rand" +#' a single trait or a function returning a vector of length nInd. +#' The function must work on a vector or matrix of \code{use} values as +#' \code{trait(pop@use, ...)} - depending on what \code{use} is. +#' See the examples and \code{\link{selIndex}}. +#' @param use the selection criterion. Either a character +#' (genetic values "gv", estimated breeding values "ebv", breeding values "bv", +#' phenotypes "pheno", or randomly "rand") or +#' a function returning a vector of length nInd. +#' The function must work on \code{pop} as \code{use(pop, trait, ...)} or +#' as \code{trait(pop@use, ...)} depending on what \code{trait} is. +#' See the examples. #' @param sex which sex to select. Use "B" for both, "F" for -#' females and "M" for males. If the simulation is not using sexes, -#' the argument is ignored. +#' females and "M" for males. If the simulation is not using sexes, +#' the argument is ignored. #' @param famType which type of family to select. Use "B" for -#' full-sib families, "F" for half-sib families on female side and "M" -#' for half-sib families on the male side. +#' full-sib families, "F" for half-sib families on female side and "M" +#' for half-sib families on the male side. #' @param selectTop selects highest values if true. -#' Selects lowest values if false. +#' Selects lowest values if false. #' @param returnPop should results be returned as a -#' \code{\link{Pop-class}}. If FALSE, only the index of selected -#' individuals is returned. +#' \code{\link{Pop-class}}. If FALSE, only the index of selected +#' individuals is returned. #' @param candidates an optional vector of eligible selection candidates. #' @param simParam an object of \code{\link{SimParam}} #' @param ... additional arguments if using a function for -#' trait +#' \code{trait} and \code{use} #' #' @return Returns an object of \code{\link{Pop-class}}, #' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} @@ -318,30 +359,35 @@ selectFam = function(pop,nFam,trait=1,use="pheno",sex="B", #' from a full-sib family if it has less than or equal to nInd individuals. #' #' @param pop and object of \code{\link{Pop-class}}, -#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} +#' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} #' @param nInd the number of individuals to select within a family #' @param trait the trait for selection. Either a number indicating -#' a single trait or a function returning a vector of length nInd. -#' The function must work on a vector or matrix of \code{use} values. -#' See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}. -#' @param use select on genetic values "gv", estimated -#' breeding values "ebv", breeding values "bv", phenotypes "pheno", -#' or randomly "rand" +#' a single trait or a function returning a vector of length nInd. +#' The function must work on a vector or matrix of \code{use} values as +#' \code{trait(pop@use, ...)} - depending on what \code{use} is. +#' See the examples and \code{\link{selIndex}}. +#' @param use the selection criterion. Either a character +#' (genetic values "gv", estimated breeding values "ebv", breeding values "bv", +#' phenotypes "pheno", or randomly "rand") or +#' a function returning a vector of length nInd. +#' The function must work on \code{pop} as \code{use(pop, trait, ...)} or +#' as \code{trait(pop@use, ...)} depending on what \code{trait} is. +#' See the examples. #' @param sex which sex to select. Use "B" for both, "F" for -#' females and "M" for males. If the simulation is not using sexes, -#' the argument is ignored. +#' females and "M" for males. If the simulation is not using sexes, +#' the argument is ignored. #' @param famType which type of family to select. Use "B" for -#' full-sib families, "F" for half-sib families on female side and "M" -#' for half-sib families on the male side. +#' full-sib families, "F" for half-sib families on female side and "M" +#' for half-sib families on the male side. #' @param selectTop selects highest values if true. -#' Selects lowest values if false. +#' Selects lowest values if false. #' @param returnPop should results be returned as a -#' \code{\link{Pop-class}}. If FALSE, only the index of selected -#' individuals is returned. +#' \code{\link{Pop-class}}. If FALSE, only the index of selected +#' individuals is returned. #' @param candidates an optional vector of eligible selection candidates. #' @param simParam an object of \code{\link{SimParam}} #' @param ... additional arguments if using a function for -#' trait +#' \code{trait} and \code{use} #' #' @return Returns an object of \code{\link{Pop-class}}, #' \code{\link{HybridPop-class}} or \code{\link{MultiPop-class}} @@ -425,25 +471,30 @@ selectWithinFam = function(pop,nInd,trait=1,use="pheno",sex="B", #' selection as occuring before or after pollination. #' #' @param pop and object of \code{\link{Pop-class}} -#' or \code{\link{MultiPop-class}} +#' or \code{\link{MultiPop-class}} #' @param nInd the number of plants to select #' @param nSeeds number of seeds per plant #' @param probSelf percentage of seeds expected from selfing. -#' Value ranges from 0 to 1. +#' Value ranges from 0 to 1. #' @param pollenControl are plants selected before pollination #' @param trait the trait for selection. Either a number indicating -#' a single trait or a function returning a vector of length nInd. -#' The function must work on a vector or matrix of \code{use} values. -#' See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}. -#' @param use select on genetic values "gv", estimated -#' breeding values "ebv", breeding values "bv", phenotypes "pheno", -#' or randomly "rand" +#' a single trait or a function returning a vector of length nInd. +#' The function must work on a vector or matrix of \code{use} values as +#' \code{trait(pop@use, ...)} - depending on what \code{use} is. +#' See the examples and \code{\link{selIndex}}. +#' @param use the selection criterion. Either a character +#' (genetic values "gv", estimated breeding values "ebv", breeding values "bv", +#' phenotypes "pheno", or randomly "rand") or +#' a function returning a vector of length nInd. +#' The function must work on \code{pop} as \code{use(pop, trait, ...)} or +#' as \code{trait(pop@use, ...)} depending on what \code{trait} is. +#' See the examples. #' @param selectTop selects highest values if true. -#' Selects lowest values if false. +#' Selects lowest values if false. #' @param candidates an optional vector of eligible selection candidates. #' @param simParam an object of \code{\link{SimParam}} #' @param ... additional arguments if using a function for -#' trait +#' \code{trait} and \code{use} #' #' @return Returns an object of \code{\link{Pop-class}} #' or \code{\link{MultiPop-class}} diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 17c81c09..aea8e1dd 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -1,7 +1,7 @@ ### For Introduction @article{MaCS, - author = {Gary K. Chen and Paul Marjoram and Jeffery D. Wall}, + author = {Chen, Gary K. and Marjoram, Paul and Wall, Jeffery D.}, title = {Fast and Flexible Simulation of DNA Sequence Data}, journal = {Genome Research}, volume = {19}, @@ -11,7 +11,7 @@ @article{MaCS } @article{AlphaSim, - author = {Anne-Michelle Faux and Gregor Gorjanc and R. Chris Gaynor and Mara Battagin and Stefan M. Edwards and David L. Wilson and Sarah J. Hearne and Serap Gonen and John M. Hickey}, + author = {Faux, Anne-Michelle andGorjanc, Gregor and Gaynor, R. Chris and Battagin, Mara and Edwards, Stefan M. and Wilson, David L. and Hearne, Sarah J. and Gonen, Serap and Hickey, John M.}, title = {AlphaSim: Software for Breeding Program Simulation}, journal = {The Plant Genome}, volume = {9}, @@ -24,7 +24,7 @@ @article{AlphaSim @book{Gallais, title = {Quantitative Genetics and Breeding Methods in Autopolyploid Plants}, - author = {A. Gallais}, + author = {Gallais, A.}, year = {2003}, publisher = {INRA}, address = {Paris}, @@ -33,7 +33,7 @@ @book{Gallais @book{FalconerMackay, title = {Introduction to Quantitative Genetics}, - author = {D. S. Falconer and Trudy F. C. Mackay}, + author = {Falconer, D. S. and Mackay, Trudy F. C.}, edition = {4}, year = {1996}, publisher = {Longman}, @@ -88,7 +88,7 @@ @article{Falconer_85 } @article{Massman13, - author = {Jon M. Massman and Andres Gordillo and Robenzon E. Lorenzana and Rex Bernardo}, + author = {Massman, Jon M. and Gordillo, Andres and Lorenzana, Robenzon E. and Bernardo, Rex}, title = {Genomewide Predictions from Maize Single-Cross Data}, journal = {Theoretical and Applied Genetics}, volume = {126}, @@ -100,7 +100,7 @@ @article{Massman13 @manual{EMMREML, title = {EMMREML: Fitting Mixed Models with Known Covariance Structures}, - author = {Deniz Akdemir and Okeke Uche Godfrey}, + author = {Akdemir, Deniz and Godfrey, Okeke Uche}, year = {2015}, note = {R package version 3.1}, url = {https://CRAN.R-project.org/package=EMMREML}, @@ -108,7 +108,7 @@ @manual{EMMREML @article{sommer, title = {Genome assisted prediction of quantitative traits using the R package sommer}, - author = {Giovanny Covarrubias-Pazaran}, + author = {Covarrubias-Pazaran, Giovanny}, journal = {PLoS ONE}, year = {2016}, volume = {11}, @@ -131,29 +131,29 @@ @article{xiang_2016 pages = {92} } -@article{fisher_1919, - title={XV.—The Correlation between Relatives on the Supposition of Mendelian Inheritance.}, - volume={52}, - DOI={10.1017/S0080456800012163}, - number={2}, - journal={Transactions of the Royal Society of Edinburgh}, - publisher={Royal Society of Edinburgh Scotland Foundation}, - author={Fisher, R. A.}, - year={1919}, +@article{fisher_1919, + title={XV.—The Correlation between Relatives on the Supposition of Mendelian Inheritance.}, + volume={52}, + DOI={10.1017/S0080456800012163}, + number={2}, + journal={Transactions of the Royal Society of Edinburgh}, + publisher={Royal Society of Edinburgh Scotland Foundation}, + author={Fisher, R. A.}, + year={1919}, pages={399–433} } ### MaCS histories -@article{cattle, - title={Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence, with Correction for Sequence Errors}, - volume={30}, +@article{cattle, + title={Inferring Demography from Runs of Homozygosity in Whole-Genome Sequence, with Correction for Sequence Errors}, + volume={30}, issue={9}, - DOI={https://doi.org/10.1093/molbev/mst125}, - journal={Molecular Biology and Evolution}, - publisher={Royal Society of Edinburgh Scotland Foundation}, - author={ Iona M. MacLeod and Denis M. Larkin and Harris A. Lewin and Ben J. Hayes and Mike E. Goddard}, - year={2013}, + DOI={https://doi.org/10.1093/molbev/mst125}, + journal={Molecular Biology and Evolution}, + publisher={Royal Society of Edinburgh Scotland Foundation}, + author={MacLeod, Iona M. and Larkin, Denis M. and Lewin, Harris A.and Hayes, Ben J. and Goddard, Mike E.}, + year={2013}, month=sept, pages={2209–2223} } @@ -165,7 +165,7 @@ @article{hickey_2014 doi = {10.2135/cropsci2013.03.0195}, language = {en}, journal = {Crop Science}, - author = {John M. Hickey and Susanne Dreisigacker and Jose Crossa and Sarah Hearne and Raman Babu and Boddupalli M. Prasanna and Martin Grondona and Andres Zambelli and Vanessa S. Windhausen and Ky Mathews and Gregor Gorjanc}, + author = {Hickey, John M.and Dreisigacker, Susanne and Crossa, Jose and Hearne, Sarah and Babu, Raman and Prasanna, Boddupalli M. and Grondona, Martin and Zambelli, Andres and Windhausen, Vanessa S. and Mathews, Ky and Gorjanc, Gregor}, year = {2014}, pages = {1476-1488}, -} \ No newline at end of file +} diff --git a/man/MultiPop-class.Rd b/man/MultiPop-class.Rd index a9f29013..52d5bbf1 100644 --- a/man/MultiPop-class.Rd +++ b/man/MultiPop-class.Rd @@ -6,6 +6,7 @@ \alias{[,MultiPop-method} \alias{[[,MultiPop-method} \alias{c,MultiPop-method} +\alias{length,MultiPop-method} \alias{isMultiPop} \title{Multi-Population} \usage{ @@ -15,6 +16,8 @@ \S4method{c}{MultiPop}(x, ...) +\S4method{length}{MultiPop}(x) + isMultiPop(x) } \arguments{ @@ -36,6 +39,8 @@ It is designed to behave like a list of populations. \item \code{c(MultiPop)}: Combine multiple MultiPops +\item \code{length(MultiPop)}: Number of pops in MultiPop + }} \section{Functions}{ \itemize{ diff --git a/man/Pop-class.Rd b/man/Pop-class.Rd index d2dd0004..0cfa3065 100644 --- a/man/Pop-class.Rd +++ b/man/Pop-class.Rd @@ -6,6 +6,7 @@ \alias{[,Pop-method} \alias{c,Pop-method} \alias{show,Pop-method} +\alias{length,Pop-method} \title{Population} \usage{ \S4method{[}{Pop}(x, i) @@ -13,6 +14,8 @@ \S4method{c}{Pop}(x, ...) \S4method{show}{Pop}(object) + +\S4method{length}{Pop}(x) } \arguments{ \item{x}{a 'Pop' object} @@ -35,6 +38,8 @@ phenotypes, and pedigrees. \item \code{show(Pop)}: Show population summary +\item \code{length(Pop)}: Number of individuals in Pop (the same as nInd()) + }} \section{Slots}{ @@ -66,13 +71,20 @@ are nInd rows and a variable number of columns.} \item{\code{fixEff}}{a fixed effect relating to the phenotype. Used by genomic selection models but otherwise ignored.} -\item{\code{misc}}{a list whose elements correspond to individuals in the -population. This list is normally empty and exists solely as an +\item{\code{misc}}{a list whose elements correspond to additional miscellaneous +nodes with the items for individuals in the population (see example in +\code{\link{newPop}}). +This list is normally empty and exists solely as an open slot available for uses to store extra information about individuals.} -\item{\code{miscPop}}{a list of any length containing optional meta data for the -population. This list is empty unless information is supplied by the user. -Note that the list is emptied every time the population is subsetted.} +\item{\code{miscPop}}{a list of any length containing optional meta data for the +population (see example in \code{\link{newPop}}). +This list is empty unless information is supplied by the user. +Note that the list is emptied every time the population is subsetted or +combined because the meta data for old population might not be valid anymore.} }} +\seealso{ +\code{\link{newPop}}, \code{\link{newEmptyPop}}, \code{\link{resetPop}} +} diff --git a/man/RRBLUP.Rd b/man/RRBLUP.Rd index 17cff96a..02ff747a 100644 --- a/man/RRBLUP.Rd +++ b/man/RRBLUP.Rd @@ -18,25 +18,25 @@ RRBLUP( \arguments{ \item{pop}{a \code{\link{Pop-class}} to serve as the training population} -\item{traits}{an integer indicating the trait or traits to model, a vector of trait names, +\item{traits}{an integer indicating the trait or traits to model, a vector of trait names, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} -\item{maxIter}{maximum number of iterations. Only used +\item{maxIter}{maximum number of iterations. Only used when number of traits is greater than 1.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ diff --git a/man/RRBLUP2.Rd b/man/RRBLUP2.Rd index 12b548ce..813d3bcf 100644 --- a/man/RRBLUP2.Rd +++ b/man/RRBLUP2.Rd @@ -23,69 +23,69 @@ RRBLUP2( \item{pop}{a \code{\link{Pop-class}} to serve as the training population} \item{traits}{an integer indicating the trait to model, a trait name, or a -function of the traits returning a single value. Unlike \code{\link{RRBLUP}}, +function of the traits returning a single value. Unlike \code{\link{RRBLUP}}, only univariate models are supported.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations.} -\item{Vu}{marker effect variance. If value is NULL, a +\item{Vu}{marker effect variance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{Ve}{error variance. If value is NULL, a +\item{Ve}{error variance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{useEM}{use EM to solve variance components. If false, +\item{useEM}{use EM to solve variance components. If false, the initial values are considered true.} \item{tol}{tolerance for EM algorithm convergence} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ -Fits an RR-BLUP model for genomic predictions. This implementation is -meant for situations where \code{\link{RRBLUP}} is too slow. Note that -RRBLUP2 is only faster in certain situations, see details below. Most +Fits an RR-BLUP model for genomic predictions. This implementation is +meant for situations where \code{\link{RRBLUP}} is too slow. Note that +RRBLUP2 is only faster in certain situations, see details below. Most users should use \code{\link{RRBLUP}}. } \details{ -The RRBLUP2 function works best when the number of markers is not -too large. This is because it solves the RR-BLUP problem by setting -up and solving Henderson's mixed model equations. Solving these equations -involves a square matrix with dimensions equal to the number of fixed -effects plus the number of random effects (markers). Whereas the \code{\link{RRBLUP}} -function solves the RR-BLUP problem using the EMMA approach. This approach involves -a square matrix with dimensions equal to the number of phenotypic records. This means -that the RRBLUP2 function uses less memory than RRBLUP when the number of markers -is approximately equal to or smaller than the number of phenotypic records. - -The RRBLUP2 function is not recommend for cases where the variance components are -unknown. This is uses the EM algorithm to solve for unknown variance components, -which is generally considerably slower than the EMMA approach of \code{\link{RRBLUP}}. -The number of iterations for the EM algorithm is set by maxIter. The default value -is typically too small for convergence. When the algorithm fails to converge a -warning is displayed, but results are given for the last iteration. These results may -be "good enough". However we make no claim to this effect, because we can not generalize +The RRBLUP2 function works best when the number of markers is not +too large. This is because it solves the RR-BLUP problem by setting +up and solving Henderson's mixed model equations. Solving these equations +involves a square matrix with dimensions equal to the number of fixed +effects plus the number of random effects (markers). Whereas the \code{\link{RRBLUP}} +function solves the RR-BLUP problem using the EMMA approach. This approach involves +a square matrix with dimensions equal to the number of phenotypic records. This means +that the RRBLUP2 function uses less memory than RRBLUP when the number of markers +is approximately equal to or smaller than the number of phenotypic records. + +The RRBLUP2 function is not recommend for cases where the variance components are +unknown. This is uses the EM algorithm to solve for unknown variance components, +which is generally considerably slower than the EMMA approach of \code{\link{RRBLUP}}. +The number of iterations for the EM algorithm is set by maxIter. The default value +is typically too small for convergence. When the algorithm fails to converge a +warning is displayed, but results are given for the last iteration. These results may +be "good enough". However we make no claim to this effect, because we can not generalize to all possible use cases. -The RRBLUP2 function can quickly solve the mixed model equations without estimating variance -components. The variance components are set by defining Vu and Ve. Estimation of components -is suppressed by setting useEM to false. This may be useful if the model is being retrained -multiple times during the simulation. You could run \code{\link{RRBLUP}} function the first -time the model is trained, and then use the variance components from this output for all -future runs with the RRBLUP2 functions. Again, we can make no claim to the general robustness +The RRBLUP2 function can quickly solve the mixed model equations without estimating variance +components. The variance components are set by defining Vu and Ve. Estimation of components +is suppressed by setting useEM to false. This may be useful if the model is being retrained +multiple times during the simulation. You could run \code{\link{RRBLUP}} function the first +time the model is trained, and then use the variance components from this output for all +future runs with the RRBLUP2 functions. Again, we can make no claim to the general robustness of this approach. } \examples{ diff --git a/man/RRBLUPMemUse.Rd b/man/RRBLUPMemUse.Rd index 484820c5..d57f080a 100644 --- a/man/RRBLUPMemUse.Rd +++ b/man/RRBLUPMemUse.Rd @@ -11,7 +11,7 @@ RRBLUPMemUse(nInd, nMarker, model = "REG") \item{nMarker}{the number of markers per individual} -\item{model}{either "REG", "GCA", or "SCA" for \code{\link{RRBLUP}} +\item{model}{either "REG", "GCA", or "SCA" for \code{\link{RRBLUP}} \code{\link{RRBLUP_GCA}} and \code{\link{RRBLUP_SCA}} respectively.} } \value{ @@ -19,7 +19,7 @@ Returns an estimate for the required gigabytes of RAM } \description{ Estimates the amount of RAM needed to run the \code{\link{RRBLUP}} -and its related functions for a given training population size. +and its related functions for a given training population size. Note that this function may underestimate total usage. } \examples{ diff --git a/man/RRBLUP_D.Rd b/man/RRBLUP_D.Rd index 0c8e7be1..e12f991b 100644 --- a/man/RRBLUP_D.Rd +++ b/man/RRBLUP_D.Rd @@ -21,26 +21,26 @@ RRBLUP_D( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} -\item{maxIter}{maximum number of iterations. Only used +\item{maxIter}{maximum number of iterations. Only used when number of traits is greater than 1.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ -Fits an RR-BLUP model for genomic predictions that includes +Fits an RR-BLUP model for genomic predictions that includes dominance effects. } \examples{ diff --git a/man/RRBLUP_D2.Rd b/man/RRBLUP_D2.Rd index 15428e1c..84747fd3 100644 --- a/man/RRBLUP_D2.Rd +++ b/man/RRBLUP_D2.Rd @@ -26,43 +26,43 @@ RRBLUP_D2( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} -\item{maxIter}{maximum number of iterations. Only used +\item{maxIter}{maximum number of iterations. Only used when number of traits is greater than 1.} -\item{Va}{marker effect variance for additive effects. If value is NULL, +\item{Va}{marker effect variance for additive effects. If value is NULL, a reasonable starting point is chosen automatically.} -\item{Vd}{marker effect variance for dominance effects. If value is NULL, +\item{Vd}{marker effect variance for dominance effects. If value is NULL, a reasonable starting point is chosen automatically.} -\item{Ve}{error variance. If value is NULL, a +\item{Ve}{error variance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{useEM}{use EM to solve variance components. If false, +\item{useEM}{use EM to solve variance components. If false, the initial values are considered true.} \item{tol}{tolerance for EM algorithm convergence} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ -Fits an RR-BLUP model for genomic predictions that includes -dominance effects. This implementation is meant for situations where -\code{\link{RRBLUP_D}} is too slow. Note that RRBLUP_D2 -is only faster in certain situations. Most users should use +Fits an RR-BLUP model for genomic predictions that includes +dominance effects. This implementation is meant for situations where +\code{\link{RRBLUP_D}} is too slow. Note that RRBLUP_D2 +is only faster in certain situations. Most users should use \code{\link{RRBLUP_D}}. } \examples{ diff --git a/man/RRBLUP_GCA.Rd b/man/RRBLUP_GCA.Rd index bf7b2e53..dee4e4d1 100644 --- a/man/RRBLUP_GCA.Rd +++ b/man/RRBLUP_GCA.Rd @@ -21,27 +21,27 @@ RRBLUP_GCA( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations for convergence.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ Fits an RR-BLUP model that estimates seperate marker effects for females and males. Useful for predicting GCA of parents -in single cross hybrids. Can also predict performance of specific +in single cross hybrids. Can also predict performance of specific single cross hybrids. } \examples{ diff --git a/man/RRBLUP_GCA2.Rd b/man/RRBLUP_GCA2.Rd index d25928c7..fb684d72 100644 --- a/man/RRBLUP_GCA2.Rd +++ b/man/RRBLUP_GCA2.Rd @@ -26,42 +26,42 @@ RRBLUP_GCA2( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations for convergence.} -\item{VuF}{marker effect variance for females. If value is NULL, a +\item{VuF}{marker effect variance for females. If value is NULL, a reasonable starting point is chosen automatically.} -\item{VuM}{marker effect variance for males. If value is NULL, a +\item{VuM}{marker effect variance for males. If value is NULL, a reasonable starting point is chosen automatically.} -\item{Ve}{error variance. If value is NULL, a +\item{Ve}{error variance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{useEM}{use EM to solve variance components. If false, +\item{useEM}{use EM to solve variance components. If false, the initial values are considered true.} \item{tol}{tolerance for EM algorithm convergence} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ Fits an RR-BLUP model that estimates seperate marker effects for -females and males. This implementation is meant for situations where -\code{\link{RRBLUP_GCA}} is too slow. Note that RRBLUP_GCA2 -is only faster in certain situations. Most users should use +females and males. This implementation is meant for situations where +\code{\link{RRBLUP_GCA}} is too slow. Note that RRBLUP_GCA2 +is only faster in certain situations. Most users should use \code{\link{RRBLUP_GCA}}. } \examples{ diff --git a/man/RRBLUP_SCA.Rd b/man/RRBLUP_SCA.Rd index 12201d48..f7b0398d 100644 --- a/man/RRBLUP_SCA.Rd +++ b/man/RRBLUP_SCA.Rd @@ -21,26 +21,26 @@ RRBLUP_SCA( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations for convergence.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ -An extention of \code{\link{RRBLUP_GCA}} that adds dominance effects. -Note that we have not seen any consistent benefit of this model over +An extention of \code{\link{RRBLUP_GCA}} that adds dominance effects. +Note that we have not seen any consistent benefit of this model over \code{\link{RRBLUP_GCA}}. } \examples{ diff --git a/man/RRBLUP_SCA2.Rd b/man/RRBLUP_SCA2.Rd index e4985642..86b023d3 100644 --- a/man/RRBLUP_SCA2.Rd +++ b/man/RRBLUP_SCA2.Rd @@ -27,45 +27,45 @@ RRBLUP_SCA2( \item{traits}{an integer indicating the trait to model, a trait name, or a function of the traits returning a single value.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations for convergence.} -\item{VuF}{marker effect variance for females. If value is NULL, a +\item{VuF}{marker effect variance for females. If value is NULL, a reasonable starting point is chosen automatically.} -\item{VuM}{marker effect variance for males. If value is NULL, a +\item{VuM}{marker effect variance for males. If value is NULL, a reasonable starting point is chosen automatically.} -\item{VuD}{marker effect variance for dominance. If value is NULL, a +\item{VuD}{marker effect variance for dominance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{Ve}{error variance. If value is NULL, a +\item{Ve}{error variance. If value is NULL, a reasonable starting point is chosen automatically.} -\item{useEM}{use EM to solve variance components. If false, +\item{useEM}{use EM to solve variance components. If false, the initial values are considered true.} \item{tol}{tolerance for EM algorithm convergence} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ Fits an RR-BLUP model that estimates seperate additive effects for -females and males and a dominance effect. This implementation is meant -for situations where \code{\link{RRBLUP_SCA}} is too slow. Note that -RRBLUP_SCA2 is only faster in certain situations. Most users should use +females and males and a dominance effect. This implementation is meant +for situations where \code{\link{RRBLUP_SCA}} is too slow. Note that +RRBLUP_SCA2 is only faster in certain situations. Most users should use \code{\link{RRBLUP_SCA}}. } \examples{ diff --git a/man/addSegSite.Rd b/man/addSegSite.Rd index 41fb2eb3..61e047a8 100644 --- a/man/addSegSite.Rd +++ b/man/addSegSite.Rd @@ -21,8 +21,8 @@ addSegSite(mapPop, siteName, chr, mapPos, haplo) an object of \code{\link{MapPop-class}} } \description{ -This function allows for adding a new -segregating site with user supplied genotypes to a MapPop. +This function allows for adding a new +segregating site with user supplied genotypes to a MapPop. The position of the site is set using a genetic map position. } \examples{ diff --git a/man/attrition.Rd b/man/attrition.Rd index fdda086d..d39e6427 100644 --- a/man/attrition.Rd +++ b/man/attrition.Rd @@ -16,8 +16,8 @@ be lost to attrition.} an object of \code{\link{Pop-class}} } \description{ -Samples individuals at random to remove from the population. -The user supplies a probability for the individuals to be +Samples individuals at random to remove from the population. +The user supplies a probability for the individuals to be removed from the population. } \examples{ diff --git a/man/editGenome.Rd b/man/editGenome.Rd index cb26bc90..f7f924b9 100644 --- a/man/editGenome.Rd +++ b/man/editGenome.Rd @@ -11,10 +11,10 @@ editGenome(pop, ind, chr, segSites, allele, simParam = NULL) \item{ind}{a vector of individuals to edit} -\item{chr}{a vector of chromosomes to edit. +\item{chr}{a vector of chromosomes to edit. Length must match length of segSites.} -\item{segSites}{a vector of segregating sites to edit. Length must +\item{segSites}{a vector of segregating sites to edit. Length must match length of chr.} \item{allele}{either 0 or 1 for desired allele} @@ -25,8 +25,8 @@ match length of chr.} Returns an object of \code{\link{Pop-class}} } \description{ -Edits selected loci of selected individuals to a homozygous -state for either the 1 or 0 allele. The gv slot is recalculated to +Edits selected loci of selected individuals to a homozygous +state for either the 1 or 0 allele. The gv slot is recalculated to reflect the any changes due to editing, but other slots remain the same. } \examples{ @@ -41,9 +41,9 @@ SP$addTraitA(10) #Create population pop = newPop(founderPop, simParam=SP) -#Change individual 1 to homozygous for the 1 allele +#Change individual 1 to homozygous for the 1 allele #at locus 1, chromosome 1 -pop2 = editGenome(pop, ind=1, chr=1, segSites=1, +pop2 = editGenome(pop, ind=1, chr=1, segSites=1, allele=1, simParam=SP) } diff --git a/man/editGenomeTopQtl.Rd b/man/editGenomeTopQtl.Rd index 8a2d8b35..667b6923 100644 --- a/man/editGenomeTopQtl.Rd +++ b/man/editGenomeTopQtl.Rd @@ -23,7 +23,7 @@ editGenomeTopQtl(pop, ind, nQtl, trait = 1, increase = TRUE, simParam = NULL) Returns an object of \code{\link{Pop-class}} } \description{ -Edits the top QTL (with the largest additive effect) to a homozygous +Edits the top QTL (with the largest additive effect) to a homozygous state for the allele increasing. Only nonfixed QTL are edited The gv slot is recalculated to reflect the any changes due to editing, but other slots remain the same. } @@ -39,7 +39,7 @@ SP$addTraitA(10) #Create population pop = newPop(founderPop, simParam=SP) -#Change up to 10 loci for individual 1 +#Change up to 10 loci for individual 1 pop2 = editGenomeTopQtl(pop, ind=1, nQtl=10, simParam=SP) - + } diff --git a/man/fastRRBLUP.Rd b/man/fastRRBLUP.Rd index 6bf7483b..8c45eb65 100644 --- a/man/fastRRBLUP.Rd +++ b/man/fastRRBLUP.Rd @@ -20,38 +20,38 @@ fastRRBLUP( \arguments{ \item{pop}{a \code{\link{Pop-class}} to serve as the training population} -\item{traits}{an integer indicating the trait to model, a trait name, -or a function of the traits returning a single value. Only univariate models +\item{traits}{an integer indicating the trait to model, a trait name, +or a function of the traits returning a single value. Only univariate models are supported.} -\item{use}{train model using phenotypes "pheno", genetic values "gv", +\item{use}{train model using phenotypes "pheno", genetic values "gv", estimated breeding values "ebv", breeding values "bv", or randomly "rand"} -\item{snpChip}{an integer indicating which SNP chip genotype +\item{snpChip}{an integer indicating which SNP chip genotype to use} -\item{useQtl}{should QTL genotypes be used instead of a SNP chip. -If TRUE, snpChip specifies which trait's QTL to use, and thus these +\item{useQtl}{should QTL genotypes be used instead of a SNP chip. +If TRUE, snpChip specifies which trait's QTL to use, and thus these QTL may not match the QTL underlying the phenotype supplied in traits.} \item{maxIter}{maximum number of iterations.} -\item{Vu}{marker effect variance. If value is NULL, a +\item{Vu}{marker effect variance. If value is NULL, a reasonable value is chosen automatically.} -\item{Ve}{error variance. If value is NULL, a +\item{Ve}{error variance. If value is NULL, a reasonable value is chosen automatically.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for traits} } \description{ -Solves an RR-BLUP model for genomic predictions given known variance -components. This implementation is meant as a fast and low memory -alternative to \code{\link{RRBLUP}} or \code{\link{RRBLUP2}}. Unlike -the those functions, the fastRRBLUP does not fit fixed effects (other +Solves an RR-BLUP model for genomic predictions given known variance +components. This implementation is meant as a fast and low memory +alternative to \code{\link{RRBLUP}} or \code{\link{RRBLUP2}}. Unlike +the those functions, the fastRRBLUP does not fit fixed effects (other than the intercept) or account for unequal replication. } \examples{ diff --git a/man/genParam.Rd b/man/genParam.Rd index 8ad6ba82..d4d5fc8e 100644 --- a/man/genParam.Rd +++ b/man/genParam.Rd @@ -42,6 +42,8 @@ genParam(pop, simParam = NULL) \item{gv_a}{a matrix of additive genetic values with dimensions nInd by nTraits} \item{gv_d}{a matrix of dominance genetic values with dimensions nInd by nTraits} \item{gv_aa}{a matrix of additive-by-additive genetic values with dimensions nInd by nTraits} +\item{alpha}{a list of average allele subsitution effects with length nTraits} +\item{alpha_HW}{a list of average allele subsitution effects at Hardy-Weinberg equilibrium with length nTraits} } } \description{ diff --git a/man/getMisc.Rd b/man/getMisc.Rd deleted file mode 100644 index 6c2babcb..00000000 --- a/man/getMisc.Rd +++ /dev/null @@ -1,57 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R -\name{getMisc} -\alias{getMisc} -\title{Get miscelaneous information in a population} -\usage{ -getMisc(x, node = NULL) -} -\arguments{ -\item{x}{\code{\link{Pop-class}}} - -\item{node}{character, name of the node to get from the \code{x@misc} slot; -if \code{NULL} the whole \code{x@misc} slot is returned} -} -\value{ -The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]} -} -\description{ -Get miscelaneous information in a population -} -\examples{ -founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 100) -SP <- SimParam$new(founderGenomes) -\dontshow{SP$nThreads = 1L} -basePop <- newPop(founderGenomes) - -basePop <- setMisc(basePop, node = "info", value = 1) -basePop@misc -getMisc(x = basePop, node = "info") - -basePop <- setMisc(basePop, node = "info2", value = c("A", "B", "C")) -basePop@misc -getMisc(x = basePop, node = "info2") - -n <- nInd(basePop) -location <- vector(mode = "list", length = n) -for (ind in seq_len(n)) { - location[[ind]] <- runif(n = 2, min = 0, max = 100) -} -location -basePop <- setMisc(basePop, node = "location", value = location) -basePop@misc -getMisc(x = basePop, node = "location") - -n <- nInd(basePop) -location <- vector(mode = "list", length = n) -for (ind in c(1, 3)) { - location[[ind]] <- runif(n = 2, min = 0, max = 100) -} -location -basePop <- setMisc(basePop, node = "location", value = location) -basePop@misc -getMisc(x = basePop, node = "location") - -getMisc(x = basePop) - -} diff --git a/man/getPed.Rd b/man/getPed.Rd index 86475e4c..96dae5e8 100644 --- a/man/getPed.Rd +++ b/man/getPed.Rd @@ -10,8 +10,8 @@ getPed(pop) \item{pop}{a population} } \description{ -Returns the population's pedigree as stored in the -id, mother and father slots. NULL is returned if the +Returns the population's pedigree as stored in the +id, mother and father slots. NULL is returned if the input population lacks the required. } \examples{ diff --git a/man/hybridCross.Rd b/man/hybridCross.Rd index 66c594a2..51524931 100644 --- a/man/hybridCross.Rd +++ b/man/hybridCross.Rd @@ -17,7 +17,7 @@ hybridCross( \item{males}{male population, an object of \code{\link{Pop-class}}} -\item{crossPlan}{either "testcross" for all possible combinantions +\item{crossPlan}{either "testcross" for all possible combinations or a matrix with two columns for designed crosses} \item{returnHybridPop}{should results be returned as @@ -27,7 +27,7 @@ or a matrix with two columns for designed crosses} \item{simParam}{an object of \code{\link{SimParam}}} } \description{ -A convience function for hybrid plant breeding simulations. Allows for +A convenient function for hybrid plant breeding simulations. Allows for easy specification of a test cross scheme and/or creation of an object of \code{\link{HybridPop-class}}. Note that the \code{\link{HybridPop-class}} should only be used if the parents were created using the \code{\link{makeDH}} diff --git a/man/meanEBV.Rd b/man/meanEBV.Rd new file mode 100644 index 00000000..a8808989 --- /dev/null +++ b/man/meanEBV.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/popSummary.R +\name{meanEBV} +\alias{meanEBV} +\title{Mean estimated breeding values} +\usage{ +meanEBV(pop) +} +\arguments{ +\item{pop}{an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}}} +} +\description{ +Returns the mean estimated breeding values for all traits +} +\examples{ +#Create founder haplotypes +founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) + +#Set simulation parameters +SP = SimParam$new(founderPop) +SP$addTraitA(10) +trtH2 = 0.5 +SP$setVarE(h2=trtH2) +\dontshow{SP$nThreads = 1L} + +#Create population +pop = newPop(founderPop, simParam=SP) +pop@ebv = trtH2 * (pop@pheno - meanP(pop)) #ind performance based EBV +meanEBV(pop) + +} diff --git a/man/mergePops.Rd b/man/mergePops.Rd index e77eb8ec..4a466ef5 100644 --- a/man/mergePops.Rd +++ b/man/mergePops.Rd @@ -27,6 +27,9 @@ SP = SimParam$new(founderPop) #Create a list of populations and merge list pop = newPop(founderPop, simParam=SP) +pop@misc$tmp = rnorm(n=10) +pop@misc$tmp2 = rnorm(n=10) + popList = list(pop, pop) pop2 = mergePops(popList) diff --git a/man/mutate.Rd b/man/mutate.Rd index 5545dd32..9430109d 100644 --- a/man/mutate.Rd +++ b/man/mutate.Rd @@ -16,17 +16,17 @@ mutate(pop, mutRate = 2.5e-08, returnPos = FALSE, simParam = NULL) \item{simParam}{an object of \code{\link{SimParam}}} } \value{ -an object of \code{\link{Pop-class}} if -returnPos=FALSE or a list containing a -\code{\link{Pop-class}} and a data.frame containing the +an object of \code{\link{Pop-class}} if +returnPos=FALSE or a list containing a +\code{\link{Pop-class}} and a data.frame containing the postions of mutations if returnPos=TRUE } \description{ -Adds random mutations to individuals in a -population. Note that any existing phenotypes -or EBVs are kept. Thus, the user will need to run -\code{\link{setPheno}} and/or \code{\link{setEBV}} -to generate new phenotypes or EBVs that reflect +Adds random mutations to individuals in a +population. Note that any existing phenotypes +or EBVs are kept. Thus, the user will need to run +\code{\link{setPheno}} and/or \code{\link{setEBV}} +to generate new phenotypes or EBVs that reflect changes introduced by the new mutations. } \examples{ diff --git a/man/newMapPop.Rd b/man/newMapPop.Rd index 9d7cb1d5..ce391ec6 100644 --- a/man/newMapPop.Rd +++ b/man/newMapPop.Rd @@ -9,7 +9,7 @@ newMapPop(genMap, haplotypes, inbred = FALSE, ploidy = 2L) \arguments{ \item{genMap}{a list of genetic maps} -\item{haplotypes}{a list of matrices or data.frames that +\item{haplotypes}{a list of matrices or data.frames that can be coerced to matrices. See details.} \item{inbred}{are individuals fully inbred} @@ -20,19 +20,19 @@ can be coerced to matrices. See details.} an object of \code{\link{MapPop-class}} } \description{ -Creates a new \code{\link{MapPop-class}} from user supplied +Creates a new \code{\link{MapPop-class}} from user supplied genetic maps and haplotypes. } \details{ -Each item of genMap must be a vector of ordered genetic lengths in -Morgans. The first value must be zero. The length of the vector +Each item of genMap must be a vector of ordered genetic lengths in +Morgans. The first value must be zero. The length of the vector determines the number of segregating sites on the chromosome. -Each item of haplotypes must be coercible to a matrix. The columns -of this matrix correspond to segregating sites. The number of rows -must match the number of individuals times the ploidy if using -inbred=FALSE. If using inbred=TRUE, the number of rows must equal -the number of individuals. The haplotypes can be stored as numeric, +Each item of haplotypes must be coercible to a matrix. The columns +of this matrix correspond to segregating sites. The number of rows +must match the number of individuals times the ploidy if using +inbred=FALSE. If using inbred=TRUE, the number of rows must equal +the number of individuals. The haplotypes can be stored as numeric, integer or raw. The underlying C++ function will use raw. } \examples{ @@ -40,7 +40,7 @@ integer or raw. The underlying C++ function will use raw. # Each chromosome contains 11 equally spaced segregating sites genMap = list(seq(0,1,length.out=11), seq(0,1,length.out=11)) - + # Create haplotypes for 10 outbred individuals chr1 = sample(x=0:1,size=20*11,replace=TRUE) chr1 = matrix(chr1,nrow=20,ncol=11) diff --git a/man/newPop.Rd b/man/newPop.Rd index 9cf11733..9844c23c 100644 --- a/man/newPop.Rd +++ b/man/newPop.Rd @@ -20,10 +20,19 @@ Returns an object of \code{\link{Pop-class}} \description{ Creates an initial \code{\link{Pop-class}} from an object of \code{\link{MapPop-class}} or \code{\link{NamedMapPop-class}}. -The function is intended for us with output from functions such +The function is intended for use with output from functions such as \code{\link{runMacs}}, \code{\link{newMapPop}}, or \code{\link{quickHaplo}}. } +\details{ +Note that \code{newPop} takes genomes from the + \code{rawPop} and uses them without recombination! Hence, if you + call \code{newPop(rawPop = founderGenomes)} twice, you will get + two sets of individuals with different id but the same genomes. + To get genetically different sets of individuals you can subset the + \code{rawPop} input, say first half for one set and the second half + for the other set. +} \examples{ #Create founder haplotypes founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) @@ -36,4 +45,11 @@ SP$addTraitA(10) pop = newPop(founderPop, simParam=SP) isPop(pop) +#Misc +pop@misc$tmp1 = rnorm(n=2) +pop@misc$tmp2 = rnorm(n=2) + +#MiscPop +pop@miscPop$tmp1 = sum(pop@misc$tmp1) +pop@miscPop$tmp2 = sum(pop@misc$tmp2) } diff --git a/man/pedigreeCross.Rd b/man/pedigreeCross.Rd index 9eb3d5d2..8dcf38c7 100644 --- a/man/pedigreeCross.Rd +++ b/man/pedigreeCross.Rd @@ -76,5 +76,4 @@ mother = c(0,0,1,3:9) father = c(0,0,2,3:9) pop2 = pedigreeCross(pop, id, mother, father, simParam=SP) - } diff --git a/man/quickHaplo.Rd b/man/quickHaplo.Rd index 55a4a9dc..8b73ab16 100644 --- a/man/quickHaplo.Rd +++ b/man/quickHaplo.Rd @@ -23,8 +23,8 @@ quickHaplo(nInd, nChr, segSites, genLen = 1, ploidy = 2L, inbred = FALSE) an object of \code{\link{MapPop-class}} } \description{ -Rapidly simulates founder haplotypes by randomly -sampling 0s and 1s. This is equivalent to having all loci with +Rapidly simulates founder haplotypes by randomly +sampling 0s and 1s. This is equivalent to having all loci with allele frequency 0.5 and being in linkage equilibrium. } \examples{ diff --git a/man/runMacs.Rd b/man/runMacs.Rd index 4cc54bff..83ba9f9c 100644 --- a/man/runMacs.Rd +++ b/man/runMacs.Rd @@ -22,7 +22,7 @@ runMacs( \item{nChr}{number of chromosomes to simulate} -\item{segSites}{number of segregating sites to keep per chromosome. A +\item{segSites}{number of segregating sites to keep per chromosome. A value of NULL results in all sites being retained.} \item{inbred}{should founder individuals be inbred} @@ -35,37 +35,37 @@ value of NULL results in all sites being retained.} \item{manualCommand}{user provided MaCS options. For advanced users only.} -\item{manualGenLen}{user provided genetic length. This must be supplied if using -manualCommand. If not using manualCommand, this value will replace the predefined -genetic length for the species. However, this the genetic length is only used by -AlphaSimR and is not passed to MaCS, so MaCS still uses the predefined genetic length. +\item{manualGenLen}{user provided genetic length. This must be supplied if using +manualCommand. If not using manualCommand, this value will replace the predefined +genetic length for the species. However, this the genetic length is only used by +AlphaSimR and is not passed to MaCS, so MaCS still uses the predefined genetic length. For advanced users only.} -\item{nThreads}{if OpenMP is available, this will allow for simulating chromosomes in parallel. +\item{nThreads}{if OpenMP is available, this will allow for simulating chromosomes in parallel. If the value is NULL, the number of threads is automatically detected.} } \value{ an object of \code{\link{MapPop-class}} } \description{ -Uses the MaCS software to produce founder haplotypes +Uses the MaCS software to produce founder haplotypes \insertCite{MaCS}{AlphaSimR}. } \details{ There are currently three species histories available: GENERIC, CATTLE, WHEAT, and MAIZE. -The GENERIC history is meant to be a reasonable all-purpose choice. It runs quickly and -models a population with an effective populations size that has gone through several historic -bottlenecks. This species history is used as the default arguments in the \code{\link{runMacs2}} +The GENERIC history is meant to be a reasonable all-purpose choice. It runs quickly and +models a population with an effective populations size that has gone through several historic +bottlenecks. This species history is used as the default arguments in the \code{\link{runMacs2}} function, so the user should examine this function for the details of how the species is modeled. The CATTLE history is based off of real genome sequence data \insertCite{cattle}{AlphaSimR}. -The WHEAT \insertCite{gaynor_2017}{AlphaSimR} and MAIZE \insertCite{hickey_2014}{AlphaSimR} -histories have been included due to their use in previous simulations. However, it should -be noted that neither faithfully simulates its respective species. This is apparent by -the low number of segregating sites simulated by each history relative to their real-world -analogs. Adjusting these histories to better represent their real-world analogs would result +The WHEAT \insertCite{gaynor_2017}{AlphaSimR} and MAIZE \insertCite{hickey_2014}{AlphaSimR} +histories have been included due to their use in previous simulations. However, it should +be noted that neither faithfully simulates its respective species. This is apparent by +the low number of segregating sites simulated by each history relative to their real-world +analogs. Adjusting these histories to better represent their real-world analogs would result in a drastic increase to runtime. } \examples{ diff --git a/man/runMacs2.Rd b/man/runMacs2.Rd index cb0007f0..f0536cfe 100644 --- a/man/runMacs2.Rd +++ b/man/runMacs2.Rd @@ -36,10 +36,10 @@ runMacs2( \item{mutRate}{per base pair mutation rate} -\item{histNe}{effective population size in previous +\item{histNe}{effective population size in previous generations} -\item{histGen}{number of generations ago for effective +\item{histGen}{number of generations ago for effective population sizes given in histNe} \item{inbred}{should founder individuals be inbred} @@ -48,25 +48,25 @@ population sizes given in histNe} \item{ploidy}{ploidy level of organism} -\item{returnCommand}{should the command passed to manualCommand in -\code{\link{runMacs}} be returned. If TRUE, MaCS will not be called and +\item{returnCommand}{should the command passed to manualCommand in +\code{\link{runMacs}} be returned. If TRUE, MaCS will not be called and the command is returned instead.} -\item{nThreads}{if OpenMP is available, this will allow for simulating chromosomes in parallel. +\item{nThreads}{if OpenMP is available, this will allow for simulating chromosomes in parallel. If the value is NULL, the number of threads is automatically detected.} } \value{ -an object of \code{\link{MapPop-class}} or if -returnCommand is true a string giving the MaCS command passed to +an object of \code{\link{MapPop-class}} or if +returnCommand is true a string giving the MaCS command passed to the manualCommand argument of \code{\link{runMacs}}. } \description{ -A wrapper function for \code{\link{runMacs}}. This wrapper is designed +A wrapper function for \code{\link{runMacs}}. This wrapper is designed to provide a more intuitive interface for writing custom commands -in MaCS \insertCite{MaCS}{AlphaSimR}. It effectively automates the creation -of an appropriate line for the manualCommand argument in \code{\link{runMacs}} -using user supplied variables, but only allows for a subset of the functionality -offered by this argument. The default arguments of this function were chosen to match +in MaCS \insertCite{MaCS}{AlphaSimR}. It effectively automates the creation +of an appropriate line for the manualCommand argument in \code{\link{runMacs}} +using user supplied variables, but only allows for a subset of the functionality +offered by this argument. The default arguments of this function were chosen to match species="GENERIC" in \code{\link{runMacs}}. } \examples{ @@ -75,6 +75,19 @@ species="GENERIC" in \code{\link{runMacs}}. # The command is equivalent to using species="GENERIC" in runMacs \dontrun{ founderPop = runMacs2(nInd=10,nChr=1,segSites=100) + +# runMacs() Implementation of the cattle demography following +# Macleod et al. (2013) https://doi.org/10.1093/molbev/mst125 +cattleChrSum = 2.8e9 # https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_002263795.3/ +(cattleChrBp = cattleChrSum / 30) +recRate = 9.26e-09 +(cattleGenLen = recRate * cattleChrBp) +mutRate = 1.20e-08 +runMacs2(nInd = 10, nChr = 1, Ne = 90, bp = cattleChrBp, + genLen = cattleGenLen, mutRate = 1.20e-08, + histNe = c(120, 250, 350, 1000, 1500, 2000, 2500, 3500, 7000, 10000, 17000, 62000), + histGen = c( 3, 6, 12, 18, 24, 154, 454, 654, 1754, 2354, 3354, 33154), + returnCommand = TRUE) } } \references{ diff --git a/man/sampleHaplo.Rd b/man/sampleHaplo.Rd index 89b8eafe..17d6a254 100644 --- a/man/sampleHaplo.Rd +++ b/man/sampleHaplo.Rd @@ -7,14 +7,14 @@ sampleHaplo(mapPop, nInd, inbred = FALSE, ploidy = NULL, replace = TRUE) } \arguments{ -\item{mapPop}{the \code{\link{MapPop-class}} used to +\item{mapPop}{the \code{\link{MapPop-class}} used to sample haplotypes} \item{nInd}{the number of individuals to create} \item{inbred}{should new individuals be fully inbred} -\item{ploidy}{new ploidy level for organism. If NULL, +\item{ploidy}{new ploidy level for organism. If NULL, the ploidy level of the mapPop is used.} \item{replace}{should haplotypes be sampled with replacement} @@ -23,7 +23,7 @@ the ploidy level of the mapPop is used.} an object of \code{\link{MapPop-class}} } \description{ -Creates a new \code{\link{MapPop-class}} from an existing +Creates a new \code{\link{MapPop-class}} from an existing \code{\link{MapPop-class}} by randomly sampling haplotypes. } \examples{ diff --git a/man/selIndex.Rd b/man/selIndex.Rd index 801f506b..47769cda 100644 --- a/man/selIndex.Rd +++ b/man/selIndex.Rd @@ -14,9 +14,9 @@ selIndex(Y, b, scale = FALSE) \item{scale}{should Y be scaled and centered} } \description{ -Calculates values of a selection index given trait values and -weights. This function is intended to be used in combination with -selection functions working on populations such as +Calculates values of a selection index given trait values and +weights. This function is intended to be used in combination with +selection functions working on populations such as \code{\link{selectInd}}. } \examples{ @@ -40,7 +40,7 @@ b = smithHazel(econWt, varG(pop), varP(pop)) #Selection 2 best individuals using Smith-Hazel index #selIndex is used as a trait -pop2 = selectInd(pop, nInd=2, trait=selIndex, +pop2 = selectInd(pop, nInd=2, trait=selIndex, simParam=SP, b=b) } diff --git a/man/selectFam.Rd b/man/selectFam.Rd index e4e47f76..2430b2e9 100644 --- a/man/selectFam.Rd +++ b/man/selectFam.Rd @@ -26,12 +26,17 @@ selectFam( \item{trait}{the trait for selection. Either a number indicating a single trait or a function returning a vector of length nInd. -The function must work on a vector or matrix of \code{use} values. -See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}.} +The function must work on a vector or matrix of \code{use} values as +\code{trait(pop@use, ...)} - depending on what \code{use} is. +See the examples and \code{\link{selIndex}}.} -\item{use}{select on genetic values "gv", estimated -breeding values "ebv", breeding values "bv", phenotypes "pheno", -or randomly "rand"} +\item{use}{the selection criterion. Either a character +(genetic values "gv", estimated breeding values "ebv", breeding values "bv", +phenotypes "pheno", or randomly "rand") or +a function returning a vector of length nInd. +The function must work on \code{pop} as \code{use(pop, trait, ...)} or +as \code{trait(pop@use, ...)} depending on what \code{trait} is. +See the examples.} \item{sex}{which sex to select. Use "B" for both, "F" for females and "M" for males. If the simulation is not using sexes, @@ -53,7 +58,7 @@ individuals is returned.} \item{simParam}{an object of \code{\link{SimParam}}} \item{...}{additional arguments if using a function for -trait} +\code{trait} and \code{use}} } \value{ Returns an object of \code{\link{Pop-class}}, diff --git a/man/selectInd.Rd b/man/selectInd.Rd index 45ccb354..459b7e41 100644 --- a/man/selectInd.Rd +++ b/man/selectInd.Rd @@ -25,12 +25,17 @@ selectInd( \item{trait}{the trait for selection. Either a number indicating a single trait or a function returning a vector of length nInd. -The function must work on a vector or matrix of \code{use} values. +The function must work on a vector or matrix of \code{use} values as +\code{trait(pop@use, ...)} - depending on what \code{use} is. See the examples and \code{\link{selIndex}}.} -\item{use}{select on genetic values "gv", estimated -breeding values "ebv", breeding values "bv", phenotypes "pheno", -or randomly "rand"} +\item{use}{the selection criterion. Either a character +(genetic values "gv", estimated breeding values "ebv", breeding values "bv", +phenotypes "pheno", or randomly "rand") or +a function returning a vector of length nInd. +The function must work on \code{pop} as \code{use(pop, trait, ...)} or +as \code{trait(pop@use, ...)} depending on what \code{trait} is. +See the examples.} \item{sex}{which sex to select. Use "B" for both, "F" for females and "M" for males. If the simulation is not using sexes, @@ -48,7 +53,7 @@ individuals is returned.} \item{simParam}{an object of \code{\link{SimParam}}} \item{...}{additional arguments if using a function for -trait} +\code{trait} or \code{use}} } \value{ Returns an object of \code{\link{Pop-class}}, @@ -73,18 +78,30 @@ pop = newPop(founderPop, simParam=SP) #Select top 5 (directional selection) pop2 = selectInd(pop, 5, simParam=SP) -hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -abline(v = pop2@pheno, col = "red", lwd = 2) +hist(pop@pheno); abline(v=pop@pheno, lwd=2) +abline(v=pop2@pheno, col="red", lwd=2) #Select 5 most deviating from an optima (disruptive selection) -squaredDeviation = function(x, optima = 0) (x - optima)^2 -pop3 = selectInd(pop, 5, simParam=SP, trait = squaredDeviation, selectTop = TRUE) -hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -abline(v = pop3@pheno, col = "red", lwd = 2) +squaredDeviation = function(x, optima=0) (x - optima)^2 +pop3 = selectInd(pop, 5, trait=squaredDeviation, selectTop=TRUE, simParam=SP) +hist(pop@pheno); abline(v=pop@pheno, lwd=2) +abline(v=pop3@pheno, col="red", lwd=2) #Select 5 least deviating from an optima (stabilising selection) -pop4 = selectInd(pop, 5, simParam=SP, trait = squaredDeviation, selectTop = FALSE) -hist(pop@pheno); abline(v = pop@pheno, lwd = 2) -abline(v = pop4@pheno, col = "red", lwd = 2) +pop4 = selectInd(pop, 5, trait=squaredDeviation, selectTop=FALSE, simParam=SP) +hist(pop@pheno); abline(v=pop@pheno, lwd=2) +abline(v=pop4@pheno, col="red", lwd=2) + +#Select 5 individuals based on miscelaneous information with use function +pop@misc = list(smth=rnorm(10), smth2=rnorm(10)) +useFunc = function(pop, trait=NULL) pop@misc$smth + pop@misc$smth2 +pop5 = selectInd(pop, 5, use=useFunc, simParam=SP) +pop5@id + +#... equivalent result with the use & trait function +useFunc2 = function(pop, trait=NULL) cbind(pop@misc$smth, pop@misc$smth2) +trtFunc = function(x) rowSums(x) +pop6 = selectInd(pop, 5, trait=trtFunc, use=useFunc2, simParam=SP) +pop6@id } diff --git a/man/selectOP.Rd b/man/selectOP.Rd index c9aeda93..153ebe29 100644 --- a/man/selectOP.Rd +++ b/man/selectOP.Rd @@ -33,12 +33,17 @@ Value ranges from 0 to 1.} \item{trait}{the trait for selection. Either a number indicating a single trait or a function returning a vector of length nInd. -The function must work on a vector or matrix of \code{use} values. -See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}.} +The function must work on a vector or matrix of \code{use} values as +\code{trait(pop@use, ...)} - depending on what \code{use} is. +See the examples and \code{\link{selIndex}}.} -\item{use}{select on genetic values "gv", estimated -breeding values "ebv", breeding values "bv", phenotypes "pheno", -or randomly "rand"} +\item{use}{the selection criterion. Either a character +(genetic values "gv", estimated breeding values "ebv", breeding values "bv", +phenotypes "pheno", or randomly "rand") or +a function returning a vector of length nInd. +The function must work on \code{pop} as \code{use(pop, trait, ...)} or +as \code{trait(pop@use, ...)} depending on what \code{trait} is. +See the examples.} \item{selectTop}{selects highest values if true. Selects lowest values if false.} @@ -48,7 +53,7 @@ Selects lowest values if false.} \item{simParam}{an object of \code{\link{SimParam}}} \item{...}{additional arguments if using a function for -trait} +\code{trait} and \code{use}} } \value{ Returns an object of \code{\link{Pop-class}} diff --git a/man/selectWithinFam.Rd b/man/selectWithinFam.Rd index 3c21ee09..e02d64ad 100644 --- a/man/selectWithinFam.Rd +++ b/man/selectWithinFam.Rd @@ -26,12 +26,17 @@ selectWithinFam( \item{trait}{the trait for selection. Either a number indicating a single trait or a function returning a vector of length nInd. -The function must work on a vector or matrix of \code{use} values. -See the examples in \code{\link{selectInd}} and \code{\link{selIndex}}.} +The function must work on a vector or matrix of \code{use} values as +\code{trait(pop@use, ...)} - depending on what \code{use} is. +See the examples and \code{\link{selIndex}}.} -\item{use}{select on genetic values "gv", estimated -breeding values "ebv", breeding values "bv", phenotypes "pheno", -or randomly "rand"} +\item{use}{the selection criterion. Either a character +(genetic values "gv", estimated breeding values "ebv", breeding values "bv", +phenotypes "pheno", or randomly "rand") or +a function returning a vector of length nInd. +The function must work on \code{pop} as \code{use(pop, trait, ...)} or +as \code{trait(pop@use, ...)} depending on what \code{trait} is. +See the examples.} \item{sex}{which sex to select. Use "B" for both, "F" for females and "M" for males. If the simulation is not using sexes, @@ -53,7 +58,7 @@ individuals is returned.} \item{simParam}{an object of \code{\link{SimParam}}} \item{...}{additional arguments if using a function for -trait} +\code{trait} and \code{use}} } \value{ Returns an object of \code{\link{Pop-class}}, diff --git a/man/setEBV.Rd b/man/setEBV.Rd index c2afd081..ae8a8462 100644 --- a/man/setEBV.Rd +++ b/man/setEBV.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GS.R \name{setEBV} \alias{setEBV} -\title{Set EBV} +\title{Set estimated breeding values (EBV)} \usage{ setEBV( pop, @@ -18,17 +18,17 @@ setEBV( \item{solution}{an object of \code{\link{RRsol-class}}} -\item{value}{the genomic value to be estimated. Can be +\item{value}{the genomic value to be estimated. Can be either "gv", "bv", "female", or "male".} -\item{targetPop}{an optional target population that can -be used when value is "bv", "female", or "male". When -supplied, the allele frequency in the targetPop is used +\item{targetPop}{an optional target population that can +be used when value is "bv", "female", or "male". When +supplied, the allele frequency in the targetPop is used to set these values.} -\item{append}{should estimated values be appended to -existing data in the EBV slot. If TRUE, a new column is -added. If FALSE, existing data is replaced with the +\item{append}{should estimated values be appended to +existing data in the EBV slot. If TRUE, a new column is +added. If FALSE, existing data is replaced with the new estimates.} \item{simParam}{an object of \code{\link{SimParam}}} @@ -37,10 +37,10 @@ new estimates.} Returns an object of \code{\link{Pop-class}} } \description{ -Adds genomic estimated values to a populations's EBV -slot using output from a genomic selection functions. -The genomic estimated values can be either estimated -breeding values, estimated genetic values, or +Adds genomic estimated values to a populations's EBV +slot using output from a genomic selection functions. +The genomic estimated values can be either estimated +breeding values, estimated genetic values, or estimated general combining values. } \examples{ diff --git a/man/setMisc.Rd b/man/setMisc.Rd deleted file mode 100644 index 32801ec8..00000000 --- a/man/setMisc.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/misc.R -\name{setMisc} -\alias{setMisc} -\title{Set miscelaneous information in a population} -\usage{ -setMisc(x, node = NULL, value = NULL) -} -\arguments{ -\item{x}{\code{\link{Pop-class}}} - -\item{node}{character, name of the node to set within the \code{x@misc} slot} - -\item{value, }{value to be saved into \code{x@misc[[*]][[node]]}; length of -\code{value} should be equal to \code{nInd(x)}; if its length is 1, then -it is repeated using \code{rep} (see examples)} -} -\value{ -\code{\link{Pop-class}} -} -\description{ -Set miscelaneous information in a population -} -\details{ -A \code{NULL} in \code{value} is ignored -} diff --git a/man/smithHazel.Rd b/man/smithHazel.Rd index 7a78c485..2080ae53 100644 --- a/man/smithHazel.Rd +++ b/man/smithHazel.Rd @@ -17,7 +17,7 @@ smithHazel(econWt, varG, varP) a vector of weight for calculating index values } \description{ -Calculates weights for Smith-Hazel index given economice weights +Calculates weights for Smith-Hazel index given economice weights and phenotypic and genotypic variance-covariance matrices. } \examples{ diff --git a/man/transMat.Rd b/man/transMat.Rd index 4d7953ed..16451077 100644 --- a/man/transMat.Rd +++ b/man/transMat.Rd @@ -10,30 +10,30 @@ transMat(R) \item{R}{a correlation matrix} } \description{ -Creates an m by m linear transformation matrix that -can be applied to n by m uncorrelated deviates +Creates an m by m linear transformation matrix that +can be applied to n by m uncorrelated deviates sampled from a standard normal distribution to produce -correlated deviates with an arbitrary correlation -of R. If R is not positive semi-definite, the function +correlated deviates with an arbitrary correlation +of R. If R is not positive semi-definite, the function returns smoothing and returns a warning (see details). } \details{ -An eigendecomposition is applied to the correlation -matrix and used to test if it is positive semi-definite. -If the matrix is not positive semi-definite, it is not a -valid correlation matrix. In this case, smoothing is -applied to the matrix (as described in the 'cor.smooth' of -the 'psych' library) to obtain a valid correlation matrix. -The resulting deviates will thus not exactly match the -desired correlation, but will hopefully be close if the -input matrix wasn't too far removed from a valid +An eigendecomposition is applied to the correlation +matrix and used to test if it is positive semi-definite. +If the matrix is not positive semi-definite, it is not a +valid correlation matrix. In this case, smoothing is +applied to the matrix (as described in the 'cor.smooth' of +the 'psych' library) to obtain a valid correlation matrix. +The resulting deviates will thus not exactly match the +desired correlation, but will hopefully be close if the +input matrix wasn't too far removed from a valid correlation matrix. } \examples{ # Create an 2x2 correlation matrix R = 0.5*diag(2) + 0.5 -# Sample 1000 uncorrelated deviates from a +# Sample 1000 uncorrelated deviates from a # bivariate standard normal distribution X = matrix(rnorm(2*1000), ncol=2) diff --git a/man/usefulness.Rd b/man/usefulness.Rd index 0c7d178c..507bbdc1 100644 --- a/man/usefulness.Rd +++ b/man/usefulness.Rd @@ -15,24 +15,24 @@ usefulness( ) } \arguments{ -\item{pop}{and object of \code{\link{Pop-class}} or +\item{pop}{and object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}}} -\item{trait}{the trait for selection. Either a number indicating +\item{trait}{the trait for selection. Either a number indicating a single trait or a function returning a vector of length nInd.} \item{use}{select on genetic values (\code{gv}, default), estimated -breeding values (\code{ebv}), breeding values (\code{bv}), +breeding values (\code{ebv}), breeding values (\code{bv}), or phenotypes (\code{pheno})} \item{p}{the proportion of individuals selected} -\item{selectTop}{selects highest values if true. +\item{selectTop}{selects highest values if true. Selects lowest values if false.} \item{simParam}{an object of \code{\link{SimParam}}} -\item{...}{additional arguments if using a function for +\item{...}{additional arguments if using a function for trait} } \value{ @@ -53,7 +53,7 @@ SP$addTraitA(10) #Create population pop = newPop(founderPop, simParam=SP) -#Determine usefulness of population +#Determine usefulness of population usefulness(pop, simParam=SP) #Should be equivalent to GV of best individual diff --git a/man/varEBV.Rd b/man/varEBV.Rd new file mode 100644 index 00000000..6a5009a0 --- /dev/null +++ b/man/varEBV.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/popSummary.R +\name{varEBV} +\alias{varEBV} +\title{Variance of estimated breeding values} +\usage{ +varEBV(pop) +} +\arguments{ +\item{pop}{an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}}} +} +\description{ +Returns variance of estimated breeding values for all traits +} +\examples{ +#Create founder haplotypes +founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) + +#Set simulation parameters +SP = SimParam$new(founderPop) +SP$addTraitA(10) +trtH2 = 0.5 +SP$setVarE(h2=trtH2) +\dontshow{SP$nThreads = 1L} + +#Create population +pop = newPop(founderPop, simParam=SP) +pop@ebv = trtH2 * (pop@pheno - meanP(pop)) #ind performance based EBV +varA(pop) +varEBV(pop) + +} diff --git a/src/algorithm.cpp b/src/algorithm.cpp index 2a6164a8..2c2675b3 100644 --- a/src/algorithm.cpp +++ b/src/algorithm.cpp @@ -1219,9 +1219,9 @@ void GraphBuilder::build(){ // gene conversion stuff GeneConversionPtr newGC; double dLogTractRatio = log((pConfig->iGeneConvTract-1.)/pConfig->iGeneConvTract); - const int XOVER=0,GCSTART=1,GCEND=2; - int xovers=0,gcstarts=0,gcends=0; - int lastRE; + // const int XOVER=0,GCSTART=1,GCEND=2; + //int xovers=0,gcstarts=0,gcends=0; + //int lastRE; int iHistoryMax = 0; do{ if (iGraphIteration==0){ @@ -1232,14 +1232,14 @@ void GraphBuilder::build(){ // at this point decide whether we invoke a plain x-over // or a new gene conversion event this->bBeginGeneConversion = false; - lastRE=XOVER; + //lastRE=XOVER; if (this->bEndGeneConversion){ - lastRE=GCEND; + //lastRE=GCEND; }else{ this->bBeginGeneConversion = pRandNumGenerator->unifRV()< (pConfig->dGeneConvRatio/(pConfig->dGeneConvRatio+1))?true:false; if (bBeginGeneConversion){ - lastRE=GCSTART; + //lastRE=GCSTART; double dTractLen = (1.+log(pRandNumGenerator->unifRV())/ dLogTractRatio)/pConfig->dSeqLength; newGC = GeneConversionPtr(new GeneConversion( @@ -1248,17 +1248,17 @@ void GraphBuilder::build(){ } } invokeRecombination(newGC); - switch(lastRE){ - case XOVER: - ++xovers; - break; - case GCSTART: - ++gcstarts; - break; - case GCEND: - ++gcends; - break; - } + // switch(lastRE){ + // case XOVER: + // ++xovers; + // break; + // case GCSTART: + // ++gcstarts; + // break; + // case GCEND: + // ++gcends; + // break; + // } // mark the graph edges as the local tree markCurrentTree(); if (!bIncrementHistory){ diff --git a/src/calcGenParam.cpp b/src/calcGenParam.cpp index 4bd6731f..93b67f08 100644 --- a/src/calcGenParam.cpp +++ b/src/calcGenParam.cpp @@ -32,6 +32,9 @@ Rcpp::List calcGenParamE(const Rcpp::S4& trait, arma::vec genicAA2(nThreads,arma::fill::zeros); // No LD and HWE arma::vec mu(nThreads,arma::fill::zeros); // Observed mean arma::vec eMu(nThreads,arma::fill::zeros); // Expected mean with HWE + arma::vec alpha(a.n_elem); + arma::vec alphaHW(a.n_elem); + arma::mat ddMat, gv_d; if(hasD){ d = Rcpp::as(trait.slot("domEff")); @@ -54,7 +57,6 @@ Rcpp::List calcGenParamE(const Rcpp::S4& trait, for(arma::uword i=0; i genoMat = getGeno(Rcpp::as > >(pop.slot("geno")), lociPerChr, lociLoc, nThreads); @@ -336,7 +344,7 @@ Rcpp::List calcGenParam(const Rcpp::S4& trait, arma::vec aEff(ploidy+1), dEff(ploidy+1), eff(ploidy+1); // Genetic values, additive and dominance arma::vec bv(ploidy+1), dd(ploidy+1), gv(ploidy+1); // Statistical values, additive and dominance arma::vec bvE(ploidy+1), ddE(ploidy+1); //Expected for random mating - double gvMu, gvEMu, genoMu, p, q, dK, alpha, alphaE; + double gvMu, gvEMu, genoMu, p, q, dK; // Compute genotype frequencies for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, // Sample the exchange point arma::vec u(1, arma::fill::randu); double exchange = u(0)*genLen; + + // Sample start point for gamma model u.randu(); double start = u(0)-10; - // Determine crossover postions + // Determine crossover positions + // Returns field with crossover positions in each arm of the quadrivalent arma::field posCO = sampleQuadChiasmata(start, exchange, genLen, v, p); - // Set chromatid configuration for chiasmata - arma::field chromatidPairs(4); + // Set chromatid configuration for each chiasmata + arma::field chromatidPairs(4); // matches posCO for(arma::uword i=0; i<4; ++i){ + // Create table for chromatid pairs on an arm chromatidPairs(i).set_size(posCO(i).n_elem,2); + + // Assign pairs if there are chiasmata if(chromatidPairs(i).n_rows>0){ + // Initializing with "0" chromatid chromatidPairs(i).zeros(); + + // Randomly switch to "1" chromatid for(arma::uword j=0; j0.5){ @@ -400,6 +410,7 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, } } } + } // Allocate output with a naive maximum number of COs @@ -412,30 +423,38 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, output(1).set_size(maxCO+1,2); // Select centromeres (which chromosome and chromatid) + // Always taking chromosome 1 (1-4) and chromatid 1 (0-1) arma::uvec chromosome(2, arma::fill::ones); arma::uvec chromatid(2, arma::fill::ones); - chromosome(1) = sampleInt(1,3)(0) + 2; - chromatid(1) = sampleInt(1,2)(0); + chromosome(1) = sampleInt(1,3)(0) + 2; // 2-4 + chromatid(1) = sampleInt(1,2)(0); // 0-1 - // Loop through each of the selected centromeres + // Find starting chromosomes and chromatids by working backwards + // from selected centromeres to start of chromosome (head) arma::uword currentChromosome, currentChromatid; for(arma::uword i=0; i<2; ++i){ // Identify starting chromosome and chromatid currentChromosome = chromosome(i); currentChromatid = chromatid(i); - if(exchangecentromere){ // Centromere is in the head + if(currentChromosome<3){ // currentChromosome is 1 or 2 - // Account for all crossovers prior to the centromere + + // Loop through all chiasmata on arm 0 for(arma::uword j=posCO(0).n_elem; j>0; --j){ + + // Check if chiasmata is before centromere, ignore if not if(posCO(0)(j-1) > findQuadrivalentCO(const arma::vec& genMap, } } } + }else{ // currentChromosome is 3 or 4 - // Account for all crossovers prior to the centromere + + // Loop through all chiasmata on arm 2 for(arma::uword j=posCO(2).n_elem; j>0; --j){ + + // Check if chiasmata is before centromere, ignore if not if(posCO(2)(j-1) > findQuadrivalentCO(const arma::vec& genMap, } } } + }else{ // Centromere is in the tail - if((currentChromosome==1) | (currentChromosome==4)){ - // Find chromosome and chromatid before transition + + if((currentChromosome==1) | (currentChromosome==4)){ // Working on arm 3 + + // Find chromosome and chromatid before transition by looping + // through chiasmata on arm 3 for(arma::uword j=posCO(3).n_elem; j>0; --j){ + + // Check if chiasmata is before centromere, ignore if not if(posCO(3)(j-1) > findQuadrivalentCO(const arma::vec& genMap, } } } - // Find starting chromosome and chromatid + + // Find starting chromosome and chromatid by working back through head switch(currentChromosome){ - case 1: + case 1: // Work through arm 0 for(arma::uword j=posCO(0).n_elem; j>0; --j){ switch(currentChromosome){ - case 1: + case 1: // Check if there's a switch between chr 1 and 2 if(chromatidPairs(0)(j-1,0) == currentChromatid){ currentChromosome = 2; currentChromatid = chromatidPairs(0)(j-1,1); } break; - case 2: + case 2: // Check if there's a switch between chr 2 and 1 if(chromatidPairs(0)(j-1,1) == currentChromatid){ currentChromosome = 1; currentChromatid = chromatidPairs(0)(j-1,0); @@ -502,16 +533,16 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, } } break; - case 4: + case 4: // Work through arm 2 for(arma::uword j=posCO(2).n_elem; j>0; --j){ switch(currentChromosome){ - case 3: + case 3: // Check if there's a switch between chr 3 and 4 if(chromatidPairs(2)(j-1,0) == currentChromatid){ currentChromosome = 4; currentChromatid = chromatidPairs(2)(j-1,1); } break; - case 4: + case 4: // Check if there's a switch between chr 4 and 3 if(chromatidPairs(2)(j-1,1) == currentChromatid){ currentChromosome = 3; currentChromatid = chromatidPairs(2)(j-1,0); @@ -520,18 +551,22 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, } } - }else{ // currentChromosome is 2 or 3 - // Find chromosome and chromatid before transition + }else{ // Working on arm 1 + + // Find chromosome and chromatid before transition by looping + // through chiasmata on arm 1 for(arma::uword j=posCO(1).n_elem; j>0; --j){ + + // Check if chiasmata is before centromere, ignore if not if(posCO(1)(j-1) > findQuadrivalentCO(const arma::vec& genMap, } } } - // Find starting chromosome and chromatid + + // Find starting chromosome and chromatid by working back through head switch(currentChromosome){ - case 2: + case 2: // Work through arm 0 for(arma::uword j=posCO(0).n_elem; j>0; --j){ switch(currentChromosome){ - case 1: + case 1: // Check if there's a switch between chr 1 and 2 if(chromatidPairs(0)(j-1,0) == currentChromatid){ currentChromosome = 2; currentChromatid = chromatidPairs(0)(j-1,1); } break; - case 2: + case 2: // Check if there's a switch between chr 2 and 1 if(chromatidPairs(0)(j-1,1) == currentChromatid){ currentChromosome = 1; currentChromatid = chromatidPairs(0)(j-1,0); @@ -558,16 +594,16 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, } } break; - case 3: + case 3: // Work through arm 2 for(arma::uword j=posCO(2).n_elem; j>0; --j){ switch(currentChromosome){ - case 3: + case 3: // Check if there's a switch between chr 3 and 4 if(chromatidPairs(2)(j-1,0) == currentChromatid){ currentChromosome = 4; currentChromatid = chromatidPairs(2)(j-1,1); } break; - case 4: + case 4: // Check if there's a switch between chr 4 and 3 if(chromatidPairs(2)(j-1,1) == currentChromatid){ currentChromosome = 3; currentChromatid = chromatidPairs(2)(j-1,0); @@ -578,15 +614,19 @@ arma::field > findQuadrivalentCO(const arma::vec& genMap, } } - // Fill in crossover map + // Fill in crossover map by working from head to tail arma::uword startPos=0, endPos, nCO=0; output(i)(0,0) = currentChromosome; output(i)(0,1) = 1; - if(currentChromosome<3){ + + if(currentChromosome<3){ // Start in arm 0 + // Fill crossovers in the head for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, } } } + // Fill crossovers in the tail - if(currentChromosome==1){ + if(currentChromosome==1){ // Move to arm 3 for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, } } } - }else{ // currentChromosome = 2 + }else{ // currentChromosome = 2, move to arm 1 for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, } } } - }else{ + }else{ // currentChromosome>2, start in arm 2 + // Fill crossovers in the head for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, } } } + // Fill crossovers in the tail - if(currentChromosome==4){ + if(currentChromosome==4){ // Move to arm 3 for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, } } } - }else{ // currentChromosome = 3 + }else{ // currentChromosome = 3, move to arm 1 for(arma::uword j=0; j > findQuadrivalentCO(const arma::vec& genMap, output(i) = output(i).rows(arma::span(0,nCO)); output(i) = removeDoubleCO(output(i)); } + return output; } diff --git a/tests/testthat/test-misc.R b/tests/testthat/test-misc.R new file mode 100644 index 00000000..02e7e41d --- /dev/null +++ b/tests/testthat/test-misc.R @@ -0,0 +1,42 @@ +context("misc") + +test_that("misc_and_miscPop",{ + founderPop = quickHaplo(nInd=2, nChr=1, segSites=10) + SP = SimParam$new(founderPop) + SP$addTraitA(10) + popOrig = newPop(founderPop, simParam=SP) + multiPop = newMultiPop(popOrig, popOrig) + + popSub = popOrig[1] + expect_equal(popSub@misc, list()) + expect_equal(popSub@miscPop, list()) + + pop = popOrig + pop@misc$vec = rnorm(n=2) + pop@misc$mtP = popOrig # hmm, should this actually be multiple pop objects or one with multiple individuals? + pop@misc$mtLP = list(popOrig, popOrig) + pop@misc$mtMP = multiPop + pop@miscPop$vec = sum(pop@misc$vec) + pop@miscPop$af = colMeans(pullSegSiteGeno(pop, simParam=SP)) + + popSub = pop[1] + expect_equal(popSub@misc$vec, pop@misc$vec[1]) + expect_equal(popSub@misc$mtP, pop@misc$mtP[1]) + expect_equal(popSub@misc$mtLP, pop@misc$mtLP[1]) + expect_equal(popSub@misc$mtMP, pop@misc$mtMP[1]) + expect_equal(popSub@miscPop, list()) + + popSub = pop[0] + expect_equal(popSub@misc$vec, numeric(0)) + expect_equal(popSub@misc$mtP, newEmptyPop(simParam=SP)) + expect_equal(popSub@misc$mtLP, list()) + expect_equal(popSub@misc$mtMP, new("MultiPop", pops=list())) + expect_equal(popSub@miscPop, list()) + + popC = c(pop, pop) + expect_equal(popC@misc$vec, c(pop@misc$vec, pop@misc$vec)) + expect_equal(popC@misc$mtP, c(pop@misc$mtP, pop@misc$mtP)) + expect_equal(popC@misc$mtLP, c(pop@misc$mtLP, pop@misc$mtLP)) + expect_equal(popC@misc$mtMP, c(pop@misc$mtMP, pop@misc$mtMP)) + expect_equal(popC@miscPop, list()) +}) diff --git a/tests/testthat/test-selection.R b/tests/testthat/test-selection.R new file mode 100644 index 00000000..3f4e1e3a --- /dev/null +++ b/tests/testthat/test-selection.R @@ -0,0 +1,37 @@ +context("selection") + +test_that("selectInd_and_getResponse",{ + founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) + SP = SimParam$new(founderPop) + SP$addTraitA(10) + SP$setVarE(h2=0.5) + pop = newPop(founderPop, simParam=SP) + + pop2 = selectInd(pop, 5, simParam=SP) + expect_equal(pop2@id, + pop[order(pop@pheno, decreasing=TRUE)[1:5]]@id) + + pop2b = selectInd(pop, 5, trait="Trait1", simParam=SP) + expect_equal(pop2@id, + pop2b@id) + + squaredDeviation = function(x, optima=0) (x - optima)^2 + pop3 = selectInd(pop, 5, trait=squaredDeviation, selectTop=TRUE, simParam=SP) + expect_equal(pop3@id, + pop[order(squaredDeviation(pop@pheno), decreasing=TRUE)[1:5]]@id) + + pop4 = selectInd(pop, 5, trait=squaredDeviation, selectTop=FALSE, simParam=SP) + expect_equal(pop4@id, + pop[order(squaredDeviation(pop@pheno), decreasing=FALSE)[1:5]]@id) + + pop@misc = list(smth=rnorm(10), smth2=rnorm(10)) + useFunc = function(pop, trait=NULL) pop@misc$smth + pop@misc$smth2 + pop5 = selectInd(pop, 5, use=useFunc, simParam=SP) + expect_equal(pop5@id, + pop[order(useFunc(pop), decreasing=TRUE)[1:5]]@id) + + useFunc2 = function(pop, trait=NULL) cbind(pop@misc$smth, pop@misc$smth2) + trtFunc = function(x) rowSums(x) + pop6 = selectInd(pop, 5, trait=trtFunc, use=useFunc2, simParam=SP) + expect_equal(pop5@id, pop6@id) +}) diff --git a/vignettes/articles/AlphaSimR.Rmd b/vignettes/articles/AlphaSimR.Rmd new file mode 100644 index 00000000..c56fbb7d --- /dev/null +++ b/vignettes/articles/AlphaSimR.Rmd @@ -0,0 +1,25 @@ +--- +title: "AlphaSimR" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Explain the package + +Point to important pages for further information (introduction) + +Include links to external sources diff --git a/vignettes/articles/DataImport.Rmd b/vignettes/articles/DataImport.Rmd new file mode 100644 index 00000000..10f55d57 --- /dev/null +++ b/vignettes/articles/DataImport.Rmd @@ -0,0 +1,21 @@ +--- +title: "Data Import" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Pull from example script diff --git a/vignettes/articles/Genome.Rmd b/vignettes/articles/Genome.Rmd new file mode 100644 index 00000000..b09d0b15 --- /dev/null +++ b/vignettes/articles/Genome.Rmd @@ -0,0 +1,26 @@ +--- +title: "Genome" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Introduction + Pull from MOOC + Explain the genome and its representation + Haplotypes / genotypes + QTL/SNP + Genetic map vs physical map diff --git a/vignettes/articles/GenomicSelection.Rmd b/vignettes/articles/GenomicSelection.Rmd new file mode 100644 index 00000000..5acee4e3 --- /dev/null +++ b/vignettes/articles/GenomicSelection.Rmd @@ -0,0 +1,39 @@ +--- +title: "Genomic Selection" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Introduction + Point to genome (QTL/SNP) and trait + +How to build a training population + Mention fixed effect + +How to fit a model + +How to get GEBVs + Accuracy of selection (vs BV or GV, train vs pred) + +How to do selections + +Explain models + Use Guilherme's paper + Explain 1 vs 2 + +How to use external software + diff --git a/vignettes/articles/Meiosis.Rmd b/vignettes/articles/Meiosis.Rmd new file mode 100644 index 00000000..eb574454 --- /dev/null +++ b/vignettes/articles/Meiosis.Rmd @@ -0,0 +1,43 @@ +--- +title: "Meiosis" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Introduction + Pull from MOOC + Explain the genome and its representation + Recombination / meiosis + Stages of meiosis + Chromosome pairing + Chiasma (crossovers) + Crossover interference + Segregation + Diploids first and then polyploids + Centromeres for polyploids + + Implementation + Gamma sprinkling model + Crossover interference + recHist + pedigree + Sex specific recombination rates + + + + + diff --git a/vignettes/articles/Misc.Rmd b/vignettes/articles/Misc.Rmd new file mode 100644 index 00000000..bcb3b982 --- /dev/null +++ b/vignettes/articles/Misc.Rmd @@ -0,0 +1,31 @@ +--- +title: "Misc" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Introduction + +Explain individual level slot + set/get + +Explain population slot + +Finalize pop? + + + diff --git a/vignettes/articles/QuantGen.Rmd b/vignettes/articles/QuantGen.Rmd new file mode 100644 index 00000000..f6bdc78f --- /dev/null +++ b/vignettes/articles/QuantGen.Rmd @@ -0,0 +1,32 @@ +--- +title: "Quantitative Genetics" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Introduction + Point to vignette on traits as first step + History + Fisher (response to selection, variance components) + Bulmer (genic variance) + +Decomposition of genetic values + Reference population (HWE vs non-HWE) + Breeding value + Dominance deviations + Epistatic deviations + Variance, genic variance, LD covariance diff --git a/vignettes/articles/Traits.Rmd b/vignettes/articles/Traits.Rmd new file mode 100644 index 00000000..41685812 --- /dev/null +++ b/vignettes/articles/Traits.Rmd @@ -0,0 +1,35 @@ +--- +title: "Traits" +author: "Vignette Author" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Vignette Title} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Placeholder for future development + +Background from existing vignette + altAddTraitAD (its own section?) + +Practical examples + How to get effects + Manual adjustments + importation + Mixture distributions + rescaleTrait + + + + + +