Merge pull request #93 from gaynorr/devel

Version 1.3.2

.gitignore

@@ -22,3 +22,4 @@ vignettes/*.pdf
 Notes.txt
 .Rproj.user/
 docs
+.DS_Store

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: AlphaSimR
 Type: Package
 Title: Breeding Program Simulations
-Version: 1.3.1
-Date: 2022-8-25
+Version: 1.3.2
+Date: 2022-11-2
 Authors@R: c(person("Chris", "Gaynor", email = "[email protected]",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")),
   person("Gregor", "Gorjanc", role = "ctb", 

NEWS.md

@@ -1,3 +1,15 @@
+# AlphaSimR 1.3.2
+
+*fixed column name bug with multiple traits in `setEBV`
+
+*fixed CTD caused by `runMacs` when too many segSites are requested
+
+*fixed missing names in GV when using `resetPop`
+
+*fixed bug in `importTrait`
+
+*`popVar` now deals with matrices having 1 row
+
 # AlphaSimR 1.3.1
 
 *updated link to Gaynor, 2017
@@ -49,10 +61,6 @@
 *added `pullMarkerGeno` and `pullMarkerHaplo` to make it easier to extract genotypes and haplotypes of specific loci without defining a trait or SNP chip
 
 *`reduceGenome`, `mergeGenome` and `doubleGenome` should really now work with pedigree and recombination tracking
-
-## Known issues
-
-*`pedigreeCross` fails without an appropriate warning for some incomplete pedigrees
 
 # AlphaSimR 1.1.2
 

R/Class-Pop.R

@@ -746,6 +746,11 @@ resetPop = function(pop,simParam=NULL){
     simParam = get("SP",envir=.GlobalEnv)
   }
   pop@nTraits = simParam$nTraits
+
+  # Extract names to add back at the end
+  traitNames = colnames(pop@gv)
+
+  # Create empty slots for traits
   pop@pheno = matrix(NA_real_,
                      nrow=pop@nInd,
                      ncol=simParam$nTraits)
@@ -756,6 +761,8 @@ resetPop = function(pop,simParam=NULL){
   pop@gv = matrix(NA_real_,nrow=pop@nInd,
                   ncol=simParam$nTraits)
   pop@fixEff = rep(1L,pop@nInd)
+
+  # Calculate genetic values
   if(simParam$nTraits>=1){
     for(i in 1:simParam$nTraits){
       tmp = getGv(simParam$traits[[i]],pop,simParam$nThreads)
@@ -765,6 +772,10 @@ resetPop = function(pop,simParam=NULL){
       }
     }
   }
+
+  # Add back trait names
+  colnames(pop@pheno) = colnames(pop@gv) = traitNames
+
   return(pop)
 }
 

R/Class-SimParam.R

@@ -1247,7 +1247,7 @@ SimParam = R6Class(
         useDom = TRUE
         domEff = as.matrix(domEff)
         stopifnot(nrow(addEff)==nrow(domEff),
-                  ncol(addEff)==nrow(domEff))
+                  ncol(addEff)==ncol(domEff))
       }
 
       # Prepare the intercept

R/GS.R

@@ -1266,6 +1266,9 @@ setEBV = function(pop, solution, value="gv", targetPop=NULL,
                nrow=pop@nInd,
                ncol=nTraits)
 
+  # Placeholder names
+  colnames(ebv) = as.character(1:nTraits)
+
   value = tolower(value)
 
   if(value=="gv"){

R/founderPop.R

@@ -218,22 +218,34 @@ runMacs = function(nInd,nChr=1, segSites=NULL, inbred=FALSE, species="GENERIC",
     }
     command = paste0(popSize," ",speciesParams,splitI," ",speciesHist,splitJ," -s ")
   }
+
   if(!is.null(manualGenLen)){
     genLen = manualGenLen
   }
+
   if(length(genLen)==1){
     genLen = rep(genLen, nChr)
   }
+
+  # Run MaCS
   macsOut = MaCS(command, segSites, inbred, ploidy, 
                  nThreads, seed)
   dim(macsOut$geno) = NULL # Account for matrix bug in RcppArmadillo
+
+  # Check if desired number of loci were obtained
   nLoci = sapply(macsOut$genMap,length)
+  isLimited = segSites>0
+  if(any(nLoci[isLimited] != segSites[isLimited])){
+    stop("MaCS did not return enough segSites, use segSites=NULL to return all sites generated by MaCS")
+  }
+
   genMap = vector("list",nChr)
   for(i in 1:nChr){
     genMap[[i]] = genLen[i]*c(macsOut$genMap[[i]]-macsOut$genMap[[i]][1])
     names(genMap[[i]]) = paste(i,1:length(genMap[[i]]),sep="_")
   }
   names(genMap) = as.character(1:nChr)
+
   output = new("MapPop",
                nInd=nInd,
                nChr=nChr,

R/misc.R

@@ -1,3 +1,9 @@
+# Converts a matrix to integer type
+# Intended for genotype matrices of raw type
+convToImat = function(X){
+  return(matrix(as.integer(X), nrow=nrow(X) ,ncol=ncol(X)))
+}
+
 #' @rdname isFemale
 #' @title Test if individuals of a population are female or male
 #'

R/pullGeno.R

@@ -1,7 +1,3 @@
-convToImat = function(X){
-  return(matrix(as.integer(X), nrow=nrow(X) ,ncol=ncol(X)))
-}
-
 selectLoci = function(chr, inLociPerChr, inLociLoc){
   if(is.null(chr)){
     return(list(lociPerChr=inLociPerChr,

src/misc.cpp

@@ -25,7 +25,11 @@ unsigned char toByte(std::bitset<8> bits){
 //' @export
 // [[Rcpp::export]]
 arma::mat popVar(const arma::mat& X) {
-  return arma::cov(X,1);
+  if(X.n_rows==1){
+    return(arma::mat(X.n_cols,X.n_cols,arma::fill::zeros));
+  }else{
+    return arma::cov(X,1);
+  }
 }
 
 // Merges geno objects, i.e. fields containing cubes of unsigned char

src/simulator.cpp

@@ -699,7 +699,7 @@ Rcpp::List MaCS(Rcpp::String args, arma::uvec maxSites, bool inbred,
   //Check input
   string t = args;
   if (t == "") {
-    Rcpp::stop("error");
+    Rcpp::stop("error passing argument string");
   }
 
   // Output objects
@@ -715,9 +715,6 @@ Rcpp::List MaCS(Rcpp::String args, arma::uvec maxSites, bool inbred,
     // Run MaCS and check for valid output
     vector<AlphaSimRReturn> macsOutput;
     macsOutput = runFromAlphaSimR(args+seed(chr));
-    if(macsOutput.empty()){
-      Rcpp::stop("Macs has failed to run.");
-    }
 
     arma::uword nSites, nBins, nHap, nInd;
     nSites = macsOutput.size();
@@ -730,7 +727,7 @@ Rcpp::List MaCS(Rcpp::String args, arma::uvec maxSites, bool inbred,
     arma::uvec selSites;
     if(maxSites(chr)>0){
       if(nSites<maxSites(chr)){
-        Rcpp::stop("Not enough segregating sites generated");
+        maxSites(chr) = nSites;
       }
       selSites = sampleInt(maxSites(chr),nSites);
       nSites = maxSites(chr);