From 13ac861444a5c82d4af582ca80c3f455336313e0 Mon Sep 17 00:00:00 2001 From: david20011999 Date: Wed, 1 Jul 2026 17:50:11 +0200 Subject: [PATCH] Implementation of imprinting simulation in new AlphaSimR version --- R/Class-LociMap.R | 62 ++++++ R/Class-SimParam.R | 275 +++++++++++++++++++++++- R/popSummary.R | 308 +++++++++++++++++++++++++-- src/calcGenParam.cpp | 370 ++++++++++++++++++++++++++++++++- src/getGv.cpp | 169 ++++++++++++++- tests/testthat/test-addTrait.R | 43 ++++ vignettes/TestAI.Rmd | 195 +++++++++++++++++ 7 files changed, 1386 insertions(+), 36 deletions(-) create mode 100644 vignettes/TestAI.Rmd diff --git a/R/Class-LociMap.R b/R/Class-LociMap.R index bb8f12e8..d4ffea72 100644 --- a/R/Class-LociMap.R +++ b/R/Class-LociMap.R @@ -170,6 +170,68 @@ isTraitAD = function(x) { return(ret) } +#TraitAI---- +#' @title Additive and imprinting trait +#' +#' @description Extends \code{\link{TraitA-class}} +#' to add imprinting +#' +#' @slot impEff imprinting effects +#' +#' @export +setClass("TraitAI", + slots=c(impEff="numeric"), + contains="TraitA") + +setValidity("TraitAI",function(object){ + errors = character() + if(object@nLoci!=length(object@impEff)){ + errors = c(errors,"nLoci!=length(impEff)") + } + if(length(errors)==0){ + return(TRUE) + }else{ + return(errors) + } +}) + +# Test if object is of a TraitAI class +isTraitAI = function(x) { + ret = is(x, class2 = "TraitAI") + return(ret) +} + +#TraitADI---- +#' @title Additive dominance and imprinting trait +#' +#' @description Extends \code{\link{TraitAD-class}} +#' to add imprinting +#' +#' @slot impEff imprinting effects +#' +#' @export +setClass("TraitADI", + slots=c(impEff="numeric"), + contains="TraitAD") + +setValidity("TraitADI",function(object){ + errors = character() + if(object@nLoci!=length(object@impEff)){ + errors = c(errors,"nLoci!=length(impEff)") + } + if(length(errors)==0){ + return(TRUE) + }else{ + return(errors) + } +}) + +# Test if object is of a TraitAI class +isTraitADI = function(x) { + ret = is(x, class2 = "TraitADI") + return(ret) +} + #TraitA2D---- #' @title Sex specific additive and dominance trait #' diff --git a/R/Class-SimParam.R b/R/Class-SimParam.R index 19449b0a..03300f6a 100644 --- a/R/Class-SimParam.R +++ b/R/Class-SimParam.R @@ -637,6 +637,202 @@ SimParam = R6Class( } invisible(self) }, + + #' @description + #' Randomly assigns eligible QTLs for one or more traits with imprinting (silencing) + #' If simulating more than one trait, all traits will be pleiotropic + #' with correlated effects. + #' + #' @param nQtlPerChr number of QTLs per chromosome. Can be a single value or nChr values. + #' @param mean a vector of desired mean genetic values for one or more traits + #' @param var a vector of desired genetic variances for one or more traits + #' @param meanID mean imprinting degree + #' @param varID variance of imprinting degree + #' @param corA a matrix of correlations between additive effects + #' @param corID a matrix of correlations between imprinting degrees + #' @param useVarA tune according to additive genetic variance if true. If + #' FALSE, tuning is performed according to total genetic variance. + #' @param gamma should a gamma distribution be used instead of normal + #' @param shape the shape parameter for the gamma distribution + #' (the rate/scale parameter of the gamma distribution is accounted + #' for via the desired level of genetic variance, the var argument) + #' @param force should the check for a running simulation be + #' ignored. Only set to TRUE if you know what you are doing. + #' @param name optional name for trait(s) + #' + #' @examples + #' #Create founder haplotypes + #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=1) + #' + #' #Set simulation parameters + #' SP = SimParam$new(founderPop) + #' \dontshow{SP$nThreads = 1L} + #' SP$addTraitAI(1, meanID=0.5) + addTraitAI = function(nQtlPerChr,mean=0,var=1,meanID=0, + varID=0,corA=NULL,corID=NULL,useVarA=TRUE, + gamma=FALSE,shape=1,force=FALSE,name=NULL){ + if(self$founderPop@ploidy > 2){ + stop("ERROR: Imprinting is not yet implemented for polyploids!") + } + if(!force){ + private$.isRunning() + } + if(length(nQtlPerChr)==1){ + nQtlPerChr = rep(nQtlPerChr,self$nChr) + } + nTraits = length(mean) + if(length(meanID)==1) meanID = rep(meanID,nTraits) + if(length(varID)==1) varID = rep(varID,nTraits) + if(length(gamma)==1) gamma = rep(gamma,nTraits) + if(length(shape)==1) shape = rep(shape,nTraits) + if(is.null(corA)) corA=diag(nTraits) + if(is.null(corID)) corID=diag(nTraits) + if(is.null(name)){ + name = paste0("Trait",1:nTraits+self$nTraits) + } + stopifnot(length(mean)==length(var), + length(meanID)==length(mean), + isSymmetric(corA), + isSymmetric(corID), + nrow(corA)==nTraits, + nrow(corID)==nTraits, + length(varID)==nTraits, + length(name)==nTraits) + qtlLoci = private$.pickLoci(nQtlPerChr) + addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, + corr=corA,gamma=gamma,shape=shape) + impEff = sampImpEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, + corID=corID,meanID=meanID,varID=varID) + for(i in 1:nTraits){ + trait = new("TraitAI", + qtlLoci, + addEff=addEff[,i], + impEff=impEff[,i], + intercept=0, + name=name[i]) + tmp = calcGenParam(trait, self$founderPop, + self$nThreads) + if(useVarA){ + scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1]) + }else{ + scale = sqrt(var[i])/sqrt(popVar(tmp$gv)[1]) + } + trait@addEff = trait@addEff*scale + trait@impEff = trait@impEff*scale + trait@intercept = mean[i]-mean(tmp$gv*scale) + if(useVarA){ + private$.addTrait(trait,var[i],popVar(tmp$gv*scale)[1]) + }else{ + private$.addTrait(trait,popVar(tmp$bv*scale)[1],var[i]) + } + } + invisible(self) + }, + + #' @description + #' Randomly assigns eligible QTLs for one or more traits with dominance and imprinting (silencing) + #' If simulating more than one trait, all traits will be pleiotropic + #' with correlated effects. + #' + #' @param nQtlPerChr number of QTLs per chromosome. Can be a single value or nChr values. + #' @param mean a vector of desired mean genetic values for one or more traits + #' @param var a vector of desired genetic variances for one or more traits + #' @param meanDD mean dominance degree + #' @param varDD variance of dominance degree + #' @param meanID mean imprinting degree + #' @param varID variance of imprinting degree + #' @param corA a matrix of correlations between additive effects + #' @param corDD a matrix of correlations between dominance degrees + #' @param corID a matrix of correlations between imprinting degrees + #' @param useVarA tune according to additive genetic variance if true. If + #' FALSE, tuning is performed according to total genetic variance. + #' @param gamma should a gamma distribution be used instead of normal + #' @param shape the shape parameter for the gamma distribution + #' (the rate/scale parameter of the gamma distribution is accounted + #' for via the desired level of genetic variance, the var argument) + #' @param force should the check for a running simulation be + #' ignored. Only set to TRUE if you know what you are doing. + #' @param name optional name for trait(s) + #' + #' @examples + #' #Create founder haplotypes + #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=1) + #' + #' #Set simulation parameters + #' SP = SimParam$new(founderPop) + #' \dontshow{SP$nThreads = 1L} + #' SP$addTraitADI(1, meanDD=0.5, meanID=0.5) + addTraitADI = function(nQtlPerChr,mean=0,var=1,meanDD=0,varDD=0,meanID=0, + varID=0,corA=NULL,corDD=NULL,corID=NULL,useVarA=TRUE, + gamma=FALSE,shape=1,force=FALSE,name=NULL){ + if(self$founderPop@ploidy > 2){ + stop("ERROR: Imprinting is not yet implemented for polyploids!") + } + if(!force){ + private$.isRunning() + } + if(length(nQtlPerChr)==1){ + nQtlPerChr = rep(nQtlPerChr,self$nChr) + } + nTraits = length(mean) + if(length(meanDD)==1) meanDD = rep(meanDD,nTraits) + if(length(meanID)==1) meanID = rep(meanID,nTraits) + if(length(varDD)==1) varDD = rep(varDD,nTraits) + if(length(varID)==1) varID = rep(varID,nTraits) + if(length(gamma)==1) gamma = rep(gamma,nTraits) + if(length(shape)==1) shape = rep(shape,nTraits) + if(is.null(corA)) corA=diag(nTraits) + if(is.null(corDD)) corDD=diag(nTraits) + if(is.null(corID)) corID=diag(nTraits) + if(is.null(name)){ + name = paste0("Trait",1:nTraits+self$nTraits) + } + stopifnot(length(mean)==length(var), + length(meanDD)==length(mean), + length(meanID)==length(mean), + isSymmetric(corA), + isSymmetric(corDD), + isSymmetric(corID), + nrow(corA)==nTraits, + nrow(corDD)==nTraits, + nrow(corID)==nTraits, + length(varDD)==nTraits, + length(varID)==nTraits, + length(name)==nTraits) + qtlLoci = private$.pickLoci(nQtlPerChr) + addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, + corr=corA,gamma=gamma,shape=shape) + domEff = sampDomEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, + corDD=corDD,meanDD=meanDD,varDD=varDD) + impEff = sampImpEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, + corID=corID,meanID=meanID,varID=varID) + for(i in 1:nTraits){ + trait = new("TraitADI", + qtlLoci, + addEff=addEff[,i], + domEff=domEff[,i], + impEff=impEff[,i], + intercept=0, + name=name[i]) + tmp = calcGenParam(trait, self$founderPop, + self$nThreads) + if(useVarA){ + scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1]) + }else{ + scale = sqrt(var[i])/sqrt(popVar(tmp$gv)[1]) + } + trait@addEff = trait@addEff*scale + trait@domEff = trait@domEff*scale + trait@impEff = trait@impEff*scale + trait@intercept = mean[i]-mean(tmp$gv*scale) + if(useVarA){ + private$.addTrait(trait,var[i],popVar(tmp$gv*scale)[1]) + }else{ + private$.addTrait(trait,popVar(tmp$bv*scale)[1],var[i]) + } + } + invisible(self) + }, #' @description #' An alternative method for adding a trait with additive and dominance effects @@ -1458,13 +1654,14 @@ SimParam = R6Class( #' formatting the trait as a \code{\link{LociMap-class}}. #' The formatting is performed automatically for the user, #' with more user friendly data.frames or matrices taken as - #' inputs. This function only works for A and AD trait types. + #' inputs. This function only works for A, AD, AI and ADI trait types. #' #' @param markerNames a vector of names for the QTL #' @param addEff a matrix of additive effects (nLoci x nTraits). #' Alternatively, a vector of length nLoci can be supplied for #' a single trait. #' @param domEff optional dominance effects for each locus + #' @param impEff optional imprinting effects for each locus #' @param intercept optional intercepts for each trait #' @param name optional name(s) for the trait(s) #' @param varE default error variance for phenotype, optional @@ -1473,6 +1670,7 @@ SimParam = R6Class( importTrait = function(markerNames, addEff, domEff=NULL, + impEff=NULL, intercept=NULL, name=NULL, varE=NULL, @@ -1493,6 +1691,16 @@ SimParam = R6Class( stopifnot(nrow(addEff)==nrow(domEff), ncol(addEff)==ncol(domEff)) } + if(is.null(impEff)){ + useImp = FALSE + }else{ + useImp = TRUE + impEff = as.matrix(impEff) + stopifnot(nrow(addEff)==nrow(impEff), + ncol(addEff)==ncol(impEff), + nrow(domEff)==nrow(impEff), + ncol(domEff)==ncol(impEff)) + } # Prepare the intercept if(is.null(intercept)){ @@ -1524,16 +1732,16 @@ SimParam = R6Class( # Create trait variables lociPerChr = integer(self$nChr) lociLoc = vector("list", self$nChr) - addEffList = domEffList = vector("list", nTraits) + addEffList = domEffList = impEffList = vector("list", nTraits) for(i in seq_len(nTraits)){ - addEffList[[i]] = domEffList[[i]] = vector("list", self$nChr) + addEffList[[i]] = domEffList[[i]] = impEffList[[i]] = vector("list", self$nChr) } # Loop through chromosomes for(i in seq_len(self$nChr)){ # Working on trait 1 # Initialize variables - addEffList[[1]][[i]] = domEffList[[1]][[i]] = numeric() + addEffList[[1]][[i]] = domEffList[[1]][[i]] = impEffList[[1]][[i]] = numeric() lociLoc[[i]] = integer() # Find matches if they exist @@ -1546,17 +1754,23 @@ SimParam = R6Class( if(useDom){ domEffList[[1]][[i]] = domEff[na.omit(take),1] } + if(useImp){ + impEffList[[1]][[i]] = impEff[na.omit(take),1] + } } # Work on additional traits? if(nTraits>1){ for(j in 2:nTraits){ - addEffList[[j]][[i]] = domEffList[[j]][[i]] = numeric() + addEffList[[j]][[i]] = domEffList[[j]][[i]] = impEffList[[j]][[i]] = numeric() if(lociPerChr[i]>0L){ addEffList[[j]][[i]] = addEff[na.omit(take),j] if(useDom){ domEffList[[j]][[i]] = domEff[na.omit(take),j] } + if(useImp){ + impEffList[[j]][[i]] = impEff[na.omit(take),j] + } } } } @@ -1570,6 +1784,18 @@ SimParam = R6Class( addEff = unlist(addEffList[[i]]) if(useDom){ domEff = unlist(domEffList[[i]]) + if(useImp){ + impEff = unlist(impEffList[[i]]) + trait = new("TraitADI", + addEff=addEff, + domEff=domEff, + impEff=impEff, + intercept=intercept[i], + nLoci=nLoci, + lociPerChr=lociPerChr, + lociLoc=lociLoc, + name=name[i]) + }else{ trait = new("TraitAD", addEff=addEff, domEff=domEff, @@ -1578,7 +1804,19 @@ SimParam = R6Class( lociPerChr=lociPerChr, lociLoc=lociLoc, name=name[i]) + } }else{ + if(useImp){ + impEff = unlist(impEffList[[i]]) + trait = new("TraitAI", + addEff=addEff, + impEff=impEff, + intercept=intercept[i], + nLoci=nLoci, + lociPerChr=lociPerChr, + lociLoc=lociLoc, + name=name[i]) + }else{ trait = new("TraitA", addEff=addEff, intercept=intercept[i], @@ -1586,6 +1824,7 @@ SimParam = R6Class( lociPerChr=lociPerChr, lociLoc=lociLoc, name=name[i]) + } } # Add trait to simParam @@ -2717,6 +2956,32 @@ sampDomEff = function(qtlLoci,nTraits,addEff,corDD, return(domEff) } +#' @title Sample imprinting effects +#' +#' @description Samples imprinting deviation effects from a normal distribution +#' and uses previously sampled additive effects to form imprinting +#' effects +#' +#' @param qtlLoci total number of loci +#' @param nTraits number of traits +#' @param addEff previously sampled additive effects +#' @param corID correlation between imprinting degrees +#' @param meanID mean value of imprinting degrees +#' @param varID variance of imprinting degrees +#' +#' @returns a matrix with dimensions qtlLoci by nTraits +#' +#' @keywords internal +sampImpEff = function(qtlLoci,nTraits,addEff,corID, + meanID,varID){ + impEff = matrix(rnorm(qtlLoci@nLoci*nTraits), + ncol=nTraits)%*%transMat(corID) + impEff = sweep(impEff,2,sqrt(varID),"*") + impEff = sweep(impEff,2,meanID,"+") + impEff = abs(addEff)*impEff + return(impEff) +} + #' @title Sample epistatic effects #' #' @description Samples epistatic effects from a normal distribution or gamma distribution diff --git a/R/popSummary.R b/R/popSummary.R index c53c8a2f..6a40b7b5 100644 --- a/R/popSummary.R +++ b/R/popSummary.R @@ -79,7 +79,8 @@ meanEBV = function(pop){ #' @title Total genetic variance #' -#' @description Returns total genetic variance for all traits +#' @description Returns total genetic variance (=variance of genetic values) +#' for all traits #' #' @param pop an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}} #' @@ -106,7 +107,8 @@ varG = function(pop){ #' @title Phenotypic variance #' -#' @description Returns phenotypic variance for all traits +#' @description Returns phenotypic variance (=variance of phenotypic values) +#' for all traits #' #' @param pop an object of \code{\link{Pop-class}} or \code{\link{HybridPop-class}} #' @@ -164,8 +166,7 @@ varEBV = function(pop){ #' @title Sumarize genetic parameters #' #' @description -#' Calculates genetic and genic additive and dominance variances -#' for an object of \code{\link{Pop-class}} +#' Calculates genetic values and variances for an object of \code{\link{Pop-class}} #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -174,32 +175,41 @@ varEBV = function(pop){ #' \describe{ #' \item{varA}{an nTrait by nTrait matrix of additive genetic variances} #' \item{varD}{an nTrait by nTrait matrix of dominance genetic variances} +#' \item{varI}{an nTrait by nTrait matrix of imprinting genetic variances} #' \item{varAA}{an nTrait by nTrait matrix of additive-by-additive genetic variances} #' \item{varG}{an nTrait by nTrait matrix of total genetic variances} #' \item{genicVarA}{an nTrait vector of additive genic variances} #' \item{genicVarD}{an nTrait vector of dominance genic variances} +#' \item{genicVarI}{an nTrait vector of imprinting genic variances} #' \item{genicVarAA}{an nTrait vector of additive-by-additive genic variances} #' \item{genicVarG}{an nTrait vector of total genic variances} #' \item{covA_HW}{an nTrait vector of additive covariances due to non-random mating} #' \item{covD_HW}{an nTrait vector of dominance covariances due to non-random mating} +#' \item{covI_HW}{an nTrait vector of imprinting covariances due to non-random mating} #' \item{covAA_HW}{an nTrait vector of additive-by-additive covariances due to non-random mating} #' \item{covG_HW}{an nTrait vector of total genic covariances due to non-random mating} #' \item{covA_L}{an nTrait vector of additive covariances due to linkage disequilibrium} #' \item{covD_L}{an nTrait vector of dominance covariances due to linkage disequilibrium} +#' \item{covI_L}{an nTrait vector of imprinting covariances due to linkage disequilibrium} #' \item{covAA_L}{an nTrait vector of additive-by-additive covariances due to linkage disequilibrium} #' \item{covAD_L}{an nTrait vector of additive by dominance covariances due to linkage disequilibrium} +#' \item{covAI_L}{an nTrait vector of additive by imprinting covariances due to linkage disequilibrium} #' \item{covAAA_L}{an nTrait vector of additive by additive-by-additive covariances due to linkage disequilibrium} #' \item{covDAA_L}{an nTrait vector of dominance by additive-by-additive covariances due to linkage disequilibrium} +#' \item{covIAA_L}{an nTrait vector of imprinting by additive-by-additive covariances due to linkage disequilibrium} #' \item{covG_L}{an nTrait vector of total genic covariances due to linkage disequilibrium} #' \item{mu}{an nTrait vector of trait means} #' \item{mu_HW}{an nTrait vector of expected trait means under random mating} #' \item{gv}{a matrix of genetic values with dimensions nInd by nTraits} #' \item{bv}{a matrix of breeding values with dimensions nInd by nTraits} #' \item{dd}{a matrix of dominance deviations with dimensions nInd by nTraits} +#' \item{idM}{a matrix of maternal imprinting deviations with dimensions nInd by nTraits +#' (paternal imprinting deviations are \code{-idM})} #' \item{aa}{a matrix of additive-by-additive epistatic deviations with dimensions nInd by nTraits} #' \item{gv_mu}{an nTrait vector of intercepts with dimensions nInd by nTraits} #' \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_i}{a matrix of imprinting 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} @@ -232,14 +242,14 @@ genParam = function(pop,simParam=NULL){ # Blank nInd x nTrait matrices gv = matrix(NA_real_, nrow=nInd, ncol=nTraits) colnames(gv) = traitNames - bv = dd = aa = gv_a = gv_d = gv_aa = gv + bv = bvM = bvP = dd = id = aa = gv_a = gv_d = gv_i = gv_aa = gv # Blank nTrait vectors genicVarA = rep(NA_real_, nTraits) names(genicVarA) = traitNames - genicVarD = genicVarAA = covA_HW = covD_HW = covAA_HW = - covG_HW = mu = mu_HW = gv_mu = covAAA_L = covDAA_L = - covAD_L = genicVarA + genicVarD = genicVarI = genicVarAA = covA_HW = covD_HW = covI_HW = covAA_HW = + covG_HW = mu = mu_HW = gv_mu = covAAA_L = covDAA_L = covIAA_L = + covAD_L = covAI_L = covDI_L = genicVarA # Average effect of an allele substitution alpha = vector("list", length=nTraits) @@ -269,6 +279,21 @@ genParam = function(pop,simParam=NULL){ dd[,i] = rep(0,pop@nInd) gv_d[,i] = rep(0,pop@nInd) } + if(.hasSlot(trait,"impEff")){ + genicVarI[i] = tmp$genicVarI2 + covI_HW[i] = tmp$genicVarI-tmp$genicVarI2 + id[,i] = tmp$id + gv_i[,i] = tmp$gv_i + bvM[,i] = tmp$bvM + bvP[,i] = tmp$bvP + }else{ + genicVarI[i] = 0 + covI_HW[i] = 0 + id[,i] = rep(0,pop@nInd) + gv_i[,i] = rep(0,pop@nInd) + bvM[,i] = rep(0,pop@nInd) + bvP[,i] = rep(0,pop@nInd) + } if(.hasSlot(trait,"epiEff")){ genicVarAA[i] = tmp$genicVarAA2 covAA_HW[i] = tmp$genicVarAA-tmp$genicVarAA2 @@ -282,60 +307,83 @@ genParam = function(pop,simParam=NULL){ } if(nInd==1){ covAD_L[i] = 0 + covAI_L[i] = 0 + covDI_L[i] = 0 covAAA_L[i] = 0 covDAA_L[i] = 0 + covIAA_L[i] = 0 } else { covAD_L[i] = popVar(cbind(bv[,i],dd[,i]))[1,2] + covAI_L[i] = popVar(cbind(bv[,i],id[,i]))[1,2] + covDI_L[i] = popVar(cbind(dd[,i],id[,i]))[1,2] covAAA_L[i] = popVar(cbind(bv[,i],aa[,i]))[1,2] covDAA_L[i] = popVar(cbind(dd[,i],aa[,i]))[1,2] + covIAA_L[i] = popVar(cbind(id[,i],aa[,i]))[1,2] + # TODO: do we need covIMA and covIPA etc.? + # TODO: do we need covIMA and covIPA etc.? } alpha[[i]] = tmp$alpha alpha_HW[[i]] = tmp$alpha_HW } + # TODO: do we need varAM and varAP? varA = popVar(bv) rownames(varA) = colnames(varA) = traitNames varD = popVar(dd) rownames(varD) = colnames(varD) = traitNames + # TODO: do we need varIM and varIP (they are the same in a "random" setting)? + varI = popVar(id) + rownames(varI) = colnames(varI) = traitNames + varAA = popVar(aa) rownames(varAA) = colnames(varAA) = traitNames varG = popVar(gv) rownames(varG) = colnames(varG) = traitNames - genicVarG = genicVarA + genicVarD + genicVarAA - covG_HW = covA_HW + covD_HW + covAA_HW + genicVarG = genicVarA + genicVarD + genicVarI + genicVarAA + covG_HW = covA_HW + covD_HW + covI_HW + covAA_HW output = list(varA=varA, varD=varD, + varI=varI, varAA=varAA, varG=varG, genicVarA=genicVarA, genicVarD=genicVarD, + genicVarI=genicVarI, genicVarAA=genicVarAA, genicVarG=genicVarG, covA_HW=covA_HW, covD_HW=covD_HW, + covI_HW=covI_HW, covAA_HW=covAA_HW, covG_HW=covG_HW, covA_L=diag(varA)-genicVarA-covA_HW, covD_L=diag(varD)-genicVarD-covD_HW, + covI_L=diag(varI)-genicVarI-covI_HW, covAA_L=diag(varAA)-genicVarAA-covAA_HW, covAD_L=covAD_L, + covAI_L=covAI_L, covAAA_L=covAAA_L, covDAA_L=covDAA_L, + covIAA_L=covIAA_L, covG_L=diag(varG)-genicVarG-covG_HW, mu=mu, mu_HW=mu_HW, gv=gv, bv=bv, + bvM=bvM, + bvP=bvP, dd=dd, + id=id, aa=aa, gv_mu=gv_mu, gv_a=gv_a, gv_d=gv_d, + gv_i=gv_i, gv_aa=gv_aa, alpha=alpha, alpha_HW=alpha_HW) @@ -344,7 +392,8 @@ genParam = function(pop,simParam=NULL){ #' @title Additive variance #' -#' @description Returns additive variance for all traits +#' @description Returns additive variance (=variance of breeding values) +#' for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -370,7 +419,8 @@ varA = function(pop,simParam=NULL){ #' @title Dominance variance #' -#' @description Returns dominance variance for all traits +#' @description Returns dominance variance (=variance of dominance deviations) +#' for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -394,10 +444,38 @@ varD = function(pop,simParam=NULL){ genParam(pop,simParam=simParam)$varD } +#' @title Imprinting variance +#' +#' @description Returns imprinting variance (=variance of imprinting deviations) +#' for all traits +#' +#' @param pop an object of \code{\link{Pop-class}} +#' @param simParam an object of \code{\link{SimParam}} +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitAI(10, meanID=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' varI(pop, simParam=SP) +#' +#' @export +varI = function(pop,simParam=NULL){ + genParam(pop,simParam=simParam)$varI +} + #' @title Additive-by-additive epistatic variance #' #' @description Returns additive-by-additive epistatic -#' variance for all traits +#' variance (=variance of additive-by-additive epistatic deviations) +#' for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -421,13 +499,19 @@ varAA = function(pop,simParam=NULL){ genParam(pop,simParam=simParam)$varAA } -#' @title Breeding value +#' @title Breeding values #' -#' @description Returns breeding values for all traits +#' @description Returns breeding values for all traits. #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} #' +#' @details +#' With imprinting the output doesn't consider imprinting effects +#' into breeding values of individuals. +#' If you must get sex-specific breeding values, +#' use \code{bvSS()}, \code{bvM()} and \code{bvP()} - see above. +#' #' @examples #' #Create founder haplotypes #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) @@ -447,6 +531,132 @@ bv = function(pop,simParam=NULL){ genParam(pop,simParam=simParam)$bv } +#' @title Sex-specific breeding values +#' +#' @description Returns breeding values for all traits having into account +#' imprinting effects. +#' +#' @param pop an object of \code{\link{Pop-class}} +#' @param simParam an object of \code{\link{SimParam}} +#' +#' @details +#' With imprinting the output depends on sex of individuals. +#' If \code{sex == "H"} (in many plants) then mean breeding value is returned, +#' which is mean of maternal and paternal breeding value. +#' If \code{sex == "F"} (female) then maternal breeding value is returned. +#' If \code{sex == "M"} (male) then paternal breeding value is returned. +#' If you must get maternal or paternal breeding value irrespective of the +#' sex of individuals, use \code{bvM()} and \code{bvP()} - see examples. +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' # --- Additive & Dominance effects --- +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitAD(10, meanDD=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' bv(pop, simParam=SP) +#' +#' # --- Additive & Imprinting effects - Bisexual individuals --- +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitAI(10, meanID=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' pop@sex +#' # Mean breeding values +#' # (individuals are bisexual) +#' (tmpBv = bvSS(pop, simParam=SP)) +#' +#' # Maternal breeding values +#' # (as if individuals are females or for their female sexual organ) +#' (tmpMBv = bvM(pop, simParam=SP)) +#' (tmpMId = tmpMBv - tmpBv) +#' idM(pop, simParam=SP) +#' +#' # Paternal breeding values +#' # (as if individuals are males or for their male sexual organ) +#' (tmpPBv = bvP(pop, simParam=SP)) +#' (tmpPId = tmpPBv - tmpBv) +#' idP(pop, simParam=SP) +#' +#' # Compare the values +#' data.frame(sex = pop@sex, bv = tmpBv[, 1], bvM = tmpMBv[, 1], bvP = tmpPBv[, 1], +#' idM = tmpMId[, 1], idP = tmpPId[, 1]) +#' +#' # --- Male and female individuals --- +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$setSexes("yes_sys") +#' SP$addTraitAI(10, meanID=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' pop@sex +#' # Individual breeding values +#' # (individuals are either females or males) +#' (tmpBv = bvSS(pop, simParam=SP)) +#' +#' # Maternal breeding values +#' # (as if individuals are females) +#' (tmpMBv = bvM(pop, simParam=SP)) +#' +#' # Paternal breeding values +#' # (as if individuals are males) +#' (tmpPBv = bvP(pop, simParam=SP)) +#' +#' # Compare the values +#' data.frame(sex = pop@sex, bv = tmpBv[, 1], bvM = tmpMBv[, 1], bvP = tmpPBv[, 1]) +#' +#' @export +bvSS = function(pop,simParam=NULL){ + + if(is.null(simParam)){ + simParam = get("SP",envir=.GlobalEnv) + } + + ret = bv(pop,simParam=simParam) + ret[,1] = rep(9,nInd(pop)) + impTest = sapply(simParam$traits, FUN = function(x) .hasSlot(x,"impEff")) + if(any(impTest)){ + retM = bvM(pop,simParam=simParam) + retP = bvP(pop,simParam=simParam) + sel = pop@sex %in% "F" + if(any(sel)){ + ret[sel, impTest] = retM[sel, impTest] # for females use mat breed val + } + sel = pop@sex %in% "M" + if(any(sel)){ + ret[sel, impTest] = retP[sel, impTest] # for males use pat breed val + } + } + return(ret) +} + +#' @describeIn id Maternal breeding values (with imprinting) +#' @export +bvM = function(pop,simParam=NULL){ + genParam(pop,simParam=simParam)$bvM +} + +#' @describeIn id Paternal breeding values (with imprinting) +#' @export +bvP = function(pop,simParam=NULL){ + genParam(pop,simParam=simParam)$bvP +} + #' @title Dominance deviations #' #' @description Returns dominance deviations for all traits @@ -473,6 +683,38 @@ dd = function(pop,simParam=NULL){ genParam(pop,simParam=simParam)$dd } +#' @title Imprinting deviations +#' +#' @description Returns imprinting deviations for all traits +#' +#' @param pop an object of \code{\link{Pop-class}} +#' @param simParam an object of \code{\link{SimParam}} +#' +#' @details +#' Imprinting was developed in AlphaSimR under an orthogonal model and imprinting +#' deviations becomes the difference between gv and (bv + dd). This differs from +#' imprinting deviation concept explained in O'Brien EK & Wolf JB (2019) even +#' both models are equivalents. +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitADI(10, meanDD=0.75, meanID=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' id(pop, simParam=SP) +#' +#' @export +id = function(pop,simParam=NULL){ + genParam(pop,simParam=simParam)$id +} + #' @title Additive-by-additive epistatic deviations #' #' @description Returns additive-by-additive epistatic @@ -502,7 +744,8 @@ aa = function(pop,simParam=NULL){ #' @title Additive genic variance #' -#' @description Returns additive genic variance for all traits +#' @description Returns additive genic variance (=sum of variances of breeding +#' values at individual loci) for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -528,7 +771,8 @@ genicVarA = function(pop,simParam=NULL){ #' @title Dominance genic variance #' -#' @description Returns dominance genic variance for all traits +#' @description Returns dominance genic variance (=sum of variances of dominance +#' deviations at individual loci) for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} @@ -552,6 +796,33 @@ genicVarD = function(pop,simParam=NULL){ genParam(pop,simParam=simParam)$genicVarD } +#' @title Imprinting genic variance +#' +#' @description Returns imprinting genic variance (=sum of variances of +#' imprinting deviations at individual loci) for all traits +#' +#' @param pop an object of \code{\link{Pop-class}} +#' @param simParam an object of \code{\link{SimParam}} +#' +#' @examples +#' #Create founder haplotypes +#' founderPop = quickHaplo(nInd=10, nChr=1, segSites=10) +#' +#' #Set simulation parameters +#' SP = SimParam$new(founderPop) +#' SP$addTraitAI(10, meanID=0.5) +#' SP$setVarE(h2=0.5) +#' \dontshow{SP$nThreads = 1L} +#' +#' #Create population +#' pop = newPop(founderPop, simParam=SP) +#' genicVarI(pop, simParam=SP) +#' +#' @export +genicVarI = function(pop,simParam=NULL){ + genParam(pop,simParam=simParam)$genicVarI +} + #' @title Additive-by-additive genic variance #' #' @description Returns additive-by-additive epistatic @@ -581,7 +852,8 @@ genicVarAA = function(pop,simParam=NULL){ #' @title Total genic variance #' -#' @description Returns total genic variance for all traits +#' @description Returns total genic variance (=sum of variances of genetic +#' values at individual loci) for all traits #' #' @param pop an object of \code{\link{Pop-class}} #' @param simParam an object of \code{\link{SimParam}} diff --git a/src/calcGenParam.cpp b/src/calcGenParam.cpp index 93b67f08..26423e53 100644 --- a/src/calcGenParam.cpp +++ b/src/calcGenParam.cpp @@ -279,6 +279,361 @@ Rcpp::List calcGenParamE(const Rcpp::S4& trait, } } +// Calculates genetic parameters for traits with imprinting +Rcpp::List calcGenParamS(const Rcpp::S4& trait, + const Rcpp::S4& pop, + int nThreads){ + //Information from pop + bool hasD = trait.hasSlot("domEff"); + arma::uword nInd = pop.slot("nInd"); + arma::uword ploidy = pop.slot("ploidy"); + double dP = double(ploidy); + //Information from trait + const arma::Col& lociPerChr = trait.slot("lociPerChr"); + arma::uvec lociLoc = trait.slot("lociLoc"); + arma::vec a = trait.slot("addEff"); + arma::vec d; + arma::vec s; + // TODO: Expand to polyploids + arma::vec x(ploidy+2); // Genotype dossage + x(0) = 0; + x(1) = 1; + x(2) = 1; + x(3) = 2; + arma::vec xa = (x-dP/2.0)*(2.0/dP); // -1, 0, 0, 1 for diploids + arma::vec xaE = (x-dP/2.0)*(2.0/dP); // -1, 0, 0, 1 for diploids + arma::vec xd = x%(dP-x)*(2.0/dP)*(2.0/dP); // 0, 1, 1, 0 for diploids + arma::vec xdE = x%(dP-x)*(2.0/dP)*(2.0/dP); // 0, 1, 1, 0 for diploids + // TODO expand to polyploids + arma::vec xs = xd; // 0, -1, 1, 0 for diploids + xs(1) = -xs(1); + arma::vec xsE = xd; // 0, -1, 1, 0 for diploids + xsE(1) = -xsE(1); + double intercept = trait.slot("intercept"); + arma::mat bvMat(nInd,nThreads,arma::fill::zeros); // "Breeding value" + arma::mat bvMatM(nInd,nThreads,arma::fill::zeros); // "Breeding value (maternal)" + arma::mat bvMatP(nInd,nThreads,arma::fill::zeros); // "Breeding value (paternal)" + arma::mat gv_t; // Total genetic value + arma::mat gv_a(nInd,nThreads,arma::fill::zeros); // Genetic value due to a + arma::vec genicA(nThreads,arma::fill::zeros); // No LD + arma::vec genicAM(nThreads,arma::fill::zeros); // No LD (maternal) + arma::vec genicAP(nThreads,arma::fill::zeros); // No LD (paternal) + arma::vec genicA2(nThreads,arma::fill::zeros); // No LD and HWE + arma::vec genicAM2(nThreads,arma::fill::zeros); // No LD and HWE (maternal) + arma::vec genicAP2(nThreads,arma::fill::zeros); // No LD and HWE (paternal) + arma::vec genicD(nThreads,arma::fill::zeros); // No LD + arma::vec genicD2(nThreads,arma::fill::zeros); // No LD and HWE + arma::vec genicS(nThreads,arma::fill::zeros); // No LD (genic imprinting devation variance is the same between sexes) + arma::vec genicS2(nThreads,arma::fill::zeros); // No LD and HWE (genic imprinting devation variance is the same between sexes) + arma::vec mu(nThreads,arma::fill::zeros); // Observed mean + arma::vec eMu(nThreads,arma::fill::zeros); // Expected mean with HWE + arma::mat ddMat, gv_d; // Dominance deviation and genetic value due to d + if(hasD){ + d = Rcpp::as(trait.slot("domEff")); + ddMat.set_size(nInd,nThreads); + ddMat.zeros(); + gv_d.set_size(nInd,nThreads); + gv_d.zeros(); + } + arma::mat sdMat, gv_s; // Imprinting deviation and genetic value due to s + s = Rcpp::as(trait.slot("impEff")); + sdMat.set_size(nInd,nThreads); + sdMat.zeros(); + gv_s.set_size(nInd,nThreads); + gv_s.zeros(); + + arma::vec alpha(a.n_elem); + arma::vec alphaHW(a.n_elem); + + arma::vec alphaM(a.n_elem); + arma::vec alphaMHW(a.n_elem); + + arma::vec alphaP(a.n_elem); + arma::vec alphaPHW(a.n_elem); + + arma::vec beta(a.n_elem); + arma::vec betaHW(a.n_elem); + + arma::vec gamma(a.n_elem); + arma::vec gammaHW(a.n_elem); + + arma::vec m(a.n_elem); + arma::vec mE(a.n_elem); + + arma::vec m_a(a.n_elem); + arma::vec m_aE(a.n_elem); + + arma::vec m_d(a.n_elem); + arma::vec m_dE(a.n_elem); + + arma::Mat genoMat = getGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + arma::Mat genoMatM = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + arma::Mat genoMatP = getPaternalGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) num_threads(nThreads) +#endif + for(arma::uword i=0; i getGvA2(const Rcpp::S4& trait, return output; } +// Calculates genetic values using parental origin +// Useful for imprinting or genomic predictions +// TODO: we should add GxE to this function too?! +arma::field getGvS(const Rcpp::S4& trait, + const Rcpp::S4& pop, + int nThreads){ + arma::field output; + bool hasD = trait.hasSlot("domEff"); + bool hasS = trait.hasSlot("impEff"); + arma::uword nInd = pop.slot("nInd"); + arma::uword ploidy = pop.slot("ploidy"); + double dP = double(ploidy); + const arma::Col& lociPerChr = trait.slot("lociPerChr"); + arma::uvec lociLoc = trait.slot("lociLoc"); + arma::vec a,d,s; + a = Rcpp::as(trait.slot("addEff")); + if(hasD){ + d = Rcpp::as(trait.slot("domEff")); + } + s = Rcpp::as(trait.slot("impEff")); + arma::mat gv(nInd,nThreads); + gv.fill(double(trait.slot("intercept"))/double(nThreads)); + output.set_size(1); + output(0).set_size(nInd); + // Half ploidy for xa + // TODO: Expand to polyploids + arma::vec x(ploidy+2); // Genotype dosage + x(0) = 0; + x(1) = 1; + x(2) = 1; + x(3) = 2; + arma::vec xa(ploidy+2); + xa = (x-dP/2.0)*(2.0/dP); // -1, 0, 0, 1 for diploids + arma::vec xd = x%(dP-x)*(2.0/dP)*(2.0/dP); // 0, 1, 1, 0 for diploids + arma::vec xs = xd; // 0, -1, 1, 0 for diploids + xs(1) = -xs(1); + + arma::Mat maternalGeno = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + arma::Mat paternalGeno = getPaternalGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) num_threads(nThreads) +#endif + for(arma::uword i=0; i getGvE(const Rcpp::S4& trait, const Rcpp::S4& pop, @@ -119,18 +192,18 @@ arma::field getGvE(const Rcpp::S4& trait, #endif for(arma::uword j=0; j getGv(const Rcpp::S4& trait, // Genomic prediction return getGvA2(trait, pop, nThreads); } + if(trait.hasSlot("impEff")){ + return getGvS(trait, pop, nThreads); + } if(trait.hasSlot("epiEff")){ return getGvE(trait, pop, nThreads); } @@ -293,6 +369,72 @@ arma::field getGvIndexStd(const Rcpp::S4& trait, return output; } +// Basic implementation of getGv for multiple traits with imprinting +arma::field getGvIndexS(const Rcpp::S4& trait, + arma::Mat& genoMat, + arma::Mat& maternalGeno, + arma::Mat& paternalGeno, + arma::Col qtlIndex, + arma::uword ploidy, + int nThreads){ + qtlIndex = qtlIndex-1; // R to C++ + arma::field output; + bool hasD = trait.hasSlot("domEff"); + bool hasGxe = trait.hasSlot("gxeEff"); + arma::uword nInd = genoMat.n_rows; + double dP = double(ploidy); + arma::vec a,d,s,g; + a = Rcpp::as(trait.slot("addEff")); + if(hasD){ + d = Rcpp::as(trait.slot("domEff")); + } + s = Rcpp::as(trait.slot("impEff")); + arma::mat gv(nInd,nThreads),gxe; + gv.fill(double(trait.slot("intercept"))/double(nThreads)); + output.set_size(1); + output(0).set_size(nInd); + + arma::vec x(ploidy+2); // Genotype dosage + x(0) = 0; + x(1) = 1; + x(2) = 1; + x(3) = 2; + arma::vec xa(ploidy+2); + xa = (x-dP/2.0)*(2.0/dP); // -1, 0, 0, 1 for diploids + arma::vec xd = x%(dP-x)*(2.0/dP)*(2.0/dP); // 0, 1, 1, 0 for diploids + arma::vec xs = xd; // 0, -1, 1, 0 for diploids + xs(1) = -xs(1); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) num_threads(nThreads) +#endif + for(arma::uword i=0; i getGvIndexE(const Rcpp::S4& trait, arma::Mat& genoMat, @@ -468,6 +610,19 @@ Rcpp::List getGvIndex(const Rcpp::S4& pop, }else if(trait.hasSlot("epiEff")){ // Trait with epistasis output = getGvIndexE(trait, genoMat, qtlIndex(i), ploidy, nThreads); + }else if(trait.hasSlot("impEff")){ + // Trait with imprinting + arma::Mat maternalGeno; + arma::Mat paternalGeno; + maternalGeno = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), + Rcpp::as >(activeQtl.slot("lociPerChr")), + Rcpp::as(activeQtl.slot("lociLoc")), + nThreads); + paternalGeno = getPaternalGeno(Rcpp::as > >(pop.slot("geno")), + Rcpp::as >(activeQtl.slot("lociPerChr")), + Rcpp::as(activeQtl.slot("lociLoc")), + nThreads); + output = getGvIndexS(trait, genoMat, maternalGeno, paternalGeno, qtlIndex(i), ploidy, nThreads); }else{ // All other traits output = getGvIndexStd(trait, genoMat, qtlIndex(i), ploidy, nThreads); diff --git a/tests/testthat/test-addTrait.R b/tests/testthat/test-addTrait.R index 9e6b1296..d141fd9c 100644 --- a/tests/testthat/test-addTrait.R +++ b/tests/testthat/test-addTrait.R @@ -68,6 +68,49 @@ test_that("addTraitADG",{ expect_equal(abs(SP$traits[[1]]@intercept),0,tolerance=1e-6) }) +#Population with 4 individuals, +#4 genotypes (0/0, 0/1, 1/0, 1/1), +#1 chromosome and 1 QTL and p=q=0.5 +founderPop = newMapPop(list(c(0)), + list(matrix(c(0,0,0,1,1,0,1,1), + nrow=8,ncol=1))) + +test_that("addTraitAI",{ + SP = SimParam$new(founderPop=founderPop) + SP$nThreads = 1L + SP$addTraitAI(nQtlPerChr=1,mean=0,var=1,meanID=1) + pop = newPop(founderPop,simParam=SP) + # Just one locus and complete alias between a and i, will give varG=2, so + # scaled effect will be sqrt(varG) + expect_equal(abs(SP$traits[[1]]@addEff),sqrt(2),tolerance=1e-6) + expect_equal(abs(SP$traits[[1]]@impEff),sqrt(2),tolerance=1e-6) + expect_equal(abs(SP$traits[[1]]@intercept),0,tolerance=1e-6) + expect_equal(SP$varA,1,tolerance=1e-6) + expect_equal(SP$varG,2,tolerance=1e-6) + ans = genParam(pop,simParam=SP) + expect_equal(unname(c(ans$varA)),1,tolerance=1e-6) + expect_equal(unname(c(ans$varI)),1,tolerance=1e-6) + expect_equal(unname(c(ans$varG)),2,tolerance=1e-6) + expect_equal(unname(ans$genicVarA),1,tolerance=1e-6) + expect_equal(unname(ans$genicVarI),1,tolerance=1e-6) + expect_equal(unname(ans$genicVarG),2,tolerance=1e-6) + gv_a = c(-SP$traits[[1]]@addEff,0,0,SP$traits[[1]]@addEff) + gv_i = c(0,SP$traits[[1]]@impEff,-SP$traits[[1]]@impEff,0) + expect_equal(unname(c(ans$gv_a)),gv_a,tolerance=1e-6) + expect_equal(unname(c(ans$gv_i)),gv_i,tolerance=1e-6) + expect_equal(unname(c(ans$gv)),gv_a+gv_i,tolerance=1e-6) + expect_equal(unname(c(ans$bv)),gv_a,tolerance=1e-6) + # Test logic of how bv() and id() interact with sex + expect_equal(bv(pop),ans$bv,tolerance=1e-6) + expect_equal(bvM(pop),ans$bvM,tolerance=1e-6) + expect_equal(bvP(pop),ans$bvP,tolerance=1e-6) + pop@sex[] = "M" + expect_equal(bvSS(pop),ans$bvP,tolerance=1e-6) + pop@sex[] = "F" + expect_equal(bvSS(pop),ans$bvM,tolerance=1e-6) +}) + + #Population with 2 individuals, 1 chromosome and 2 QTL #Population is fully inbred and p=q=0.5 founderPop = newMapPop(list(c(0,0)), diff --git a/vignettes/TestAI.Rmd b/vignettes/TestAI.Rmd new file mode 100644 index 00000000..88bc4566 --- /dev/null +++ b/vignettes/TestAI.Rmd @@ -0,0 +1,195 @@ +--- +title: "Imprinting implementation checking" +author: "López-Carbonell, D." +date: "2023-12-01" +output: html_document +--- + +```{r setup, include=FALSE} +devtools::load_all() +``` + +## TraitAI testing + +This is an example to test an Additive and Imprinting effect. + +Firstly we create a founder Population + +```{r} +founderPop = quickHaplo(nInd=10, nChr=1, segSites=1) +``` + +SP object will be created from this founder population and the new addTraitAi function will be used. + +```{r pressure, echo=FALSE} +SP = SimParam$new(founderPop) +SP$addTraitAI(1, meanID=0.5) +``` + +The output are the different parameters and values that carry out the simulation. + +Also newPop function has been modified to create imprinted values if SP object has this info. We are going to use founderPop as a new pop. Then, values obtained in SP\$addTraitAI will be obtained (but without intercept). + +```{r pressure, } +pop = newPop(founderPop) +pop@gv +pop@gv - SP$traits[[1]]@intercept +``` + +From SP we can obtain the imprinting model values. Remember we are not including dominance. + +```{r pressure, } +addEff <- SP$traits[[1]]@addEff +impEff <- SP$traits[[1]]@impEff +intercept <- SP$traits[[1]]@intercept + +while (impEff < 0 || addEff < 0) { + SP = SimParam$new(founderPop) + SP$addTraitAI(1, meanID = 0.5) + pop = newPop(founderPop) + addEff <- SP$traits[[1]]@addEff + impEff <- SP$traits[[1]]@impEff + intercept <- SP$traits[[1]]@intercept +} +``` + +For calculating genomic values following an imprinted model Maternal and Paternal haplotypes will be required. Genomic data can also be obtained from haplotypes. + +```{r pressure, } +MatHaplo <- pullQtlHaplo(pop = pop, trait = 1, haplo = 1) +PatHaplo <- pullQtlHaplo(pop = pop, trait = 1, haplo = 2) +Geno <- pullQtlGeno(pop = pop, trait = 1) #it is also MatHaplo + PatHaplo +``` + +Imprinting is only affecting to heterocygous, then they have to be detected and homozygous dosage should be removed in haplotypic data. + +```{r pressure, } +Homozygous <- Geno %%2 == 0 + +MatHaploHet <- MatHaplo +MatHaploHet[Homozygous] <- 0 + +PatHaploHet <- PatHaplo +PatHaploHet[Homozygous] <- 0 +``` + +GV = mu + a + maternal imprinting + paternal imprinting Geno have to be codified as -1, 0, 1 + +```{r pressure, } +(calculated_gv <- intercept + (Geno -1) %*% addEff + MatHaploHet %*% -impEff + PatHaploHet %*% impEff) +``` + +This plot shows the Genomic values using newPop function and our results are the same. + +```{r non.finished.plotting, eval=FALSE} +plot(calculated_gv, pop@gv, abline(a = 0, b = 1)) + +``` +Alternatively, this can be also done using: +```{r} +calculated_gv_alt <- intercept + ((MatHaplo + PatHaplo)-1) * addEff + (PatHaplo - MatHaplo) * impEff + +``` + +This plot shows the Genomic values using newPop function and our alternative results are the same. + +```{r non.finished.plotting, eval=FALSE} +plot(calculated_gv_alt, pop@gv, abline(a = 0, b = 1)) + +``` + +We are going to prepare a graphical representation of the Imprinting Model, following . To do that, Genotypic Values will be stored in `table_gvs`. Furthermore, genotypes will be presented in 00, 01, 10 & 11 and A2A2, A1A2, A2A1 & A1A1 notation. + +```{r} +Geno_cod <- paste(MatHaplo, PatHaplo, sep = "") #codifiyng genotypes as 00, 01, 10, 11 +gvs = as.data.frame(cbind(gv(pop), Geno_cod)) +colnames(gvs) <- c("GV", "Genotype") +table_gvs <- unique(gvs) +table_gvs$GV <- as.numeric(table_gvs$GV) +(table_gvs <- table_gvs[order(table_gvs$GV),]) +table_gvs$Genotype_2 <- c("A2A2", "A1A2", "A2A1", "A1A1") +table_gvs$Value <- c("-a", "d-i", "d+i", "a") +``` + +Once values storaged, the plot is represented. + +```{r} +plot(x = table_gvs$GV, y = rep(1,4), type = "l", ylab = "", xlab = "Genotypic value", yaxt = "n") + + points(x = table_gvs$GV,y = rep(1,4), col = "salmon", pch = 18, cex = 2.5) + + text(x = table_gvs$GV, y = rep(1.2,4), labels = table_gvs$Genotype, cex = 1.5) + + text(x = table_gvs$GV, y = rep(0.85,4), labels = table_gvs$Genotype_2, cex = 1) + + text(x = table_gvs$GV, y = rep(0.7,4), labels = table_gvs$Value, cex = 1) +``` + +Also Falconer's Fig. 7.2. is recreated including under Imprinting model. To do that, Genotypic Values and Breeding Values are stored. + +```{r} +bvs = as.data.frame(cbind(gv(pop), bv(pop), bvM(pop), bvP(pop))) +colnames(bvs) <- c("GV","BV","BVM","BVP") +table_bvs <- unique(bvs) +table_bvs$GV <- as.numeric(table_bvs$GV) +table_bvs$BV <- as.numeric(table_bvs$BV) +table_bvs$BVM <- as.numeric(table_bvs$BVM) +table_bvs$BVP <- as.numeric(table_bvs$BVP) +(table_bvs <- table_bvs[order(table_bvs$BV), ]) +``` + +```{r} +y_range<- range(c(table_bvs$GV,table_bvs$BV,table_bvs$BVM,table_bvs$BVP)) + +plot(x = c(1:3), y = unique(table_bvs$BV), type = "l", xlab = "Genotypic value", lwd = 2.5, col= "darkolivegreen3", ylim = y_range) + points(x = c(1:3), y = unique(table_bvs$BV), col = "darkolivegreen3", pch = 18, cex = 2.5) + lines(x = c(1:3),y = unique(table_bvs$BVM), col = "salmon", lwd = 2.5) + points(x = c(1:3), y = unique(table_bvs$BVM), col = "salmon", pch = 18, cex = 2.5) + lines(x = c(1:3),y = unique(table_bvs$BVP), col = "cornflowerblue", lwd = 2.5) + points(x = c(1:3), y = unique(table_bvs$BVP), col = "cornflowerblue", pch = 18, cex = 2.5) + lines(x = c(1,1),y = c(unique(table_bvs$BV)[1],unique(table_bvs$BVP)[1]), col = "deeppink", lwd = 2.5) + lines(x = c(2,2),y = c(unique(table_bvs$BV)[2],unique(table_bvs$BVP)[2]), col = "deeppink", lwd = 2.5) + lines(x = c(3,3),y = c(unique(table_bvs$BV)[3],unique(table_bvs$BVP)[3]), col = "deeppink", lwd = 2.5) + lines(x = c(1,2),y = c(unique(table_bvs$BV)[1],unique(table_bvs$BV)[1]), col = "cyan4", lwd = 2.5, lty = 5) + lines(x = c(2,2),y = c(unique(table_bvs$BV)[1],unique(table_bvs$BV)[2]), col = "cyan4", lwd = 2.5, lty = 5) + lines(x = c(3,2.01),y = c(unique(table_bvs$BVM)[3],unique(table_bvs$BVM)[3]), col = "darkgoldenrod", lwd = 2.5, lty = 5) + lines(x = c(2.01,2.01),y = c(unique(table_bvs$BVM)[2],unique(table_bvs$BVM)[3]), col = "darkgoldenrod", lwd = 2.5, lty = 5) + lines(x = c(3,1.99),y = c(unique(table_bvs$BVP)[3],unique(table_bvs$BVP)[3]), col = "aquamarine3", lwd = 2.5, lty = 5) + lines(x = c(1.99,1.99),y = c(unique(table_bvs$BVP)[2],unique(table_bvs$BVP)[3]), col = "aquamarine3", lwd = 2.5, lty = 5) + legend("topleft", legend=c("BV_all", "BV_M", "BV_P", "Imp dev", + expression(paste(alpha,"_all",sep = "")), + expression(paste(alpha,"_M",sep = "")), + expression(paste(alpha,"_P",sep = ""))), + fill = c("darkolivegreen3","salmon","cornflowerblue", + "deeppink","cyan4","darkgoldenrod","aquamarine3"), + bty="n" +) + +``` + +Now we are going to repeat this procedure with a bigger example and it will work properly again. + +```{r} +founderPop = quickHaplo(nInd = 10000, + nChr = 1, + segSites = 200) +SP = SimParam$new(founderPop) +SP$addTraitAI(200, meanID = 0.5) +pop = newPop(founderPop) +addEff <- SP$traits[[1]]@addEff +impEff <- SP$traits[[1]]@impEff +intercept <- SP$traits[[1]]@intercept +MatHaplo <- pullQtlHaplo(pop = pop, trait = 1, haplo = 1) +PatHaplo <- pullQtlHaplo(pop = pop, trait = 1, haplo = 2) +Geno <- pullQtlGeno(pop = pop, trait = 1) + +Homozygous <- Geno %% 2 == 0 + +MatHaploHet <- MatHaplo +MatHaploHet[Homozygous] <- 0 + +PatHaploHet <- PatHaplo +PatHaploHet[Homozygous] <- 0 + +calculated_gv <- + intercept + (Geno - 1) %*% addEff + + MatHaploHet %*% -impEff + PatHaploHet %*% impEff + +plot(calculated_gv, pop@gv, abline(a = 0, b = 1)) + +```