Adding latest changes:

- Improved functions with faster benchmarks - New function to read data (accomodates mappoly2 datasets) - Small changes in code to meet CRAN policies and guidelines - Minor bug fixes
gabrielgesteira · Mar 20, 2024 · 5c9ec33 · 5c9ec33
1 parent 47faed3
commit 5c9ec33
Show file tree

Hide file tree

Showing 14 changed files with 33 additions and 45 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -44,7 +44,7 @@ LazyDataCompression: xz
 Depends:
     R (>= 4.0)
 Imports:
-     ggplot2 (>= 3.1), abind (>= 1.4), MASS (>= 7.3), gtools (>= 3.9.2), CompQuadForm, Matrix, RLRsim, mvtnorm, nlme, quadprog, parallel, doParallel, foreach, stats, methods, Rcpp (>= 0.12.19)
+     ggplot2 (>= 3.1), abind (>= 1.4), MASS (>= 7.3), gtools (>= 3.9.2), CompQuadForm, Matrix, RLRsim, mvtnorm, nlme, quadprog, parallel, doParallel, foreach, stats, methods, mappoly, Rcpp (>= 0.12.19)
 LinkingTo: Rcpp, RcppArmadillo, RcppProgress
-Suggests: mappoly, rmarkdown, devtools, knitr
+Suggests: rmarkdown, devtools, knitr
 RoxygenNote: 7.2.3
diff --git a/NAMESPACE b/NAMESPACE
@@ -53,6 +53,7 @@ importFrom(abind,abind)
 importFrom(graphics,abline)
 importFrom(graphics,plot)
 importFrom(gtools,combinations)
+importFrom(mappoly,calc_genoprob)
 importFrom(mvtnorm,pmvnorm)
 importFrom(nlme,fixed.effects)
 importFrom(nlme,fixef)

