From 78ac6f751e1e0caa75601d478cbf9b6cc8d3687a Mon Sep 17 00:00:00 2001 From: david20011999 Date: Mon, 4 Dec 2023 16:58:16 +0100 Subject: [PATCH 1/2] add imprinting vignette --- vignettes/TestAI.Rmd | 183 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100644 vignettes/TestAI.Rmd diff --git a/vignettes/TestAI.Rmd b/vignettes/TestAI.Rmd new file mode 100644 index 00000000..03bdd264 --- /dev/null +++ b/vignettes/TestAI.Rmd @@ -0,0 +1,183 @@ +--- +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)) + +``` + +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)) + +``` From 154dc0a1af13227d9a776435b04a97b907d1316d Mon Sep 17 00:00:00 2001 From: david20011999 Date: Mon, 29 Jan 2024 14:23:41 +0100 Subject: [PATCH 2/2] add imprinting following Chris regression --- R/Class-LociMap.R | 31 ++ R/Class-SimParam.R | 166 ++++++++++- R/importData.R | 2 +- R/popSummary.R | 104 ++----- src/calcGenParam.cpp | 512 +++++++++++++++++++++------------ src/getGv.cpp | 108 +++++-- tests/testthat/test-addTrait.R | 9 - 7 files changed, 611 insertions(+), 321 deletions(-) diff --git a/R/Class-LociMap.R b/R/Class-LociMap.R index 66245168..d4ffea72 100644 --- a/R/Class-LociMap.R +++ b/R/Class-LociMap.R @@ -201,6 +201,37 @@ isTraitAI = function(x) { return(ret) } +#TraitADI---- +#' @title Additive dominance and imprinting trait +#' +#' @description Extends \code{\link{TraitAD-class}} +#' to add imprinting +#' +#' @slot impEff imprinting effects +#' +#' @export +setClass("TraitADI", + slots=c(impEff="numeric"), + contains="TraitAD") + +setValidity("TraitADI",function(object){ + errors = character() + if(object@nLoci!=length(object@impEff)){ + errors = c(errors,"nLoci!=length(impEff)") + } + if(length(errors)==0){ + return(TRUE) + }else{ + return(errors) + } +}) + +# Test if object is of a TraitAI class +isTraitADI = function(x) { + ret = is(x, class2 = "TraitADI") + return(ret) +} + #TraitA2D---- #' @title Sex specific additive and dominance trait #' diff --git a/R/Class-SimParam.R b/R/Class-SimParam.R index 7a5ef8d2..9fe52885 100644 --- a/R/Class-SimParam.R +++ b/R/Class-SimParam.R @@ -671,10 +671,13 @@ SimParam = R6Class( name = paste0("Trait",1:nTraits+self$nTraits) } stopifnot(length(mean)==length(var), + length(meanID)==length(mean), isSymmetric(corA), isSymmetric(corID), - length(mean)==nrow(corA), - length(mean)==length(name)) + nrow(corA)==nTraits, + nrow(corID)==nTraits, + length(varID)==nTraits, + length(name)==nTraits) qtlLoci = private$.pickLoci(nQtlPerChr) addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, corr=corA,gamma=gamma,shape=shape) @@ -689,7 +692,6 @@ SimParam = R6Class( name=name[i]) tmp = calcGenParam(trait, self$founderPop, self$nThreads) - if(useVarA){ scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1]) }else{ @@ -707,6 +709,111 @@ SimParam = R6Class( invisible(self) }, + #' @description + #' Randomly assigns eligible QTLs for one or more traits with dominance and imprinting (silencing) + #' If simulating more than one trait, all traits will be pleiotropic + #' with correlated effects. + #' + #' @param nQtlPerChr number of QTLs per chromosome. Can be a single value or nChr values. + #' @param mean a vector of desired mean genetic values for one or more traits + #' @param var a vector of desired genetic variances for one or more traits + #' @param meanDD mean dominance degree + #' @param varDD variance of dominance degree + #' @param meanID mean imprinting degree + #' @param varID variance of imprinting degree + #' @param corA a matrix of correlations between additive effects + #' @param corDD a matrix of correlations between dominance degrees + #' @param corID a matrix of correlations between imprinting degrees + #' @param useVarA tune according to additive genetic variance if true. If + #' FALSE, tuning is performed according to total genetic variance. + #' @param gamma should a gamma distribution be used instead of normal + #' @param shape the shape parameter for the gamma distribution + #' (the rate/scale parameter of the gamma distribution is accounted + #' for via the desired level of genetic variance, the var argument) + #' @param force should the check for a running simulation be + #' ignored. Only set to TRUE if you know what you are doing. + #' @param name optional name for trait(s) + #' + #' @examples + #' #Create founder haplotypes + #' founderPop = quickHaplo(nInd=10, nChr=1, segSites=1) + #' + #' #Set simulation parameters + #' SP = SimParam$new(founderPop) + #' \dontshow{SP$nThreads = 1L} + #' SP$addTraitADI(1, meanDD=0.5, meanID=0.5) + addTraitADI = function(nQtlPerChr,mean=0,var=1,meanDD=0,varDD=0,meanID=0, + varID=0,corA=NULL,corDD=NULL,corID=NULL,useVarA=TRUE, + gamma=FALSE,shape=1,force=FALSE,name=NULL){ + if(self$founderPop@ploidy > 2){ + stop("ERROR: Imprinting is not yet implemented for polyploids!") + } + if(!force){ + private$.isRunning() + } + if(length(nQtlPerChr)==1){ + nQtlPerChr = rep(nQtlPerChr,self$nChr) + } + nTraits = length(mean) + if(length(meanDD)==1) meanDD = rep(meanDD,nTraits) + if(length(meanID)==1) meanID = rep(meanID,nTraits) + if(length(varDD)==1) varDD = rep(varDD,nTraits) + if(length(varID)==1) varID = rep(varID,nTraits) + if(length(gamma)==1) gamma = rep(gamma,nTraits) + if(length(shape)==1) shape = rep(shape,nTraits) + if(is.null(corA)) corA=diag(nTraits) + if(is.null(corDD)) corDD=diag(nTraits) + if(is.null(corID)) corID=diag(nTraits) + if(is.null(name)){ + name = paste0("Trait",1:nTraits+self$nTraits) + } + stopifnot(length(mean)==length(var), + length(meanDD)==length(mean), + length(meanID)==length(mean), + isSymmetric(corA), + isSymmetric(corDD), + isSymmetric(corID), + nrow(corA)==nTraits, + nrow(corDD)==nTraits, + nrow(corID)==nTraits, + length(varDD)==nTraits, + length(varID)==nTraits, + length(name)==nTraits) + qtlLoci = private$.pickLoci(nQtlPerChr) + addEff = sampAddEff(qtlLoci=qtlLoci,nTraits=nTraits, + corr=corA,gamma=gamma,shape=shape) + domEff = sampDomEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, + corDD=corDD,meanDD=meanDD,varDD=varDD) + impEff = sampImpEff(qtlLoci=qtlLoci,nTraits=nTraits,addEff=addEff, + corID=corID,meanID=meanID,varID=varID) + for(i in 1:nTraits){ + trait = new("TraitADI", + qtlLoci, + addEff=addEff[,i], + domEff=domEff[,i], + impEff=impEff[,i], + intercept=0, + name=name[i]) + tmp = calcGenParam(trait, self$founderPop, + self$nThreads) + if(useVarA){ + scale = sqrt(var[i])/sqrt(popVar(tmp$bv)[1]) + }else{ + scale = sqrt(var[i])/sqrt(popVar(tmp$gv)[1]) + } + trait@addEff = trait@addEff*scale + trait@domEff = trait@domEff*scale + trait@impEff = trait@impEff*scale + trait@intercept = mean[i]-mean(tmp$gv*scale) + if(useVarA){ + private$.addTrait(trait,var[i],popVar(tmp$gv*scale)[1]) + }else{ + private$.addTrait(trait,popVar(tmp$bv*scale)[1],var[i]) + } + } + invisible(self) + }, + #' @description #' An alternative method for adding a trait with additive and dominance effects #' to an AlphaSimR simulation. The function attempts to create a trait matching @@ -1527,13 +1634,14 @@ SimParam = R6Class( #' formatting the trait as a \code{\link{LociMap-class}}. #' The formatting is performed automatically for the user, #' with more user friendly data.frames or matrices taken as - #' inputs. This function only works for A and AD trait types. + #' inputs. This function only works for A, AD, AI and ADI trait types. #' #' @param markerNames a vector of names for the QTL #' @param addEff a matrix of additive effects (nLoci x nTraits). #' Alternatively, a vector of length nLoci can be supplied for #' a single trait. #' @param domEff optional dominance effects for each locus + #' @param impEff optional imprinting effects for each locus #' @param intercept optional intercepts for each trait #' @param name optional name(s) for the trait(s) #' @param varE default error variance for phenotype, optional @@ -1542,6 +1650,7 @@ SimParam = R6Class( importTrait = function(markerNames, addEff, domEff=NULL, + impEff=NULL, intercept=NULL, name=NULL, varE=NULL, @@ -1562,6 +1671,16 @@ SimParam = R6Class( stopifnot(nrow(addEff)==nrow(domEff), ncol(addEff)==ncol(domEff)) } + if(is.null(impEff)){ + useImp = FALSE + }else{ + useImp = TRUE + impEff = as.matrix(impEff) + stopifnot(nrow(addEff)==nrow(impEff), + ncol(addEff)==ncol(impEff), + nrow(domEff)==nrow(impEff), + ncol(domEff)==ncol(impEff)) + } # Prepare the intercept if(is.null(intercept)){ @@ -1593,16 +1712,16 @@ SimParam = R6Class( # Create trait variables lociPerChr = integer(self$nChr) lociLoc = vector("list", self$nChr) - addEffList = domEffList = vector("list", nTraits) + addEffList = domEffList = impEffList = vector("list", nTraits) for(i in 1:nTraits){ - addEffList[[i]] = domEffList[[i]] = vector("list", self$nChr) + addEffList[[i]] = domEffList[[i]] = impEffList[[i]] = vector("list", self$nChr) } # Loop through chromosomes for(i in 1:self$nChr){ # Working on trait 1 # Initialize variables - addEffList[[1]][[i]] = domEffList[[1]][[i]] = numeric() + addEffList[[1]][[i]] = domEffList[[1]][[i]] = impEffList[[1]][[i]] = numeric() lociLoc[[i]] = integer() # Find matches if they exist @@ -1615,17 +1734,23 @@ SimParam = R6Class( if(useDom){ domEffList[[1]][[i]] = domEff[na.omit(take),1] } + if(useImp){ + impEffList[[1]][[i]] = impEff[na.omit(take),1] + } } # Work on additional traits? if(nTraits>1){ for(j in 2:nTraits){ - addEffList[[j]][[i]] = domEffList[[j]][[i]] = numeric() + addEffList[[j]][[i]] = domEffList[[j]][[i]] = impEffList[[j]][[i]] = numeric() if(lociPerChr[i]>0L){ addEffList[[j]][[i]] = addEff[na.omit(take),j] if(useDom){ domEffList[[j]][[i]] = domEff[na.omit(take),j] } + if(useImp){ + impEffList[[j]][[i]] = impEff[na.omit(take),j] + } } } } @@ -1639,6 +1764,18 @@ SimParam = R6Class( addEff = unlist(addEffList[[i]]) if(useDom){ domEff = unlist(domEffList[[i]]) + if(useImp){ + impEff = unlist(impEffList[[i]]) + trait = new("TraitADI", + addEff=addEff, + domEff=domEff, + impEff=impEff, + intercept=intercept[i], + nLoci=nLoci, + lociPerChr=lociPerChr, + lociLoc=lociLoc, + name=name[i]) + }else{ trait = new("TraitAD", addEff=addEff, domEff=domEff, @@ -1647,7 +1784,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], @@ -1655,6 +1804,7 @@ SimParam = R6Class( lociPerChr=lociPerChr, lociLoc=lociLoc, name=name[i]) + } } # Add trait to simParam diff --git a/R/importData.R b/R/importData.R index 34056927..b41dc7bc 100644 --- a/R/importData.R +++ b/R/importData.R @@ -276,7 +276,7 @@ importHaplo = function(haplo, genMap, ploidy=2L, ped=NULL){ genMap[[i]] = genMap[[i]] - genMap[[i]]-genMap[[i]][1] take = na.omit(take) } - haplotypes[[i]] = haplo[,take] + haplotypes[[i]] = as.matrix(haplo[,take]) } founderPop = newMapPop(genMap=genMap, diff --git a/R/popSummary.R b/R/popSummary.R index 6033c58e..fddd530b 100644 --- a/R/popSummary.R +++ b/R/popSummary.R @@ -187,7 +187,7 @@ genParam = function(pop,simParam=NULL){ # Blank nInd x nTrait matrices gv = matrix(NA_real_, nrow=nInd, ncol=nTraits) colnames(gv) = traitNames - bv = bvM = bvP = dd = idM = aa = gv_a = gv_d = gv_i = 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) @@ -222,14 +222,14 @@ genParam = function(pop,simParam=NULL){ if(.hasSlot(trait,"impEff")){ genicVarI[i] = tmp$genicVarI2 covI_HW[i] = tmp$genicVarI-tmp$genicVarI2 - idM[,i] = tmp$idM # taking just mat dev since pat dev = - mat dev + 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 - idM[,i] = rep(0,pop@nInd) + 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) @@ -254,11 +254,11 @@ genParam = function(pop,simParam=NULL){ covIAA_L[i] = 0 } else { covAD_L[i] = popVar(cbind(bv[,i],dd[,i]))[1,2] - covAI_L[i] = popVar(cbind(bv[,i],idM[,i]))[1,2] - covDI_L[i] = popVar(cbind(dd[,i],idM[,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(idM[,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.? } @@ -272,7 +272,7 @@ genParam = function(pop,simParam=NULL){ rownames(varD) = colnames(varD) = traitNames # TODO: do we need varIM and varIP (they are the same in a "random" setting)? - varI = popVar(idM) + varI = popVar(id) rownames(varI) = colnames(varI) = traitNames varAA = popVar(aa) @@ -316,7 +316,7 @@ genParam = function(pop,simParam=NULL){ bvM=bvM, bvP=bvP, dd=dd, - idM=idM, + id=id, aa=aa, gv_mu=gv_mu, gv_a=gv_a, @@ -581,102 +581,34 @@ dd = function(pop,simParam=NULL){ #' @title Imprinting deviations #' -#' @description Returns imprinting deviations for all traits according to the -#' sex of individuals. +#' @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 -#' If \code{sex == "H"} (in many plants) then mean imprinting deviation is returned, -#' which is 0 be definition. -#' If \code{sex == "F"} (female) then maternal imprinting deviation is returned. -#' If \code{sex == "M"} (male) then paternal imprinting deviation is returned, -#' which is -(maternal imprinting deviation). -#' -#' If you must get maternal or paternal imprinting deviation irrespective of the -#' sex of individuals, use \code{idM()} and \code{idP()} - see examples. +#' 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) #' -#' # --- Bisexual individuals --- #' #Set simulation parameters #' SP = SimParam$new(founderPop) -#' SP$addTraitAI(10, meanID=0.5) +#' 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) -#' pop@sex -#' # Mean imprinting deviations -#' # (individuals are bisexual) -#' (tmp = id(pop, simParam=SP)) -#' -#' # Maternal imprinting deviations -#' # (as if individuals are females or for their female sexual organ) -#' (tmpM = idM(pop, simParam=SP)) -#' -#' # Paternal imprinting deviations -#' # (as if individuals are males or for their male sexual organ) -#' (tmpP = idP(pop, simParam=SP)) -#' -#' # Compare the values -#' data.frame(sex = pop@sex, id = tmp[, 1], idM = tmpM[, 1], idP = tmpP[, 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 imprinting deviations -#' # (individuals are either females or males) -#' (tmp = id(pop, simParam=SP)) -#' -#' # Maternal imprinting deviations -#' # (as if individuals are females) -#' (tmpM = idM(pop, simParam=SP)) -#' -#' # Paternal imprinting deviations -#' # (as if individuals are males) -#' (tmpP = idP(pop, simParam=SP)) -#' -#' # Compare the values -#' data.frame(sex = pop@sex, id = tmp[, 1], idM = tmpM[, 1], idP = tmpP[, 1]) -#' +#' id(pop, simParam=SP) +#' #' @export id = function(pop,simParam=NULL){ - ret = idM(pop,simParam=simParam) - sel = pop@sex %in% "H" - if(any(sel)){ - ret[sel, ] = 0 # for bisexual use mean imp dev = 0 - } - sel = pop@sex %in% "M" # for males use pat imp dev - if(any(sel)){ - ret[sel, ] = -ret[sel, ] # pat imp dev = - mat imp dev - } - return(ret) -} - -#' @describeIn id Maternal imprinting deviations -#' @export -idM = function(pop,simParam=NULL){ - genParam(pop,simParam=simParam)$idM -} - -#' @describeIn id Paternal imprinting deviations -#' @export -idP = function(pop,simParam=NULL){ - -genParam(pop,simParam=simParam)$idM + genParam(pop,simParam=simParam)$id } #' @title Additive-by-additive epistatic deviations diff --git a/src/calcGenParam.cpp b/src/calcGenParam.cpp index b01f000e..620e97ca 100644 --- a/src/calcGenParam.cpp +++ b/src/calcGenParam.cpp @@ -273,6 +273,321 @@ 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"); + bool hasS = trait.hasSlot("impEff"); // Imprinting (=silencing) but using s since i is an iterator + 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 x0(ploidy+2); // Genotype dosage + x0(0) = 0; + x0(1) = 1; + x0(2) = 1; + x0(3) = 2; + arma::vec xa0(ploidy+2); + // arma::vec xd(ploidy+2); + // arma::vec xs(ploidy+2); + for(arma::uword i=0; i(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::Mat genoMatM; + genoMatM = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), + lociPerChr, lociLoc, nThreads); + arma::Mat genoMatP; + 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(trait.slot("impEff")); - sdMatM.set_size(nInd,nThreads); - sdMatP.set_size(nInd,nThreads); - sdMatM.zeros(); - sdMatP.zeros(); - gv_s.set_size(nInd,nThreads); - gv_s.zeros(); - } arma::Mat genoMat = getGeno(Rcpp::as > >(pop.slot("geno")), lociPerChr, lociLoc, nThreads); - arma::Mat genoMatM; - if(hasS){ - genoMatM = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), - lociPerChr, lociLoc, nThreads); - } #ifdef _OPENMP #pragma omp parallel for schedule(static) num_threads(nThreads) @@ -362,44 +651,22 @@ for(arma::uword i=0; i0) - sEffP(genoMat(j,i)) * (1 - genoMatM(j,i)); - // +i for paternal het (=01 for diploids, maternal allele is silenced when i>0) - // no need for sdMat since it would be zero (by definition) - sdMatM(j,tid) += sdM(genoMat(j,i)); - sdMatP(j,tid) += sdP(genoMat(j,i)); - } } } - if(hasD && !hasS){ + if(hasD){ gv_t = gv_a + gv_d; return Rcpp::List::create(Rcpp::Named("gv")=sum(gv_t,1)+intercept, Rcpp::Named("bv")=sum(bvMat,1), @@ -533,72 +736,6 @@ for(arma::uword i=0; i names; - names[0] = "gv"; - result[0] = sum(gv_t,1)+intercept; - names[1] = "bv"; - result[1] = sum(bvMat,1); - names[2] = "bvM"; - result[2] = sum(bvMatM,1); - names[3] = "bvP"; - result[3] = sum(bvMatP,1); - names[4] = "idM"; - result[4] = sum(sdMatM,1); - names[5] = "idP"; - result[5] = sum(sdMatP,1); - names[6] = "genicVarA"; - result[6] = accu(genicA); - names[7] = "genicVarAM"; - result[7] = accu(genicAM); - names[8] = "genicVarAP"; - result[8] = accu(genicAP); - names[9] = "genicVarI"; - result[9] = accu(genicS); - names[10] = "genicVarA2"; - result[10] = accu(genicA2); - names[11] = "genicVarAM2"; - result[11] = accu(genicAM2); - names[12] = "genicVarAP2"; - result[12] = accu(genicAP2); - names[13] = "genicVarI2"; - result[13] = accu(genicS2); - names[14] = "mu"; - result[14] = accu(mu)+intercept; - names[15] = "mu_HWE"; - result[15] = accu(eMu)+intercept; - names[16] = "gv_a"; - result[16] = sum(gv_a,1); - names[17] = "gv_i"; - result[17] = sum(gv_s,1); - names[18] = "gv_mu"; - result[18] = intercept; - result.attr("names") = Rcpp::wrap(names); }else{ return Rcpp::List::create(Rcpp::Named("gv")=sum(gv_a,1)+intercept, Rcpp::Named("bv")=sum(bvMat,1), @@ -610,3 +747,4 @@ for(arma::uword i=0; i getGvA2(const Rcpp::S4& trait, 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); @@ -16,18 +15,10 @@ arma::field getGvA2(const Rcpp::S4& trait, arma::uvec lociLoc = trait.slot("lociLoc"); arma::vec a1,a2,d,s; a1 = Rcpp::as(trait.slot("addEff")); - if(trait.hasSlot("addEffMale")){ - a2 = Rcpp::as(trait.slot("addEffMale")); - }else{ - a1 = a1 / 2; // the a1 and a2 in addEffMale are for a full dosage, but not with impEff - a2 = a1; - } + a2 = Rcpp::as(trait.slot("addEffMale")); if(hasD){ d = Rcpp::as(trait.slot("domEff")); } - if(hasS){ - s = Rcpp::as(trait.slot("impEff")); - } arma::mat gv(nInd,nThreads); gv.fill(double(trait.slot("intercept"))/double(nThreads)); output.set_size(1); @@ -40,10 +31,6 @@ arma::field getGvA2(const Rcpp::S4& trait, arma::vec xd(ploidy+1); for(arma::uword i=0; i maternalGeno = getMaternalGeno(Rcpp::as > >(pop.slot("geno")), lociPerChr, lociLoc, nThreads); @@ -67,26 +54,84 @@ arma::field getGvA2(const Rcpp::S4& trait, if(hasD){ dEff = xd*d(i); } - if(hasS){ - sEffM = xsM*s(i); - sEffP = xsP*s(i); - } for(arma::uword j=0; j0) - sEffP(tmp) * double(tmpP); - // +i for paternal het (=01 for diploids, maternal allele is silenced when i>0) - } - } + } + } + output(0) = sum(gv,1); + 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 + arma::vec x(ploidy+2); // Genotype dosage + x(0) = 0; + x(1) = 1; + x(2) = 1; + x(3) = 2; + arma::vec xa(ploidy+2); + for(arma::uword i=0; i 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, arma::field getGv(const Rcpp::S4& trait, const Rcpp::S4& pop, int nThreads){ - if(trait.hasSlot("impEff") || trait.hasSlot("addEffMale")){ + if(trait.hasSlot("addEffMale")){ // Imprinting or 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); } diff --git a/tests/testthat/test-addTrait.R b/tests/testthat/test-addTrait.R index 745142c0..155da74d 100644 --- a/tests/testthat/test-addTrait.R +++ b/tests/testthat/test-addTrait.R @@ -94,27 +94,18 @@ test_that("addTraitAI",{ 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) - idM = c(+SP$traits[[1]]@impEff, 0, 0, -SP$traits[[1]]@impEff) - # +2pi, (p-q)i, (p-q)i, -2qi --> with p=q=0.5 --> +i, 0, 0, -i - idP = c(-SP$traits[[1]]@impEff, 0, 0, +SP$traits[[1]]@impEff) - # -2pi, (q-p)i, (q-p)i, +2qi --> with p=q=0.5 --> -i, 0, 0, +i 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) - expect_equal(unname(c(ans$idM)),idM,tolerance=1e-6) - expect_equal(unname(c(ans$bvM)),unname(c(ans$bv))+idM,tolerance=1e-6) - expect_equal(unname(c(ans$bvP)),unname(c(ans$bv))+idP,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(bv(pop),ans$bvP,tolerance=1e-6) - expect_equal(unname(c(id(pop))),idP,tolerance=1e-6) pop@sex[] = "F" expect_equal(bv(pop),ans$bvM,tolerance=1e-6) - expect_equal(unname(c(id(pop))),idM,tolerance=1e-6) }) #Population with 2 individuals, 1 chromosome and 2 QTL