Merge branch 'devel' of https://github.com/HighlanderLab/SIMplyBee into devel

Author: gregorgorjanc

.gitignore

@@ -40,3 +40,5 @@ vignettes/*.pdf
 
 # Other
 .DS_Store
+src/*.o
+src/*.so

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: SIMplyBee
 Type: Package
 Title: 'AlphaSimR' Extension for Simulating Honeybee Populations and Breeding Programmes
-Version: 0.4.1
+Version: 0.4.2
 Authors@R: c(
   person("Jana", "Obšteter", email = "obsteter.jana@gmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1511-3916")),
@@ -27,7 +27,7 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Imports: methods, R6, stats, utils, extraDistr (>= 1.9.1), RANN, Rcpp (>= 0.12.7)
-Depends: R (>= 3.3.0), AlphaSimR (>= 1.5.3)
+Depends: R (>= 3.3.0), AlphaSimR (>= 1.6.1)
 LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH
 RoxygenNote: 7.3.2
 Suggests: 

NAMESPACE

@@ -69,7 +69,6 @@ export(getGv)
 export(getIbdHaplo)
 export(getId)
 export(getLocation)
-export(getMisc)
 export(getPheno)
 export(getPooledGeno)
 export(getQtlGeno)
@@ -192,7 +191,6 @@ export(replaceWorkers)
 export(resetEvents)
 export(selectColonies)
 export(setLocation)
-export(setMisc)
 export(setQueensYearOfBirth)
 export(simulateHoneyBeeGenomes)
 export(split)

NEWS.md

@@ -4,6 +4,28 @@ editor_options:
     wrap: 72
 ---
 
+# SIMplyBee version 0.4.2
+
+-   20??-??-??
+
+## Major changes
+-   TODO
+
+## New features
+-   TODO
+
+## Bug fixes
+-   editCsdLocus() works now also with just 2 alleles #591
+
+## Backgrond/package/etc. work
+-   calcBeeGRMIbs() can now return centred genotype matrix, allele frequencies,
+and scale factor used to calculate the GRM #594
+-   Improved default for getPooledGeno() (to type="mean"") and added an example
+on how to collect pooled workers' genotype accross colonies (but pooling is
+done within a colony!) #592
+-   We now removed setMisc() and getMisc() because we now use the new AlphaSimR
+structure of the misc slot that is much easier to use. #584
+
 # SIMplyBee version 0.4.1
 
 -   2024-09-19
@@ -26,7 +48,7 @@ which caused an error. We now read in the locations from a csv file.
     now c(0, 0) PR#500
 
 
-## New features ##
+## New features
 -   In setLocation(MultiColony) we can set one location (numeric) or
     multiple (list or data.frame) PR#500
 -   getLocation(MultiColony) got the collapse argument
@@ -57,8 +79,6 @@ which caused an error. We now read in the locations from a csv file.
 -   Bug fix - get\*Haplo() functions were returning diploid drones when
     input was a Pop-class
 
--   
-
 # SIMplyBee version 0.3.0
 
 -   2022-12-05 First public/CRAN version of the package

R/Functions_L0_auxilary.R

@@ -1276,8 +1276,6 @@ getCasteId <- function(x, caste = "all", collapse = FALSE, simParamBee = NULL) {
 #'   vector with sex information
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getCaste}}
-#'
 #' @return when \code{x} is \code{\link[AlphaSimR]{Pop-class}} for \code{caste != "all"}
 #'  or list for \code{caste == "all"} with sex nodes named by caste;
 #'    when \code{x} is \code{\link[SIMplyBee]{Colony-class}} return is a named list of
@@ -2869,7 +2867,7 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) {
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getIbdHaplo}} and \code{\link[AlphaSimR]{pullIbdHaplo}}
+#' @seealso \code{\link[AlphaSimR]{pullIbdHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -3124,7 +3122,7 @@ getDronesIbdHaplo <- function(x, nInd = NULL, chr = NULL, snpChip = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getQtlHaplo}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
+#' @seealso \code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
 #'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
@@ -3374,7 +3372,9 @@ getDronesQtlHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlGeno}} as well as
+#' @seealso \code{\link[SIMplyBee]{getQtlHaplo}},
+#'   \code{\link[AlphaSimR]{pullQtlGeno}}, and
+#'   \code{\link[SIMplyBee]{getPooledGeno}}, as well as
 #'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
@@ -3624,7 +3624,7 @@ getDronesQtlGeno <- function(x, nInd = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSegSiteHaplo}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
+#' @seealso \code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -3866,7 +3866,9 @@ getDronesSegSiteHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteGeno}}
+#' @seealso \code{\link[SIMplyBee]{getSegSiteHaplo}},
+#'   \code{\link[AlphaSimR]{pullSegSiteGeno}}, and
+#'   \code{\link[SIMplyBee]{getPooledGeno}}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
 #'   list of matrices with genotypes when \code{x} is
@@ -4106,7 +4108,7 @@ getDronesSegSiteGeno <- function(x, nInd = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSnpHaplo}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
+#' @seealso \code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -4347,7 +4349,9 @@ getDronesSnpHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpGeno}}
+#' @seealso \code{\link[SIMplyBee]{getSnpHaplo}},
+#'  \code{\link[AlphaSimR]{pullSnpGeno}}, and
+#'  \code{\link[SIMplyBee]{getPooledGeno}}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
 #'   list of matrices with genotypes when \code{x} is
@@ -4573,14 +4577,14 @@ getDronesSnpGeno <- function(x, nInd = NULL,
 #'   genotypes to mimic genotyping of a pool of colony members.
 #'
 #' @param x matrix, true genotypes with individuals in rows and sites in columns
-#' @param type character, "mean" for average genotype or "count" for the counts
-#'   of reference and alternative alleles
+#' @param type character, \code{"mean"} for average genotype (default) or
+#'   \code{"count"} for the counts of reference and alternative alleles
 #' @param sex character, vector of "F" and "M" to denote the sex of individuals
 #'   in \code{x}
 #'
 #' @return a numeric vector with average allele dosage when \code{type = "mean"}
-#'   and a two-row matrix with the counts of reference (1st row) and
-#'   alternative (2nd row) alleles
+#'   or a two-row matrix with the counts of reference (1st row) and
+#'   alternative (2nd row) alleles when \code{type = "count"}
 #'
 #' @examples
 #' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 50)
@@ -4632,11 +4636,19 @@ getDronesSnpGeno <- function(x, nInd = NULL,
 #'
 #' # As an exercise you could repeat the above with different numbers of workers!
 #'
+#' # How to get pooled genotypes of workers across multiple colonies?
+#' tmp <- getSegSiteGeno(x = apiary, caste = "workers")
+#' (tmp2 = lapply(X = tmp, FUN = getPooledGeno)) # as a list of one row matrices
+#' t(sapply(X = tmp, FUN = getPooledGeno)) # as one matrix - option A
+#' do.call(what = rbind, args = tmp2) # as one matrix - option B
 #' @export
-getPooledGeno <- function(x, type = NULL, sex = NULL) {
+getPooledGeno <- function(x, type = "mean", sex = NULL) {
   if (!is.matrix(x)) {
     stop("Argument x must be a matrix class object!")
   }
+  if (is.null(type) | !(type %in% c("mean", "count"))) {
+    stop("Argument type must be specified as either mean or count!")
+  }
   n <- nrow(x)
   if (is.null(sex)) {
     warning("Argument sex is NULL. Assuming that all individuals are female/diploid!")
@@ -4654,8 +4666,6 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
   } else if (type == "count") {
     ret <- rbind(nPloids - ret, ret)
     rownames(ret) <- c("0", "1")
-  } else {
-    stop("Argument type must be mean or count!")
   }
   return(ret)
 }
@@ -4676,9 +4686,16 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' @param sex character vector denoting sex for individuals with genotypes in
 #'   \code{x} - \code{"F"} for female and \code{"M"} for male
 #' @param alleleFreq numeric, vector of allele frequencies for the sites in
-#'   \code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}} is used
+#'   \code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}}
+#'   is used
+#' @param returnComponents logical, return GRM as well as the components used
+#'   to compute it (useful for GWAS by GBLUP)
 #'
-#' @return matrix of genomic relatedness coefficients
+#' @return if \code{returnComponents = FALSE} (default) return a matrix of
+#'   genomic relatedness coefficients; if \code{returnComponents = TRUE} return
+#'   a list with the GRM, centred genotype matrix, allele frequencies, and
+#'   scaling factor used to scale the crossproduct of centred genotype matrix
+#'   to get the GRM.
 #'
 #' @references Druet and Legarra (2020) Theoretical and empirical comparisons of
 #'   expected and realized relationships for the X-chromosome. Genetics
@@ -4704,12 +4721,12 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' GRM <- calcBeeGRMIbs(x = geno, sex = sex)
 #' # You can visualise this matrix with the function image() from the package 'Matrix'
 #'
-#' #Look at the diagonal at the relationship matrix
+#' # Look at the diagonal at the relationship matrix
 #' x <- diag(GRM)
 #' hist(x)
 #' summary(x)
 #'
-#' #Look at the off-diagonal at the relationship matrix
+#' # Look at the off-diagonal at the relationship matrix
 #' x <- GRM[lower.tri(x = GRM, diag = FALSE)]
 #' hist(x)
 #' summary(x)
@@ -4745,8 +4762,12 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"))
 #' # You can now compare how this compare to relationships between the queen
 #' # individual workers!
+#'
+#' # Return components
+#' calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"),
+#'               returnComponents = TRUE)
 #' @export
-calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL) {
+calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL, returnComponents = FALSE) {
   if (!is.matrix(x)) {
     stop("Argument x must be a matrix class object!")
   }
@@ -4776,8 +4797,13 @@ calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL) {
     #       This would overwrite x only once, at expense of doubling RAM
     x[, site] <- x[, site] - ploidy * alleleFreq[site]
   }
-  G <- tcrossprod(x) / (2 * sum(alleleFreq * (1 - alleleFreq)))
-  return(G)
+  scale <- 2 * sum(alleleFreq * (1 - alleleFreq))
+  G <- tcrossprod(x) / scale
+  if (returnComponents) {
+    return(list(G = G, x = x, alleleFreq = alleleFreq, scale = scale))
+  } else {
+    return(G)
+  }
 }
 
 #' @describeIn calcBeeGRMIbs Calculate allele frequencies from honeybee genotypes
@@ -6322,6 +6348,16 @@ calcColonyAa <- function(x, FUN = mapCasteToColonyAa, simParamBee = NULL, ...) {
 #'   in \code{\link[SIMplyBee]{SimParamBee}}. The two csd alleles must be different to
 #'   ensure heterozygosity at the csd locus.
 #' @param simParamBee global simulation parameters.
+#' @examples
+#' founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 10)
+#'
+#' SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 2)
+#' tmp <- createVirginQueens(founderGenomes)
+#' getCsdAlleles(tmp)
+#'
+#' SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4)
+#' tmp <- createVirginQueens(founderGenomes)
+#' getCsdAlleles(tmp)
 #'
 #' @return Returns an object of \code{\link[AlphaSimR]{Pop-class}}
 editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) {
@@ -6334,7 +6370,10 @@ editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) {
     alleles <- expand.grid(as.data.frame(matrix(rep(0:1, length(csdSites)), nrow = 2, byrow = FALSE)))
     # Sample two different alleles (without replacement) for each individual
     nAlleles <- simParamBee$nCsdAlleles
-    alleles <- sapply(seq_len(pop@nInd), FUN = function(x) list(alleles[sample(nAlleles, size = 2, replace = F), ]))
+    alleles <- sapply(X = seq_len(pop@nInd),
+                      FUN = function(x) {
+                        list(alleles[sample(nAlleles, size = 2, replace = FALSE), , drop = FALSE])
+                      })
   }
 
   if (pop@nInd != length(alleles)) {

R/Functions_L1_Pop.R

@@ -423,10 +423,8 @@ createCastePop <- function(x, caste = NULL, nInd = NULL,
       stop("MapPop-class can only be used to create virgin queens!")
     }
     ret <- newPop(x, simParam = simParamBee)
-    if (!is.null(simParamBee$csdChr)) {
-      if (editCsd) {
-        ret <- editCsdLocus(ret, alleles = csdAlleles, simParamBee = simParamBee)
-      }
+    if (!is.null(simParamBee$csdChr) && editCsd) {
+      ret <- editCsdLocus(pop = ret, alleles = csdAlleles, simParamBee = simParamBee)
     }
     ret@sex[] <- "F"
     simParamBee$changeCaste(id = ret@id, caste = "virginQueens")