diff --git a/R/fit_model.R b/R/fit_model.R
@@ -288,8 +288,6 @@ summary.qtlpoly.fitted <- function(object, pheno.col=NULL, ...) {
 #' Adapts genomic incidence and relationship (varcov) matrices to run using sommer's C++ core function (v. 4.0 or higher)
 #' Function adapted from sommer v. 3.6 (Author: Giovanny Covarrubias-Pazaran)
 #'
-#' @param void internal function to be documented
-#'
 #' @keywords internal
 #'
 #' @author Gabriel de Siqueira Gesteira, \email{gdesiqu@@ncsu.edu}

diff --git a/R/fit_model2.R b/R/fit_model2.R
@@ -14,13 +14,9 @@
 #'
 #' @param verbose if \code{TRUE} (default), current progress is shown; if
 #'     \code{FALSE}, no output is produced.
-#' 
-#' @param object an object of class \code{qtlpoly.fitted} to be summarized.
 #'
 #' @param pheno.col a numeric vector with the phenotype column numbers to be summarized; if \code{NULL} (default), all phenotypes from \code{'data'} will be included.
 #'
-#' @param ... currently ignored
-#' 
 #' @return An object of class \code{qtlpoly.fitted} which contains a list of \code{results} for each trait with the following components:
 #'
 #'     \item{pheno.col}{a phenotype column number.}

diff --git a/R/null_model2.R b/R/null_model2.R
@@ -58,6 +58,7 @@ null_model2 <- function(data, offset.data = NULL, pheno.col = NULL, n.clusters =
   results <- vector("list", length(pheno.col))
   names(results) <- colnames(data$pheno)[pheno.col]
   markers <- c(1:data$nmrk)
+  m = 1
 
   for(p in 1:length(results)) {
 

diff --git a/R/read_data2.R b/R/read_data2.R
@@ -18,12 +18,6 @@
 #'
 #' @param verbose if \code{TRUE} (default), current progress is shown; if \code{FALSE}, no output is produced.
 #'
-#' @param x an object of class \code{qtlpoly.data} to be printed.
-#'
-#' @param detailed if \code{TRUE}, detailed information on linkage groups and phenotypes in shown; if \code{FALSE}, no details are printed.
-#'
-#' @param ... currently ignored
-#'
 #' @return An object of class \code{qtlpoly.data} which is a list containing the following components:
 #'
 #'     \item{ploidy}{a scalar with ploidy level.}
@@ -62,10 +56,12 @@
 #' @export read_data2
 #' @importFrom abind abind
 #' @importFrom gtools combinations
+#' @importFrom mappoly calc_genoprob
 
 read_data2 <- function(ploidy = 6, geno.prob, geno.dose = NULL, double.reduction = FALSE, pheno, weights = NULL, step = 1, verbose = TRUE) {
 
-  if (class(geno.prob) == "mappoly2.sequence"){
+  if(inherits(geno.prob, "mappoly2.sequence")){
+  ## if (class(geno.prob) == "mappoly2.sequence"){
 
   if(is.null(step)) step <- 1e-10
 
@@ -78,11 +74,12 @@ read_data2 <- function(ploidy = 6, geno.prob, geno.dose = NULL, double.reduction
     probs = x$map.genome$phase[[1]]$haploprob
     a = split(1:nrow(probs), ceiling(seq_along(1:nrow(probs)) / (ploidy*2)))
     b = lapply(a, function(y) return(as.matrix(probs[y,-c(1:3)])))
-    c = abind::abind(b, along = 3)
+    c = abind(b, along = 3)
     dimnames(c)[[1]] = letters[1:(ploidy*2)]
     dimnames(c)[[2]] = rownames(x$map.genome$phase[[1]]$p1)
     dimnames(c)[[3]] = raw.individual.names
-    map = c(0, cumsum(mappoly::imf_h(x$map.genome$phase[[1]]$rf)))
+    mpgpt = calc_genoprob # to ensure mappoly's function is required in the package
+    map = c(0, cumsum(imf_h(x$map.genome$phase[[1]]$rf)))
     names(map) = rownames(x$map.genome$phase[[1]]$p1)
     return(list(probs = c, map = map))
   })
@@ -523,3 +520,9 @@ read_data2 <- function(ploidy = 6, geno.prob, geno.dose = NULL, double.reduction
  }
 
 }
+
+## Support function from mappoly
+imf_h <- function(r) {
+  r[r >= 0.5] <- 0.5 - 1e-14
+  -50 * log(1 - 2 * r)
+}
diff --git a/R/remim2.R b/R/remim2.R
@@ -111,6 +111,7 @@ remim2 <- function (data, pheno.col = NULL, w.size = 15, sig.fwd = 0.01,
       weight <- data$weights[ind, pheno.col[p]]
     } else weight <- rep(1, length(ind))
     markers <- c(1:data$nmrk)
+    m = 1
     G <- data$G[ind,ind,markers]
     for(i in markers){
       G[,,i] = G[,,i]/mean(diag(G[,,i]))

diff --git a/R/varComp.R b/R/varComp.R
@@ -2083,17 +2083,17 @@ varComp.test.varComp = function(object, object2,
 			alt=object
 		}
 		if(!all(nul$random.labels %in% alt$random.labels)) stop("Two model ares not nested.")
-		varComp.test.2modelDoTest(null.fit=nul, alt.fit=alt, test=test, control=control, ...)
+		varCompTest2modelDoTest(null.fit=nul, alt.fit=alt, test=test, control=control, ...)
 	}else if(!missing(additional.varcov)){  ## null fit is available
 		if(!missing(null)) warning("'null' is ignored when 'additional.varcov' is provided.")
-		varComp.test.nulDoTest(null.fit=object, additional.varcov = additional.varcov, test=test, control=control, ...)
+		varCompTestNulDoTest(null.fit=object, additional.varcov = additional.varcov, test=test, control=control, ...)
 	}else if(!missing(null)){ ##  alt fit is available
-		varComp.test.altDoTest(alt.fit=object, null=null, test=test,  control=control, ...)
+		varCompTestAltDoTest(alt.fit=object, null=null, test=test,  control=control, ...)
 	}else {  ## default testing all components of object
-		varComp.test.altDoTest(alt.fit=object, null=integer(0L), test=test,  control=control, ...)
+		varCompTestAltDoTest(alt.fit=object, null=integer(0L), test=test,  control=control, ...)
 	}
 }
-varComp.test.2modelDoTest = function(null.fit, alt.fit, test='LinScore', control=varCompTest.control(test), ...)
+varCompTest2modelDoTest = function(null.fit, alt.fit, test='LinScore', control=varCompTest.control(test), ...)
 {
   # information=match.arg(information, informationTypes)
   varCompScoreTests = c(score='LinScore',score='VM03',score='SS95')
@@ -2136,7 +2136,7 @@ varComp.test.2modelDoTest = function(null.fit, alt.fit, test='LinScore', control
 	environment(varComp.test.Common)  = sys.frame(sys.nframe())
 	varComp.test.Common()
 }
-varComp.test.nulDoTest = function(null.fit, additional.varcov, test='LinScore', control=varCompTest.control(test), alt.fit=NULL, ...)
+varCompTestNulDoTest = function(null.fit, additional.varcov, test='LinScore', control=varCompTest.control(test), alt.fit=NULL, ...)
 {
   # information=match.arg(information, informationTypes)
   varCompScoreTests = c(score='LinScore',score='VM03',score='SS95')
@@ -2199,7 +2199,7 @@ varComp.test.nulDoTest = function(null.fit, additional.varcov, test='LinScore',
 	varComp.test.Common()
 }
 
-varComp.test.altDoTest = function(alt.fit, null=integer(0L), test='LinScore', control=varCompTest.control(test), null.fit=NULL, ...)
+varCompTestAltDoTest = function(alt.fit, null=integer(0L), test='LinScore', control=varCompTest.control(test), null.fit=NULL, ...)
 {
   # information=match.arg(information, informationTypes)
   varCompScoreTests = c(score='LinScore',score='VM03',score='SS95')
@@ -2425,7 +2425,8 @@ function(control, infoMat, tau.idx, LIkLI, tr1, n, LIy, all.scores)# , ...)
     repeat{
       i=i+1
       tmp=try( solve.QP(Dmat=invNegHess/fact, dvec=drop(invNegHess%*%all.scores[tau.idx])/fact, Amat=Amat, bvec=bvec, meq=0L, factorized=FALSE)$value )
-      if(class(tmp)!='try-error'){
+      if(!inherits(tmp, "try-error")){
+      ## if(class(tmp)!='try-error'){
         qp.term=tmp*fact
         break
       }else if(i==20L){

diff --git a/README.md b/README.md
@@ -1,7 +1,7 @@
-[![Travis Build Status](https://app.travis-ci.com/gabrielgesteira/QTLpoly.svg?branch=main)](https://app.travis-ci.com/gabrielgesteira/QTLpoly)
+r[![Travis Build Status](https://app.travis-ci.com/gabrielgesteira/QTLpoly.svg?branch=main)](https://app.travis-ci.com/gabrielgesteira/QTLpoly)
  [![Development](https://img.shields.io/badge/development-active-blue.svg)](https://img.shields.io/badge/development-active-blue.svg)
 [![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/qtlpoly)](https://cran.r-project.org/package=qtlpoly)
+[![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/qtlpoly)](https://cran.r-project.org/package=qtlpoly)
 [![R-universe PolyVerse Status Badge](https://polyploids.r-universe.dev/badges/qtlpoly)](https://polyploids.r-universe.dev/badges/qtlpoly)
 [![CRAN_monthly_downloads](https://cranlogs.r-pkg.org/badges/qtlpoly)](https://cranlogs.r-pkg.org/badges/qtlpoly?color=orange)
 
@@ -47,7 +47,7 @@ Tutorials as well as simulated and real data set analyses will be listed here in
 
 # Related software
 
-* [Polyverse](https://polyploids.r-universe.dev/ui#builds) - the polyploid R universe (a Lindsay Clark's initiative)
+* [Polyverse](https://polyploids.r-universe.dev/builds) - the polyploid R universe (a Lindsay Clark's initiative)
 
 * Variant Calling
   *  [GBSapp: An automated pipeline for variant calling and filtering.](https://github.com/bodeolukolu/GBSapp)
@@ -67,7 +67,7 @@ Tutorials as well as simulated and real data set analyses will be listed here in
   * [MAPpoly: Genetic maps in complex autopolyploids with even ploidy levels](https://CRAN.R-project.org/package=mappoly)
   * [MDSMap: High Density Genetic Linkage Mapping using Multidimensional Scaling](https://CRAN.R-project.org/package=MDSMap)
   * [polymapR: Linkage Analysis in Outcrossing Polyploids](https://CRAN.R-project.org/package=polymapR)
-  * [TetraploidSNPMap: Linkage maps and mapping QTLs for autotetraploid species, using SNP dosage data.](https://www.bioss.ac.uk/knowledge/tetraploidmap/)
+  * [TetraploidSNPMap: Linkage maps and mapping QTLs for autotetraploid species, using SNP dosage data.](https://www.bioss.ac.uk/knowledge-exchange/software/TetraploidSNPMap)
 
 
 * Haplotype reconstruction
@@ -96,4 +96,4 @@ Pereira GS, Gemenet DC, Mollinari M, Olukolu BA, Wood JC, Mosquera V, Gruneberg
 
 Qu L, Guennel T, Marshall SL. 2013. “Linear score tests for variance components in linear mixed models and applications to genetic association studies.” Biometrics 69 (4): 883-892.
 
-Wickham H. 2016. “ggplot2: Elegant Graphics for Data Analysis.” Springer. [doi:10.1007/978-0-387-98141-3](https://www.springer.com/gp/book/9780387981413).
+Wickham H. 2016. “ggplot2: Elegant Graphics for Data Analysis.” Springer. [doi:10.1007/978-0-387-98141-3](https://link.springer.com/book/10.1007/978-0-387-98141-3).
diff --git a/man/fit_model2.Rd b/man/fit_model2.Rd
diff --git a/man/mmer_adapted.Rd b/man/mmer_adapted.Rd
diff --git a/man/read_data2.Rd b/man/read_data2.Rd
diff --git a/src/Makevars b/src/Makevars
@@ -8,7 +8,7 @@
 ##
 ## And with R 3.4.0, and RcppArmadillo 0.7.960.*, we turn C++11 on as OpenMP
 ## support within Armadillo prefers / requires it
-CXX_STD = CXX11
+# CXX_STD = CXX11
 
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) 
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
diff --git a/src/Makevars.win b/src/Makevars.win
@@ -8,7 +8,7 @@
 ##
 ## And with R 3.4.0, and RcppArmadillo 0.7.960.*, we turn C++11 on as OpenMP
 ## support within Armadillo prefers / requires it
-CXX_STD = CXX11
+# CXX_STD = CXX11
 
 PKG_CXXFLAGS = $(SHLIB_OPENMP_CXXFLAGS) 
 PKG_LIBS = $(SHLIB_OPENMP_CXXFLAGS) $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)