diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION deleted file mode 100644 index b91a0268..00000000 --- a/CRAN-SUBMISSION +++ /dev/null @@ -1,3 +0,0 @@ -Version: 0.3.3 -Date: 2023-01-02 05:12:51 UTC -SHA: 1ae3aae45632b59d6240e724dec4002cef41f64d diff --git a/DESCRIPTION b/DESCRIPTION index 09f34a47..e4f3f8ba 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: mappoly Type: Package Title: Genetic Linkage Maps in Autopolyploids -Version: 0.4.0 +Version: 0.4.1 Authors@R: c(person(given = "Marcelo", family = "Mollinari", role = c("aut", "cre"), diff --git a/MAPpoly.Rproj b/MAPpoly.Rproj index 9b14df21..29e86a75 100755 --- a/MAPpoly.Rproj +++ b/MAPpoly.Rproj @@ -14,7 +14,7 @@ LaTeX: pdfLaTeX BuildType: Package PackageUseDevtools: Yes -PackageInstallArgs: --no-multiarch --with-keep.source --resave-data +PackageInstallArgs: --no-multiarch --with-keep.source --resave-data --use-LTO PackageBuildArgs: --resave-data PackageCheckArgs: --no-vignettes --as-cran --install-args=--build PackageRoxygenize: rd,collate,namespace,vignette diff --git a/NAMESPACE b/NAMESPACE index 25f5c5a4..7aed1ad2 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -113,7 +113,6 @@ export(split_and_rephase) export(summary_maps) export(update_framework_map) export(update_map) -export(update_missing) import(RCurl) import(Rcpp) import(fields) diff --git a/NEWS.md b/NEWS.md index a89dbd2b..a21b37eb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +# MAPpoly version 0.4.1 + - Addressed the `--use-LTO` installation failure issue + # MAPpoly version 0.4.0 - Functions to build maps in individual parents - Functions to merge individual maps diff --git a/R/RcppExports.R b/R/RcppExports.R index da90effa..a591069e 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,22 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +calc_genoprob_cpp <- function(m, geno, ph1, ph2, rf, probs, verbose) { + .Call('_mappoly_calc_genoprob_cpp', PACKAGE = 'mappoly', m, geno, ph1, ph2, rf, probs, verbose) +} + +calc_genprob_haplo_cpp <- function(m, n_mrk, n_ind, haplo, emit, rf, probs, verbose) { + .Call('_mappoly_calc_genprob_haplo_cpp', PACKAGE = 'mappoly', m, n_mrk, n_ind, haplo, emit, rf, probs, verbose) +} + +calc_genprob_haplo_highprec_cpp <- function(m, n_mrk, n_ind, haplo, emit, rf, probs, verbose) { + .Call('_mappoly_calc_genprob_haplo_highprec_cpp', PACKAGE = 'mappoly', m, n_mrk, n_ind, haplo, emit, rf, probs, verbose) +} + +loglike_hmm_cpp <- function(m, geno, ph1, ph2, rf, verbose) { + .Call('_mappoly_loglike_hmm_cpp', PACKAGE = 'mappoly', m, geno, ph1, ph2, rf, verbose) +} + .vcf_get_probabilities <- function(mat, pl_pos) { .Call('_mappoly_vcf_get_probabilities', PACKAGE = 'mappoly', mat, pl_pos) } diff --git a/R/calc_genoprob.R b/R/calc_genoprob.R index 84c1c503..dd271f50 100755 --- a/R/calc_genoprob.R +++ b/R/calc_genoprob.R @@ -96,17 +96,17 @@ calc_genoprob <- function(input.map, step = 0, phase.config = "best", verbose = phQ[rownames(Dtemp)] <- phQtemp seq.rf.pseudo <- mf_h(diff(map.pseudo)) } - for (j in 1:nrow(D)) - D[j, D[j, ] == input.map$info$ploidy + 1] <- dp[j] + dq[j] + 1 + as.numeric(dp[j] == 0 || dq[j] == 0) - res.temp <- .Call("calc_genoprob", - ploidy, - t(D), - phP, - phQ, - seq.rf.pseudo, - as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), - verbose = verbose, - PACKAGE = "mappoly") + for (j in 1:nrow(D)){ + D[j, D[j, ] == input.map$info$ploidy + 1] <- dp[j] + dq[j] + 1 + + as.numeric(dp[j] == 0 || dq[j] == 0) + } + res.temp <- calc_genoprob_cpp(ploidy, + t(D), + phP, + phQ, + seq.rf.pseudo, + as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), + verbose = verbose) if(verbose) cat("\n") dim(res.temp[[1]]) <- c(choose(ploidy,ploidy/2)^2,n.mrk,n.ind) dimnames(res.temp[[1]]) <- list(kronecker(apply(combn(letters[1:ploidy],ploidy/2),2, paste, collapse = ""), diff --git a/R/haplotype_map_utils.R b/R/haplotype_map_utils.R index 98374fd3..334a2eaa 100644 --- a/R/haplotype_map_utils.R +++ b/R/haplotype_map_utils.R @@ -258,27 +258,23 @@ calc_genoprob_haplo <- function(ploidy, n.mrk, n.ind, haplo, emit = NULL, } mrk.names <- names(haplo) if(highprec){ - res.temp <- .Call("calc_genprob_haplo_highprec", - ploidy, - n.mrk, - n.ind, - haplo, - emit, - rf_vec, - as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), - verbose, - PACKAGE = "mappoly") + res.temp <- calc_genprob_haplo_highprec_cpp(ploidy, + n.mrk, + n.ind, + haplo, + emit, + rf_vec, + as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), + verbose) } else{ - res.temp <- .Call("calc_genprob_haplo", - ploidy, - n.mrk, - n.ind, - haplo, - emit, - rf_vec, - as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), - verbose, - PACKAGE = "mappoly") + res.temp <- calc_genprob_haplo_cpp(ploidy, + n.mrk, + n.ind, + haplo, + emit, + rf_vec, + as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)), + verbose) } if(verbose) cat("\n") dim(res.temp[[1]]) <- c(choose(ploidy,ploidy/2)^2,n.mrk,n.ind) diff --git a/R/loglike_hmm.R b/R/loglike_hmm.R index 7806115a..40658608 100644 --- a/R/loglike_hmm.R +++ b/R/loglike_hmm.R @@ -43,14 +43,12 @@ loglike_hmm <- function(input.map, input.data = NULL, verbose = FALSE) D[j, D[j, ] == input.map$info$ploidy + 1] <- dp[j] + dq[j] + 1 + as.numeric(dp[j] == 0 || dq[j] == 0) for(i in 1:length(input.map$maps)){ rf.temp <- input.map$maps[[i]]$seq.rf - res.temp <- .Call("loglike_hmm", - input.map$info$ploidy, - t(D), - lapply(input.map$maps[[i]]$seq.ph$P, function(x) x-1), - lapply(input.map$maps[[i]]$seq.ph$Q, function(x) x-1), - rf.temp, - verbose, - PACKAGE = "mappoly") + res.temp <- loglike_hmm_cpp(input.map$info$ploidy, + t(D), + lapply(input.map$maps[[i]]$seq.ph$P, function(x) x-1), + lapply(input.map$maps[[i]]$seq.ph$Q, function(x) x-1), + rf.temp, + verbose) input.map$maps[[i]]$loglike <- res.temp[[1]] } return(input.map) diff --git a/R/utils.R b/R/utils.R index 4f52af82..36221c97 100755 --- a/R/utils.R +++ b/R/utils.R @@ -478,21 +478,7 @@ mp_pallet2 <- colorRampPalette(c("#911eb4", "#000075","#4363d8","#42d4f4","#4699 mp_pallet3 <- colorRampPalette(c("#ffe119", "#f58231","#e6194b","#808000","#9a6324", "#800000","#911eb4", "#000075","#4363d8","#42d4f4","#469990", "#3cb44b")) #' Update missing information -#' -#' Updates the missing data in the dosage matrix of an object of class -#' \code{mappoly.data} given a new probability threshold -#' @param input.data an object of class \code{mappoly.data} -#' @param prob.thres probability threshold to associate a marker call to a -#' dosage. Markers with maximum genotype probability smaller than 'prob.thres' -#' are considered as missing data for the dosage calling purposes -#' @examples -#' \donttest{ -#' data.updated = update_missing(hexafake.geno.dist, prob.thres = 0.5) -#' print(hexafake.geno.dist) -#' print(data.updated) -#' } -#' @author Marcelo Mollinari, \email{mmollin@ncsu.edu} -#' @export +#' @keywords internal update_missing <- function(input.data, prob.thres = 0.95){ geno.dose <- dist_prob_to_class(geno = input.data$geno, prob.thres = prob.thres) if(geno.dose$flag) diff --git a/cran-comments.md b/cran-comments.md index e538781e..e05ff5df 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,32 +1,51 @@ ## Re-submission of MAPpoly package -This is a re-submission of MAPpoly package. This version (0.4.0), contains the following changes +This is a re-submission of the MAPpoly package. This version includes changes to address the issues that led to its archival on CRAN, specifically the failure in installing the package using `--use-LTO`. We have made corrections to ensure compatibility and successful installation with `--use-LTO`. - - Update DESCRIPTION file - - Functions to build maps in individual parents - - Functions to merge individual maps - - Imputation functions based on map - - Functions to edit order interactively - - Fix minor bugs +Changes in this submission include: + + - Addressed the `--use-LTO` installation failure issue Thank you for reviewing our re-submission! ## Test environments -* local R installation (macOS Somona 14.3.1), R 4.3.2 -* local R installation (Ubuntu 22.04), R 4.3.2 +* local R installation (macOS Sonoma 14.3.1), R 4.3.3 + - Apple clang version 14.0.0 (clang-1400.0.29.202) + - GNU Fortran (GCC) 12.2.0 +* local R installation (Ubuntu 22.04), R 4.3.3 + - gcc (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0 + - GNU Fortran (Ubuntu 9.4.0-1ubuntu1~20.04.2) 9.4.0 + - Compiled using the following configuration: + + ./configure \ + --prefix=/opt/R/${R_VERSION} \ + --enable-R-shlib \ + --enable-memory-profiling \ + --with-blas \ + --with-lapack\ + --enable-lto=R + + And installed using `R CMD INSTALL --use-LTO mappoly_0.4.1.tar.gz` +* local R installation (Ubuntu 22.04), R devel + - Intel(R) oneAPI DPC++/C++ Compiler 2024.0.2 (2024.0.2.20231213) + - ifx (IFX) 2024.0.2 20231213 + * Win-builder (release and devel) + ## R CMD check results 0 errors | 0 warnings | 2 notes - - on local macOS: R 4.3.2 + - on local macOS: R 4.3.3 * installed size is 8.2Mb sub-directories of 1Mb or more: * R: 2.8Mb * data: 3.0Mb - * GNU make is a SystemRequirements. - + * GNU make is a SystemRequirements. + + - Checked with `--use-LTO` flag: Installation successful without any errors or warnings. + ## Downstream dependencies We checked 3 reverse dependencies (3 from CRAN + 0 from Bioconductor), comparing R CMD check results across CRAN and dev versions of this package. diff --git a/man/update_missing.Rd b/man/update_missing.Rd index ed972158..2b773166 100644 --- a/man/update_missing.Rd +++ b/man/update_missing.Rd @@ -6,24 +6,7 @@ \usage{ update_missing(input.data, prob.thres = 0.95) } -\arguments{ -\item{input.data}{an object of class \code{mappoly.data}} - -\item{prob.thres}{probability threshold to associate a marker call to a -dosage. Markers with maximum genotype probability smaller than 'prob.thres' -are considered as missing data for the dosage calling purposes} -} \description{ -Updates the missing data in the dosage matrix of an object of class -\code{mappoly.data} given a new probability threshold -} -\examples{ -\donttest{ -data.updated = update_missing(hexafake.geno.dist, prob.thres = 0.5) -print(hexafake.geno.dist) -print(data.updated) -} -} -\author{ -Marcelo Mollinari, \email{mmollin@ncsu.edu} +Update missing information } +\keyword{internal} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f10427fb..5853cd03 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -10,6 +10,75 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// calc_genoprob_cpp +List calc_genoprob_cpp(int m, NumericMatrix geno, List ph1, List ph2, NumericVector rf, std::vector probs, int verbose); +RcppExport SEXP _mappoly_calc_genoprob_cpp(SEXP mSEXP, SEXP genoSEXP, SEXP ph1SEXP, SEXP ph2SEXP, SEXP rfSEXP, SEXP probsSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type m(mSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type geno(genoSEXP); + Rcpp::traits::input_parameter< List >::type ph1(ph1SEXP); + Rcpp::traits::input_parameter< List >::type ph2(ph2SEXP); + Rcpp::traits::input_parameter< NumericVector >::type rf(rfSEXP); + Rcpp::traits::input_parameter< std::vector >::type probs(probsSEXP); + Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(calc_genoprob_cpp(m, geno, ph1, ph2, rf, probs, verbose)); + return rcpp_result_gen; +END_RCPP +} +// calc_genprob_haplo_cpp +List calc_genprob_haplo_cpp(int m, int n_mrk, int n_ind, List haplo, List emit, NumericVector rf, std::vector probs, int verbose); +RcppExport SEXP _mappoly_calc_genprob_haplo_cpp(SEXP mSEXP, SEXP n_mrkSEXP, SEXP n_indSEXP, SEXP haploSEXP, SEXP emitSEXP, SEXP rfSEXP, SEXP probsSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type m(mSEXP); + Rcpp::traits::input_parameter< int >::type n_mrk(n_mrkSEXP); + Rcpp::traits::input_parameter< int >::type n_ind(n_indSEXP); + Rcpp::traits::input_parameter< List >::type haplo(haploSEXP); + Rcpp::traits::input_parameter< List >::type emit(emitSEXP); + Rcpp::traits::input_parameter< NumericVector >::type rf(rfSEXP); + Rcpp::traits::input_parameter< std::vector >::type probs(probsSEXP); + Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(calc_genprob_haplo_cpp(m, n_mrk, n_ind, haplo, emit, rf, probs, verbose)); + return rcpp_result_gen; +END_RCPP +} +// calc_genprob_haplo_highprec_cpp +List calc_genprob_haplo_highprec_cpp(int m, int n_mrk, int n_ind, List haplo, List emit, NumericVector rf, std::vector probs, int verbose); +RcppExport SEXP _mappoly_calc_genprob_haplo_highprec_cpp(SEXP mSEXP, SEXP n_mrkSEXP, SEXP n_indSEXP, SEXP haploSEXP, SEXP emitSEXP, SEXP rfSEXP, SEXP probsSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type m(mSEXP); + Rcpp::traits::input_parameter< int >::type n_mrk(n_mrkSEXP); + Rcpp::traits::input_parameter< int >::type n_ind(n_indSEXP); + Rcpp::traits::input_parameter< List >::type haplo(haploSEXP); + Rcpp::traits::input_parameter< List >::type emit(emitSEXP); + Rcpp::traits::input_parameter< NumericVector >::type rf(rfSEXP); + Rcpp::traits::input_parameter< std::vector >::type probs(probsSEXP); + Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(calc_genprob_haplo_highprec_cpp(m, n_mrk, n_ind, haplo, emit, rf, probs, verbose)); + return rcpp_result_gen; +END_RCPP +} +// loglike_hmm_cpp +List loglike_hmm_cpp(int m, NumericMatrix geno, List ph1, List ph2, NumericVector rf, int verbose); +RcppExport SEXP _mappoly_loglike_hmm_cpp(SEXP mSEXP, SEXP genoSEXP, SEXP ph1SEXP, SEXP ph2SEXP, SEXP rfSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type m(mSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type geno(genoSEXP); + Rcpp::traits::input_parameter< List >::type ph1(ph1SEXP); + Rcpp::traits::input_parameter< List >::type ph2(ph2SEXP); + Rcpp::traits::input_parameter< NumericVector >::type rf(rfSEXP); + Rcpp::traits::input_parameter< int >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(loglike_hmm_cpp(m, geno, ph1, ph2, rf, verbose)); + return rcpp_result_gen; +END_RCPP +} // vcf_get_probabilities Rcpp::List vcf_get_probabilities(Rcpp::StringMatrix& mat, int pl_pos); RcppExport SEXP _mappoly_vcf_get_probabilities(SEXP matSEXP, SEXP pl_posSEXP) { @@ -59,32 +128,29 @@ BEGIN_RCPP END_RCPP } -RcppExport SEXP calc_genoprob(void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP calc_genoprob_prior(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP calc_genprob_haplo(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP calc_genprob_haplo_highprec(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP calc_genprob_single_parent(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP est_haplotype_map(void *, void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP est_haplotype_map_highprec(void *, void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP est_hmm_map_single_parent(void *, void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP est_map_hmm(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP est_map_hmm_highprec(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP get_counts_single_parent_cpp(void *, void *, void *, void *, void *, void *); -RcppExport SEXP loglike_hmm(void *, void *, void *, void *, void *, void *); -RcppExport SEXP pairwise_rf_estimation_disc(void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP pairwise_rf_estimation_disc_rcpp(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP pairwise_rf_estimation_prob(void *, void *, void *, void *, void *, void *, void *, void *); -RcppExport SEXP poly_hmm_est_CPP(void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *, void *); +RcppExport SEXP calc_genoprob_prior(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP calc_genprob_single_parent(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP est_haplotype_map(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP est_haplotype_map_highprec(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP est_hmm_map_single_parent(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP est_map_hmm(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP est_map_hmm_highprec(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP get_counts_single_parent_cpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP pairwise_rf_estimation_disc(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP pairwise_rf_estimation_disc_rcpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP pairwise_rf_estimation_prob(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +RcppExport SEXP poly_hmm_est_CPP(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); static const R_CallMethodDef CallEntries[] = { + {"_mappoly_calc_genoprob_cpp", (DL_FUNC) &_mappoly_calc_genoprob_cpp, 7}, + {"_mappoly_calc_genprob_haplo_cpp", (DL_FUNC) &_mappoly_calc_genprob_haplo_cpp, 8}, + {"_mappoly_calc_genprob_haplo_highprec_cpp", (DL_FUNC) &_mappoly_calc_genprob_haplo_highprec_cpp, 8}, + {"_mappoly_loglike_hmm_cpp", (DL_FUNC) &_mappoly_loglike_hmm_cpp, 6}, {"_mappoly_vcf_get_probabilities", (DL_FUNC) &_mappoly_vcf_get_probabilities, 2}, {"_mappoly_vcf_transform_dosage", (DL_FUNC) &_mappoly_vcf_transform_dosage, 2}, {"_mappoly_vcf_get_ploidy", (DL_FUNC) &_mappoly_vcf_get_ploidy, 2}, {"_mappoly_vcf_get_depth", (DL_FUNC) &_mappoly_vcf_get_depth, 2}, - {"calc_genoprob", (DL_FUNC) &calc_genoprob, 7}, {"calc_genoprob_prior", (DL_FUNC) &calc_genoprob_prior, 12}, - {"calc_genprob_haplo", (DL_FUNC) &calc_genprob_haplo, 8}, - {"calc_genprob_haplo_highprec", (DL_FUNC) &calc_genprob_haplo_highprec, 8}, {"calc_genprob_single_parent", (DL_FUNC) &calc_genprob_single_parent, 8}, {"est_haplotype_map", (DL_FUNC) &est_haplotype_map, 9}, {"est_haplotype_map_highprec", (DL_FUNC) &est_haplotype_map_highprec, 9}, @@ -92,7 +158,6 @@ static const R_CallMethodDef CallEntries[] = { {"est_map_hmm", (DL_FUNC) &est_map_hmm, 8}, {"est_map_hmm_highprec", (DL_FUNC) &est_map_hmm_highprec, 8}, {"get_counts_single_parent_cpp", (DL_FUNC) &get_counts_single_parent_cpp, 6}, - {"loglike_hmm", (DL_FUNC) &loglike_hmm, 6}, {"pairwise_rf_estimation_disc", (DL_FUNC) &pairwise_rf_estimation_disc, 7}, {"pairwise_rf_estimation_disc_rcpp", (DL_FUNC) &pairwise_rf_estimation_disc_rcpp, 12}, {"pairwise_rf_estimation_prob", (DL_FUNC) &pairwise_rf_estimation_prob, 8}, diff --git a/src/calc_genoprob.cpp b/src/calc_genoprob.cpp index 24c53038..2c178bd2 100755 --- a/src/calc_genoprob.cpp +++ b/src/calc_genoprob.cpp @@ -51,23 +51,15 @@ using namespace std; using namespace Rcpp; -RcppExport SEXP calc_genoprob(SEXP ploidyR, - SEXP genoR, - SEXP phPR, - SEXP phQR, - SEXP rfR, - SEXP probsR, - SEXP verboseR) +// [[Rcpp::export]] +List calc_genoprob_cpp(int m, + NumericMatrix geno, + List ph1, + List ph2, + NumericVector rf, + std::vector probs, + int verbose) { - //*convert input to C++ types - int m = Rcpp::as(ploidyR); - Rcpp::NumericMatrix geno(genoR); //(n.ind x n.col) - Rcpp::List ph1(phPR); - Rcpp::List ph2(phQR); - Rcpp::NumericVector rf(rfR); - std::vector probs = Rcpp::as >(probsR); - int verbose = Rcpp::as(verboseR); - //*Initializing some variables int g = nChoosek(m, m/2), k, k1; int n_mrk = geno.ncol(); // markers are disposed in columns diff --git a/src/calc_genoprob.h b/src/calc_genoprob.h index 16455270..668d6a11 100755 --- a/src/calc_genoprob.h +++ b/src/calc_genoprob.h @@ -17,10 +17,13 @@ For a copy of the GNU General Public License, please visit . */ -RcppExport SEXP calc_genoprob(SEXP ploidyR, - SEXP genoR, - SEXP phPR, - SEXP phQR, - SEXP rfR, - SEXP probsR, - SEXP verboseR); +using namespace std; +using namespace Rcpp; + +List calc_genoprob_cpp(int m, + NumericMatrix geno, + List ph1, + List ph2, + NumericVector rf, + std::vector probs, + int verbose); diff --git a/src/calc_genoprob_based_on_phased_marker_blocks.cpp b/src/calc_genoprob_based_on_phased_marker_blocks.cpp index 2242ab7a..d8978aae 100644 --- a/src/calc_genoprob_based_on_phased_marker_blocks.cpp +++ b/src/calc_genoprob_based_on_phased_marker_blocks.cpp @@ -51,25 +51,16 @@ using namespace std; using namespace Rcpp; -RcppExport SEXP calc_genprob_haplo(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP probsR, - SEXP verboseR) +// [[Rcpp::export]] +List calc_genprob_haplo_cpp(int m, + int n_mrk, + int n_ind, + List haplo, + List emit, + NumericVector rf, + std::vector probs, + int verbose) { - //convert input to C++ types - int m = Rcpp::as(ploidyR); - int n_mrk = Rcpp::as(n_mrkR); - int n_ind = Rcpp::as(n_indR); - Rcpp::List haplo(haploR); - Rcpp::List emit(emitR); - Rcpp::NumericVector rf(rfR); - int verbose = Rcpp::as(verboseR); - std::vector probs = Rcpp::as >(probsR); - //Initializing some variables int g = nChoosek(m, m/2), k, k1, count = 0; std::vector term(n_ind); @@ -175,27 +166,16 @@ RcppExport SEXP calc_genprob_haplo(SEXP ploidyR, List z = List::create(probs); return z ; } - - -RcppExport SEXP calc_genprob_haplo_highprec(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP probsR, - SEXP verboseR) +// [[Rcpp::export]] +List calc_genprob_haplo_highprec_cpp(int m, + int n_mrk, + int n_ind, + List haplo, + List emit, + NumericVector rf, + std::vector probs, + int verbose) { - //convert input to C++ types - int m = Rcpp::as(ploidyR); - int n_mrk = Rcpp::as(n_mrkR); - int n_ind = Rcpp::as(n_indR); - Rcpp::List haplo(haploR); - Rcpp::List emit(emitR); - Rcpp::NumericVector rf(rfR); - int verbose = Rcpp::as(verboseR); - std::vector probs = Rcpp::as >(probsR); - //Initializing some variables int g = nChoosek(m, m/2), k, k1, count = 0; std::vector term(n_ind); diff --git a/src/calc_genoprob_based_on_phased_marker_blocks.h b/src/calc_genoprob_based_on_phased_marker_blocks.h index a1ded9e3..bb66cc75 100644 --- a/src/calc_genoprob_based_on_phased_marker_blocks.h +++ b/src/calc_genoprob_based_on_phased_marker_blocks.h @@ -31,20 +31,20 @@ using namespace std; using namespace Rcpp; -RcppExport SEXP calc_genprob_haplo(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP probsR, - SEXP verboseR); +List calc_genprob_haplo_cpp(int m, + int n_mrk, + int n_ind, + List haplo, + List emit, + NumericVector rf, + std::vector probs, + int verbose); -RcppExport SEXP calc_genprob_haplo_highprec(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP probsR, - SEXP verboseR); \ No newline at end of file +List calc_genprob_haplo_highprec(int m, + int n_mrk, + int n_ind, + List haplo, + List emit, + NumericVector rf, + std::vector probs, + int verbose); \ No newline at end of file diff --git a/src/calc_loglike_given_map.cpp b/src/calc_loglike_given_map.cpp index 18a914e6..0c00d584 100644 --- a/src/calc_loglike_given_map.cpp +++ b/src/calc_loglike_given_map.cpp @@ -49,21 +49,14 @@ using namespace std; using namespace Rcpp; -RcppExport SEXP loglike_hmm(SEXP ploidyR, - SEXP genoR, - SEXP phPR, - SEXP phQR, - SEXP rfR, - SEXP verboseR) +// [[Rcpp::export]] +List loglike_hmm_cpp(int m, + NumericMatrix geno, + List ph1, + List ph2, + NumericVector rf, + int verbose) { - //*convert input to C++ types - int m = Rcpp::as(ploidyR); - Rcpp::NumericMatrix geno(genoR); //(n.ind x n.col) - Rcpp::List ph1(phPR); - Rcpp::List ph2(phQR); - Rcpp::NumericVector rf(rfR); - int verbose = Rcpp::as(verboseR); - //*Initializing some variables int g = nChoosek(m, m/2); int n_mar = geno.ncol(); // markers are disposed in columns diff --git a/src/calc_loglike_given_map.h b/src/calc_loglike_given_map.h index 452533c6..3ae6effc 100644 --- a/src/calc_loglike_given_map.h +++ b/src/calc_loglike_given_map.h @@ -28,9 +28,12 @@ #include #include "combinatorial.h" -RcppExport SEXP loglike_hmm(SEXP ploidyR, - SEXP genoR, - SEXP phPR, - SEXP phQR, - SEXP rfR, - SEXP verboseR); \ No newline at end of file +using namespace std; +using namespace Rcpp; + +List loglike_hmm_cpp(int m, + NumericMatrix geno, + List ph1, + List ph2, + NumericVector rf, + int verbose); \ No newline at end of file diff --git a/src/hmm_elements.cpp b/src/hmm_elements.cpp index 27c41f48..1fd5ce1a 100755 --- a/src/hmm_elements.cpp +++ b/src/hmm_elements.cpp @@ -259,15 +259,17 @@ std::vector > index_func(int m, do { s=0; - for(int j=0; (unsigned)j < p.size(); j++) - if(p[j]>=0) s=s+vp[p[j]]; - for(int j=0; (unsigned)j < q.size(); j++) - if(q[j]>=0) s=s+vq[q[j]]; - v1[s].push_back(ip); - v2[s].push_back(iq); - v1[v1.size()-1].push_back(ip); - v2[v2.size()-1].push_back(iq); - iq++; + for(int j=0; (unsigned)j < p.size(); j++){ + if(p[j]>=0) s=s+vp[p[j]]; + } + for(int j=0; (unsigned)j < q.size(); j++){ + if(q[j]>=0) s=s+vq[q[j]]; + } + v1[s].push_back(ip); + v2[s].push_back(iq); + v1[v1.size()-1].push_back(ip); + v2[v2.size()-1].push_back(iq); + iq++; } while (std::prev_permutation(vq.begin(), vq.end())); ip++; @@ -400,11 +402,11 @@ std::vector backward_emit(int m, Classical forward equation presented in Rabiner 1989. */ std::vector forward_emit_single_parent(int m, - std::vector& fk, - std::vector& ik, - std::vector& ik1, - std::vector& emit, - std::vector >& T) + std::vector& fk, + std::vector& ik, + std::vector& ik1, + std::vector& emit, + std::vector >& T) { int ngenk = ik.size(); int ngenk1 = ik1.size(); @@ -425,11 +427,11 @@ std::vector forward_emit_single_parent(int m, Classical backward equation presented in Rabiner 1989. */ std::vector backward_emit_single_parent(int m, - std::vector& fk1, - std::vector& ik, - std::vector& ik1, - std::vector& emit, - std::vector >& T) + std::vector& fk1, + std::vector& ik, + std::vector& ik1, + std::vector& emit, + std::vector >& T) { int ngenk = ik.size(); int ngenk1 = ik1.size(); diff --git a/src/hmm_single_parent.cpp b/src/hmm_single_parent.cpp index 5f686303..372cc6e3 100644 --- a/src/hmm_single_parent.cpp +++ b/src/hmm_single_parent.cpp @@ -53,14 +53,14 @@ using namespace std; using namespace Rcpp; RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP verboseR, - SEXP tolR, - SEXP ret_H0R) + SEXP n_mrkR, + SEXP n_indR, + SEXP haploR, + SEXP emitR, + SEXP rfR, + SEXP verboseR, + SEXP tolR, + SEXP ret_H0R) { //convert input to C++ types int m = Rcpp::as(ploidyR); @@ -126,7 +126,7 @@ RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, //Initializing recombination number matrix std::vector< std::vector > R; R = rec_num(m); - + //begin EM algorithm for(int it=0; it 0) // Verify theoretical implications of this condition - rf[k] += nr * gamma[i][j]/s; + rf[k] += nr * gamma[i][j]/s; } } } @@ -203,7 +203,7 @@ RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, term[ind] += alpha[ind][n_mrk-1][j]; } } // loop over individuals - + //Likelihood using a specific recombination fraction vector //Usually, this is used to compute LOD Score under H0: rf=0.5 if(ret_H0 == 1) @@ -212,10 +212,12 @@ RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, for(int i=0; (unsigned)i < alpha.size(); i++) { temp=0.0; - for(int j=0; (unsigned)j < alpha[i][n_mrk-1].size(); j++) - temp += alpha[i][n_mrk-1][j]; - if(temp > 0) - loglike += log10(temp); + for(int j=0; (unsigned)j < alpha[i][n_mrk-1].size(); j++){ + temp += alpha[i][n_mrk-1][j]; + } + if(temp > 0){ + loglike += log10(temp); + } } if(verbose) Rcpp::Rcout << "\n"; @@ -251,15 +253,17 @@ RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, if(!flag) break; }//end of EM algorithm if(flag && verbose) Rcpp::Rcout << "Didn't converge!\n"; - + //Loglike computation for(int i=0; (unsigned)i < alpha.size(); i++) { temp=0.0; - for(int j=0; (unsigned)j < alpha[i][n_mrk-1].size(); j++) - temp += alpha[i][n_mrk-1][j]; - if(temp > 0) - loglike += log10(temp); + for(int j=0; (unsigned)j < alpha[i][n_mrk-1].size(); j++){ + temp += alpha[i][n_mrk-1][j]; + } + if(temp > 0){ + loglike += log10(temp); + } } @@ -281,13 +285,13 @@ RcppExport SEXP est_hmm_map_single_parent(SEXP ploidyR, } RcppExport SEXP calc_genprob_single_parent(SEXP ploidyR, - SEXP n_mrkR, - SEXP n_indR, - SEXP haploR, - SEXP emitR, - SEXP rfR, - SEXP probsR, - SEXP verboseR) + SEXP n_mrkR, + SEXP n_indR, + SEXP haploR, + SEXP emitR, + SEXP rfR, + SEXP probsR, + SEXP verboseR) { //convert input to C++ types int m = Rcpp::as(ploidyR); @@ -376,7 +380,7 @@ RcppExport SEXP calc_genprob_single_parent(SEXP ploidyR, } } // loop over individuals - + // for(int k = 0; k < n_mrk; k++) // { // for(int j = 0; j < alpha[1][k].size(); j++){ diff --git a/src/pairwise_estimation_rcppparallel.cpp b/src/pairwise_estimation_rcppparallel.cpp index 1ec3dcc8..713e4a0e 100755 --- a/src/pairwise_estimation_rcppparallel.cpp +++ b/src/pairwise_estimation_rcppparallel.cpp @@ -123,7 +123,7 @@ struct rf_two_point_parallel : public Worker { int start = 0; int rn = 1; // Looping through string - for(rn = 1; rn < mystring.size(); rn++){ + for(std::string::size_type rn = 1; rn < mystring.size(); rn++) { if( mystring[rn] == split){ std::string temp = mystring.substr(start, (rn - start)); z.push_back(temp); @@ -142,7 +142,7 @@ struct rf_two_point_parallel : public Worker { start = 0; rn = 1; // Looping through string - for(rn = 1; rn < mystring.size(); rn++){ + for(std::string::size_type rn = 1; rn < mystring.size(); rn++) { if( mystring[rn] == split){ std::string temp = mystring.substr(start, (rn - start)); phases.push_back(temp); @@ -154,7 +154,7 @@ struct rf_two_point_parallel : public Worker { phases.push_back(temp); std::vector phase1(phases.size()), phase2(phases.size()); std::string delimiter2 = "-"; - for(int j=0; j < phases.size(); j++) + for(std::string::size_type j=0; j < phases.size(); j++) { std::string lnames = phases[j]; phase1[j] = std::stoi(lnames.substr(0,lnames.find(delimiter2))); @@ -172,7 +172,7 @@ struct rf_two_point_parallel : public Worker { ((count_matrix_pos[(id-1)]-1) + ((i+1)*len_ac - 1))); std::vector dk(z.size()), dk1(z.size()); std::string delimiter = " "; - for(int j=0; j < z.size(); j++) + for(std::string::size_type j=0; j < z.size(); j++) { std::string lnames = z[j]; dk[j] = std::stoi(lnames.substr(0,lnames.find(delimiter))); @@ -211,7 +211,7 @@ struct rf_two_point_parallel : public Worker { double temp=0.0; std::vector Tr((m+2) * (m+2)); std::fill(Tr.begin(), Tr.end(), 1); - for(int i = 0; i < dk.size(); i++){ + for(std::string::size_type i = 0; i < dk.size(); i++){ //Rcpp::Rcout << dk(i) << " " << dk1(i) << std::endl;; count=0; Tr[dk[i] + (dk1[i] * (m+2))]=0.0; @@ -293,7 +293,7 @@ struct rf_two_point_parallel : public Worker { count2=0; temp=0.0; std::fill(Tr.begin(), Tr.end(), 1); - for(int i = 0; i < dk.size(); i++){ + for(std::string::size_type i = 0; i < dk.size(); i++){ //Rcpp::Rcout << dk(i) << " " << dk1(i) << std::endl;; count=0; Tr[dk[i] + (dk1[i] * (m+2))]=0; @@ -338,7 +338,7 @@ struct rf_two_point_parallel : public Worker { count2=0; temp=0.0; std::fill(Tr.begin(), Tr.end(), 1); - for(int i = 0; i < dk.size(); i++){ + for(std::string::size_type i = 0; i < dk.size(); i++){ //Rcpp::Rcout << dk(i) << " " << dk1(i) << std::endl;; count=0; Tr[dk[i] + (dk1[i] * (m+2))]=0; @@ -369,7 +369,7 @@ struct rf_two_point_parallel : public Worker { count2=0; temp=0.0; std::fill(Tr.begin(), Tr.end(), 1); - for(int i = 0; i < dk.size(); i++){ + for(std::string::size_type i = 0; i < dk.size(); i++){ //Rcpp::Rcout << dk(i) << " " << dk1(i) << std::endl;; count=0; Tr[dk[i] + (dk1[i] * (m+2))]=0; diff --git a/src/read_mappoly_vcf.cpp b/src/read_mappoly_vcf.cpp index 8cda7a3e..138102a3 100755 --- a/src/read_mappoly_vcf.cpp +++ b/src/read_mappoly_vcf.cpp @@ -170,7 +170,7 @@ int get_depth(std::string mystring, int dp_pos){ std::string temp = mystring.substr(start, i - start); vec_o_strings.push_back(temp); // Returning DP - if (vec_o_strings.size() < dp_pos) { + if (vec_o_strings.size() < static_cast::size_type>(dp_pos)){ start = 0; } else if (vec_o_strings[dp_pos-1] == ".") { start = 0; @@ -206,7 +206,7 @@ std::vector get_probabilities(std::string mystring, int pl_pos){ std::string temp = mystring.substr(start, i - start); vec_o_strings.push_back(temp); // Checking PL information - if (vec_o_strings.size() < pl_pos){ + if (vec_o_strings.size() < static_cast::size_type>(pl_pos)){ final_vec_out2[0] = -1.0; return final_vec_out2; } @@ -273,7 +273,7 @@ Rcpp::List vcf_get_probabilities(Rcpp::StringMatrix& mat, int pl_pos){ for (l=0; l < colmat; l++){ // Rcout << "The value of l : " << l << "\n"; out = get_probabilities(as(mat(k,l)), pl_pos); - if (out.size() == size){ + if (out.size() == static_cast::size_type>(size)){ for (s=0; s < size; s++){ // Rcout << "The value of s : " << s << "\n"; partial_results(l,s) = out[s]; diff --git a/tests/testthat.R b/tests/testthat.R deleted file mode 100644 index 3292308b..00000000 --- a/tests/testthat.R +++ /dev/null @@ -1,5 +0,0 @@ -if(requireNamespace("testthat", quietly = TRUE)){ - library(testthat) - library(mappoly) - test_check("mappoly") -} diff --git a/tests/testthat/test-2pt.R b/tests/testthat/test-2pt.R deleted file mode 100644 index 6a64ffaa..00000000 --- a/tests/testthat/test-2pt.R +++ /dev/null @@ -1,37 +0,0 @@ -context("Two point estimates") -test_that("estimate two-points rf correctly", { - dat <- filter_missing(tetra.solcap.geno.dist, inter = FALSE) - expect_equal(mean(as.matrix(dat$geno.dose)), 1.92626, tolerance = 1e-4) - dat <- filter_missing(tetra.solcap.geno.dist, inter = FALSE, type = "individual") - expect_equal(mean(as.matrix(dat$geno.dose)), 1.93118, tolerance = 1e-4) - f <- filter_segregation(dat, inter = FALSE) - expect_equivalent(length(f$keep), 3678) - s <- make_seq_mappoly(tetra.solcap.geno.dist, f$keep) - s <- make_seq_mappoly(tetra.solcap.geno.dist, s$seq.mrk.names[seq(1,length(s$seq.num), 50)]) - tpt1 <- est_pairwise_rf(s) - expect_is(plot(tpt1, first.mrk = 1, second.mrk = 101), "ggplot") - tpt11 <- est_pairwise_rf(s, est.type = "prob") - expect_is(plot(tpt11, first.mrk = 1, second.mrk = 101), "ggplot") - m <- rf_list_to_matrix(tpt1) - m2 <- make_mat_mappoly(m, make_seq_mappoly(tetra.solcap.geno.dist, names(which(s$chrom==1)))) - expect_equal(sum(m2$rec.mat, na.rm = TRUE), 12.51488, tolerance = 10e-5) - m11 <- rf_list_to_matrix(tpt11) - m12 <- make_mat_mappoly(m11, make_seq_mappoly(tetra.solcap.geno.dist, names(which(s$chrom==1)))) - expect_equal(sum(m12$rec.mat, na.rm = TRUE), 12.54115, tolerance = 10e-5) - o <- mds_mappoly(m) - expect_is(o, "mappoly.pcmap") - expect_null(plot(m)) - expect_equivalent(round(mean(m[[4]], na.rm = T),6), 0.396882) - skip_on_os("mac") - tpt2 <- est_pairwise_rf(s, ncpus = 2) - expect_equal(tpt1, tpt2, tolerance = 1e-4) - tpt13 <- est_pairwise_rf(s, ncpus = 2, est.type = "prob") - expect_equal(tpt11, tpt13, tolerance = 1e-4) - sf<-rf_snp_filter(tpt2) - sf2<-rf_snp_filter(tpt13) - expect_is(sf, "mappoly.sequence") - expect_is(sf2, "mappoly.sequence") - expect_equal(sum(sf$seq.num), 135192) - tpt2 <- est_pairwise_rf2(s, ncpus = 2) - expect_equal(round(sum(tpt2$pairwise$rf - m$rec.mat, na.rm = TRUE),3), 0) -}) \ No newline at end of file diff --git a/tests/testthat/test-calc_genoprobs.R b/tests/testthat/test-calc_genoprobs.R deleted file mode 100644 index 3cbb5198..00000000 --- a/tests/testthat/test-calc_genoprobs.R +++ /dev/null @@ -1,17 +0,0 @@ -context("Conditional probabilities") -test_that("computes genotype probabilities correctly", { - x1 <- get_submap(solcap.dose.map[[1]], 1:20, reestimate.rf = F) - probs.t1<-calc_genoprob(input.map = x1, - verbose = TRUE) - expect_equal(var(probs.t1$probs), 0.0153, tolerance = 1e-3) - probs.t2<-calc_genoprob_error(input.map = x1, - verbose = TRUE, - error = 0.2) - expect_equal(var(probs.t2$probs), 0.0131, tolerance = 1e-3) - expect_error(calc_genoprob_dist(input.map = x1, - verbose = TRUE)) - x3 <- get_submap(solcap.prior.map[[1]], 1:20, reestimate.rf = F) - probs.t3<-calc_genoprob_dist(input.map = x3, - verbose = TRUE) - expect_equal(var(probs.t3$probs), 0.0153, tolerance = 1e-3) -}) \ No newline at end of file diff --git a/tests/testthat/test-est_hmm_map.R b/tests/testthat/test-est_hmm_map.R deleted file mode 100644 index b469094a..00000000 --- a/tests/testthat/test-est_hmm_map.R +++ /dev/null @@ -1,61 +0,0 @@ -context("Estimate HMM map") -test_that("map contructed correctly", { - ##### Tetraploid - s<-make_seq_mappoly(tetra.solcap, 1:2) - map <- est_rf_hmm(input.seq = make_seq_mappoly(tetra.solcap, 1:2), - twopt =est_pairwise_rf(s), - thres = 2, tol = 10e-4) - map <- loglike_hmm(map) - map <- est_full_hmm_with_global_error(map, error = .05) - map <- est_full_hmm_with_prior_prob(map, dat.prob = tetra.solcap.geno.dist) - s<-make_seq_mappoly(tetra.solcap, 3:7) - expect_is(s, "mappoly.sequence") - expect_output(str(s), "List of 14") - tpt<-est_pairwise_rf(s) - expect_is(tpt, "mappoly.twopt") - expect_output(str(tpt), "List of 7") - expect_output(str(tpt$pairwise), "List of 10") - expect_is(map, "mappoly.map") - map <- est_rf_hmm(input.seq = s, twopt = tpt, thres = 2, tol = 10e-4) - print(map, detailed = TRUE) - expect_is(map, "mappoly.map") - map <- loglike_hmm(map) - expect_output(str(map$info), "List of 13") - expect_equivalent(map$maps[[1]]$seq.rf, - c(6.787204e-03, 1.283112e-03, 1.137237e-03, 3.272807e-05)) - expect_equivalent(map$maps[[1]]$seq.ph, - list(P=list('3' = c(1, 2, 3, 4), '4' = c(1, 2, 3), '5' = c(1, 2, 3, 4), '6' = c(1, 2, 3), '7' = 0), - Q=list('3' = c(1, 2, 3), '4' = c(1, 2), '5' = c(1, 2, 3), '6' = c(1, 2), '7' = 4))) - expect_output(print(map), "This is an object of class 'mappoly.map'\\n Ploidy level:\\t 4 \\n No. individuals:\\t 160 \\n No. markers:\\t 5 \\n No. linkage phases:\\t 1 \\n\\n ---------------------------------------------\\n Number of linkage phase configurations: 1\\n ---------------------------------------------\\n Linkage phase configuration: 1\\n map length:\\t 0.93\\n log-likelihood:\\t -123.7\\n LOD:\\t\\t 0\\n ~~~~~~~~~~~~~~~~~~") - expect_length(plot(map), 2) - expect_length(plot(map, left.lim = .2, right.lim = .9), 2) - expect_equal(round(mean(plot_map_list(list(map, map))[,3]),6), 0.669977) - expect_equal(round(mean(plot_map_list(list(map, map), horiz = F, col = "ggstyle")[,3]),6), 0.669977) - expect_is(plot_genome_vs_map(list(map, map)), "ggplot") - expect_is(plot_genome_vs_map(list(map, map), same.ch.lg = TRUE), "ggplot") - expect_equivalent(summary_maps(list(map, map))[,3], c("0.93", "0.93", "1.86")) -}) -test_that("sequential map contructed correctly", { - s<-make_seq_mappoly(tetra.solcap, 1:2) - map <- est_rf_hmm_sequential(input.seq = s, - twopt = est_pairwise_rf(s), - verbose = FALSE) - ##### Tetraploid - s<-make_seq_mappoly(tetra.solcap, 3:7) - expect_is(s, "mappoly.sequence") - expect_output(str(s), "List of 14") - tpt<-est_pairwise_rf(s) - expect_is(tpt, "mappoly.twopt") - expect_output(str(tpt), "List of 7") - expect_output(str(tpt$pairwise), "List of 10") - map <- est_rf_hmm_sequential(input.seq = s, twopt = tpt, thres.twopt = 2, tol.final = 10e-4, detailed.verbose = TRUE) - map <- est_rf_hmm_sequential(input.seq = s, twopt = tpt, thres.twopt = 2, tol.final = 10e-4, high.prec = TRUE) - expect_is(map, "mappoly.map") - expect_output(str(map$info), "List of 13") - expect_equivalent(map$maps[[1]]$seq.rf, - c(6.787204e-03, 1.283112e-03, 1.137237e-03, 3.272807e-05)) - expect_equivalent(map$maps[[1]]$seq.ph, - list(P=list('3' = c(1, 2, 3, 4), '4' = c(1, 2, 3), '5' = c(1, 2, 3, 4), '6' = c(1, 2, 3), '7' = 0), - Q=list('3' = c(1, 2, 3), '4' = c(1, 2), '5' = c(1, 2, 3), '6' = c(1, 2), '7' = 4))) - expect_output(print(map), "This is an object of class 'mappoly.map'\\n Ploidy level:\\t 4 \\n No. individuals:\\t 160 \\n No. markers:\\t 5 \\n No. linkage phases:\\t 1 \\n\\n ---------------------------------------------\\n Number of linkage phase configurations: 1\\n ---------------------------------------------\\n Linkage phase configuration: 1\\n map length:\\t 0.93\\n log-likelihood:\\t -123.7\\n LOD:\\t\\t 0\\n ~~~~~~~~~~~~~~~~~~") -}) diff --git a/tests/testthat/test-est_with_probs.R b/tests/testthat/test-est_with_probs.R deleted file mode 100644 index 648e91f8..00000000 --- a/tests/testthat/test-est_with_probs.R +++ /dev/null @@ -1,14 +0,0 @@ -context("Estimate map with probabilities") -test_that("map estimated correctly", { - map <- get_submap(solcap.prior.map[[1]], c(1,3,5,7,9)) - mape <- est_full_hmm_with_global_error(map, to = 10e-4) - mapp <- est_full_hmm_with_prior_prob(map, to = 10e-4) - expect_equivalent(map$maps[[1]]$seq.rf, - c(0.037783070, 0.001634628, 0.015620733, 0.016022936)) - expect_equivalent(mape$maps[[1]]$seq.rf, - c(0.02923187, 0.00162774, 0.01649081, 0.01506356)) - expect_equivalent(mapp$maps[[1]]$seq.rf, - c(0.02923187, 0.00162774, 0.01649081, 0.01506356)) - map <- get_submap(solcap.dose.map[[1]], 1:5) - expect_error(est_full_hmm_with_prior_prob(map)) -}) \ No newline at end of file diff --git a/tests/testthat/test-export_map_list.R b/tests/testthat/test-export_map_list.R deleted file mode 100644 index 24fbd7fd..00000000 --- a/tests/testthat/test-export_map_list.R +++ /dev/null @@ -1,5 +0,0 @@ -context("Export map list") -test_that("export map list correctly", { - ##### Tetraploid - expect_invisible(export_map_list(solcap.dose.map[1:2], "")) -}) diff --git a/tests/testthat/test-filters.R b/tests/testthat/test-filters.R deleted file mode 100644 index 7cd6a4a7..00000000 --- a/tests/testthat/test-filters.R +++ /dev/null @@ -1,28 +0,0 @@ -context("Filter functions") -test_that("test filter functions", { - a1<-sample_data(tetra.solcap, n = 20, type = "marker") - expect_equal(check_data_sanity(a1), 0) - a2<-sample_data(tetra.solcap.geno.dist, n = 20, type = "marker") - expect_equal(check_data_sanity(a2), 0) - a3<-filter_missing(input.data = a1, type = "marker", - filter.thres = 0.1, inter = FALSE) - expect_equal(check_data_sanity(a3), 0) - a4<-filter_missing(input.data = a2, type = "marker", - filter.thres = 0.1, inter = FALSE) - expect_equal(check_data_sanity(a4), 0) - a5<-filter_missing(input.data = a1, type = "individual", - filter.thres = 0.1, inter = FALSE) - expect_equal(check_data_sanity(a5), 0) - a6<-filter_missing(input.data = a2, type = "marker", - filter.thres = 0.1, inter = FALSE) - expect_equal(check_data_sanity(a5), 0) - a1$geno.dose[10, 10] <- 8 - a7<-filter_non_conforming_classes(input.data = a1) - expect_equal(a7$geno.dose[10, 10],5) - w1<-filter_segregation(tetra.solcap, chisq.pval.thres = 0.001, inter = FALSE) - w2<-make_seq_mappoly(w1) - expect_equal(length(w2$seq.num), 3659) - w3<-filter_segregation(tetra.solcap, chisq.pval.thres = 1e-10, inter = FALSE) - w4<-make_seq_mappoly(w3) - expect_equal(length(w4$seq.num), length(tetra.solcap$mrk.names)) -}) \ No newline at end of file diff --git a/tests/testthat/test-get_submap.R b/tests/testthat/test-get_submap.R deleted file mode 100644 index 7c9b6294..00000000 --- a/tests/testthat/test-get_submap.R +++ /dev/null @@ -1,9 +0,0 @@ -context("Get submap") -test_that("sub-map extracted correctly", { - ##### Hexaploid - sub.map1<-get_submap(solcap.dose.map[[1]], 1:5) - expect_is(sub.map1, "mappoly.map") - expect_output(str(sub.map1$info), "List of 11") - expect_equivalent(sub.map1$maps[[1]]$seq.rf, - c(3.113371e-02, 1.000000e-05, 1.838496e-05, 8.084546e-03)) -}) \ No newline at end of file diff --git a/tests/testthat/test-group.R b/tests/testthat/test-group.R deleted file mode 100644 index 63c4b3ee..00000000 --- a/tests/testthat/test-group.R +++ /dev/null @@ -1,10 +0,0 @@ -context("Groupping") -test_that("assemble linkage groups correctly", { - x <- est_pairwise_rf(make_seq_mappoly(tetra.solcap, 1:50)) - M <- rf_list_to_matrix(x) - gr<-group_mappoly(M, expected.groups = 12, verbose = FALSE, inter = FALSE) - y<-plot(gr) - expect_is(gr, "mappoly.group") - expect_null(y, NULL) - expect_equal(round(mean(gr$groups.snp), 6), 4.32) -}) diff --git a/tests/testthat/test-homolog_prob.R b/tests/testthat/test-homolog_prob.R deleted file mode 100644 index d124187f..00000000 --- a/tests/testthat/test-homolog_prob.R +++ /dev/null @@ -1,15 +0,0 @@ -context("Homolog probability") -test_that("computes homolog probabilities correctly", { - x1 <- get_submap(solcap.dose.map[[1]], 1:20, reestimate.rf = F) - probs.t1<-calc_genoprob(input.map = x1, - verbose = TRUE) - expect_equal(var(probs.t1$probs), 0.0153, tolerance = 1e-3) - hom.t1 <- calc_homologprob(input.genoprobs = probs.t1) - expect_equal(var(hom.t1$homoprob$probability), 0.187, tolerance = 1e-3) - skip_if_not(capabilities("long.double")) - expect_is(plot(hom.t1), "plotly") - skip_if_not(capabilities("long.double")) - expect_is(plot(hom.t1, stack = T), "plotly") - skip_if_not(capabilities("long.double")) - expect_is(plot(hom.t1, use.plotly = FALSE), "ggplot") -}) \ No newline at end of file diff --git a/tests/testthat/test-merge_maps.R b/tests/testthat/test-merge_maps.R deleted file mode 100644 index f5c0aade..00000000 --- a/tests/testthat/test-merge_maps.R +++ /dev/null @@ -1,22 +0,0 @@ -context("Merge maps") -test_that("merging maps correctly", { - ##### Tetraploid - map1<-get_submap(solcap.dose.map[[1]], 1:5, reestimate.rf = FALSE, verbose = FALSE) - map2<-get_submap(solcap.dose.map[[1]], 6:15, reestimate.rf = FALSE, verbose = FALSE) - map3<-get_submap(solcap.dose.map[[1]], 16:30, reestimate.rf = FALSE, verbose = FALSE) - full.map<-get_submap(solcap.dose.map[[1]], 1:30, reestimate.rf = FALSE, verbose = FALSE) - s<-make_seq_mappoly(tetra.solcap, full.map$maps[[1]]$seq.num) - twopt <- est_pairwise_rf(input.seq = s) - merged.maps<-merge_maps(map.list = list(map1, map2, map3), - twopt = twopt, - thres.twopt = 3) - expect_equal(round(mean(unlist(merged.maps$maps[[1]]$seq.ph)),6), 2.495868) - best.phase <- merged.maps$maps[[1]]$seq.ph - names.id<-names(best.phase$P) - x1 <- compare_haplotypes(ploidy = 4, best.phase$P[names.id], - full.map$maps[[1]]$seq.ph$P[names.id]) - x2 <- compare_haplotypes(ploidy = 4, best.phase$Q[names.id], - full.map$maps[[1]]$seq.ph$Q[names.id]) - expect_true(x1$is.same.haplo) - expect_true(x2$is.same.haplo) -}) diff --git a/tests/testthat/test-preferential_pairing.R b/tests/testthat/test-preferential_pairing.R deleted file mode 100644 index a14f5eda..00000000 --- a/tests/testthat/test-preferential_pairing.R +++ /dev/null @@ -1,11 +0,0 @@ -context("Preferential Pairing") -test_that("computes pairing probabilities correctly", { - x1 <- get_submap(solcap.dose.map[[1]], 1:20, reestimate.rf = F) - probs.t1<-calc_genoprob(input.map = x1, - verbose = TRUE) - expect_equivalent(round(var(probs.t1$probs), 6), 0.015261) - pref.t1 <- calc_prefpair_profiles(input.genoprobs = probs.t1) - expect_equivalent(round(apply(pref.t1$prefpair.psi.pval, 2, mean),6), - c(0.822695, 0.876216, 1.000000, 7.917049)) - expect_is(plot(pref.t1), "ggplot") -}) \ No newline at end of file diff --git a/tests/testthat/test-read_data.R b/tests/testthat/test-read_data.R deleted file mode 100644 index 77c79da0..00000000 --- a/tests/testthat/test-read_data.R +++ /dev/null @@ -1,40 +0,0 @@ -context("Read data") -test_that("read data from VCF file correctly", { - fl = "https://github.com/mmollina/MAPpoly_vignettes/raw/master/data/sweet_sample_ch3.vcf.gz" - tempfl <- tempfile(pattern = 'chr1_', fileext = '.vcf.gz') - download.file(fl, destfile = tempfl) - dat.dose.vcf = read_vcf(file = tempfl, parent.1 = "PARENT1", parent.2 = "PARENT2") - expect_equal(check_data_sanity(dat.dose.vcf), 0) - expect_null(print(dat.dose.vcf, detailed = TRUE)) -}) -test_that("read data from dosage file correctly", { - fl1 = "https://raw.githubusercontent.com/mmollina/MAPpoly_vignettes/master/data/SolCAP_dosage" - tempfl <- tempfile() - download.file(fl1, destfile = tempfl) - SolCAP.dose <- read_geno(file.in = tempfl) - expect_equal(check_data_sanity(SolCAP.dose), 0) -}) -test_that("read data from probability file correctly", { - ft="https://raw.githubusercontent.com/mmollina/MAPpoly_vignettes/master/data/hexa_sample" - tempfl <- tempfile() - download.file(ft, destfile = tempfl) - SolCAP.dose.prob <- read_geno_prob(file.in = tempfl) - expect_equal(check_data_sanity(SolCAP.dose.prob), 0) -}) -test_that("read data from CSV file correctly", { - ft="https://raw.githubusercontent.com/mmollina/MAPpoly_vignettes/master/data/tetra_solcap.csv" - tempfl <- tempfile() - download.file(ft, destfile = tempfl) - SolCAP.dose <- read_geno_csv(file.in = tempfl, ploidy = 4) - expect_equal(check_data_sanity(SolCAP.dose), 0) -}) -test_that("read data from fitpoly file correctly", { - fl <- "https://raw.githubusercontent.com/mmollina/MAPpoly_vignettes/master/data/fitpoly.dat" - tempfl <- tempfile() - download.file(fl, destfile = tempfl) - fitpoly.dat <- read_fitpoly(file.in = tempfl, ploidy = 4, - parent1 = "P1", parent2 = "P2", - verbose = TRUE) - expect_equal(check_data_sanity(fitpoly.dat), 0) - expect_null(print(fitpoly.dat, detailed = TRUE)) -}) diff --git a/tests/testthat/test-reestmap.R b/tests/testthat/test-reestmap.R deleted file mode 100644 index 0cf170bf..00000000 --- a/tests/testthat/test-reestmap.R +++ /dev/null @@ -1,17 +0,0 @@ -context("Map reestimation") -test_that("reestimate genetic maps correctly correctly", { - x1 <- get_submap(solcap.dose.map[[1]], 1:5, reestimate.rf = F) - s<-make_seq_mappoly(x1) - tpt<-est_pairwise_rf(s) - m <- rf_list_to_matrix(tpt) - x2 <- reest_rf(x1) - expect_equivalent(mean(x2$maps[[1]]$seq.rf), 0.00973, tolerance = 1e-4) - x3 <- reest_rf(x1, input.mat = m, method = "ols") - expect_equivalent(mean(x3$maps[[1]]$seq.rf), 0.0248, tolerance = 1e-4) - x4 <- est_full_hmm_with_global_error(x1) - expect_equivalent(mean(x4$maps[[1]]$seq.rf), 0.00993, tolerance = 1e-4) - expect_error(x5 <- est_full_hmm_with_prior_prob(x1)) - x5 <- get_submap(solcap.prior.map[[1]], 1:10, reestimate.rf = FALSE) - x6 <- est_full_hmm_with_prior_prob(x5, tol = 10e-3) - expect_equivalent(mean(x6$maps[[1]]$seq.rf), 0.008775, tolerance = 1e-4) -}) \ No newline at end of file diff --git a/tests/testthat/test-simulation.R b/tests/testthat/test-simulation.R deleted file mode 100644 index 49a0873a..00000000 --- a/tests/testthat/test-simulation.R +++ /dev/null @@ -1,6 +0,0 @@ -context("Simulate datasets") -test_that("simulate datasets correctly", { - h.temp<-sim_homologous(ploidy=6, n.mrk=20) - dat<-cross_simulate(h.temp, 100, n.ind = 100) - expect_is(dat, "mappoly.data") -}) diff --git a/tests/testthat/test-split_and_rephase.R b/tests/testthat/test-split_and_rephase.R deleted file mode 100644 index e3f3be02..00000000 --- a/tests/testthat/test-split_and_rephase.R +++ /dev/null @@ -1,8 +0,0 @@ -context("Split and rephase") -test_that("split and rephase the map correctly", { - ##### Tetraploid - map<-get_submap(solcap.err.map[[1]], 1:20) - tpt<-est_pairwise_rf(make_seq_mappoly(map)) - map2<-split_and_rephase(map, tpt, gap.threshold = 2) - expect_is(map2, "mappoly.map") -}) diff --git a/tests/testthat/test-utility_func.R b/tests/testthat/test-utility_func.R deleted file mode 100644 index 1d57b31d..00000000 --- a/tests/testthat/test-utility_func.R +++ /dev/null @@ -1,90 +0,0 @@ -context("Utility functions") -test_that("test several utility functions", { - x1<-plot_mrk_info(tetra.solcap, 1) - x2<-plot_mrk_info(tetra.solcap, "solcap_snp_c2_41437") - x3<-plot_mrk_info(tetra.solcap.geno.dist, 1) - x4<-plot_mrk_info(tetra.solcap.geno.dist, "solcap_snp_c2_41437") - expect_is(x1, "list") - expect_is(x2, "list") - expect_is(x3, "list") - expect_is(x4, "list") - expect_equal(get_LOD(solcap.err.map[[1]]),0) - tpt<-est_pairwise_rf(make_seq_mappoly(tetra.solcap, 1:30)) - x1<-make_seq_mappoly(tetra.solcap, 1:20) - x2<-make_seq_mappoly(tetra.solcap, 21:40) - x1<-as.numeric(check_if_rf_is_possible(x1)) - x2<-as.numeric(check_if_rf_is_possible(x2)) - expect_equal(as.numeric(crossprod(x1,x2)), 14) - M<-rf_list_to_matrix(tpt, shared.alleles = TRUE) - expect_equal(round(sum(get_rf_from_mat(M$rec.mat), na.rm = TRUE), 6), 3.913633) - expect_equal(get_w_m(6), 15) - expect_error(get_w_m(0)) - expect_error(get_w_m(3)) - rm<-rev_map(maps.hexafake[[1]]) - expect_equal(sum(rm$maps[[1]]$seq.rf), sum(maps.hexafake[[1]]$maps[[1]]$seq.rf), tolerance = 10e-5) - a1<-sample_data(tetra.solcap.geno.dist, n = 20, type = "marker") - a2<-sample_data(a1, n = 20, type = "individual") - plot(a2) - a3<-dist_prob_to_class(a2$geno) - expect_equal(prod(dim(a3$geno.dose)), nrow(a3$geno)) - #w1<-export_data_to_polymapR(hexafake) - #expect_equal(dim(w1) , c(1500, 302)) - #expect_equal(sum(w1), 461530) - #w2<-export_data_to_polymapR(tetra.solcap) - #expect_equal(dim(w2) , c(3679, 162)) - #expect_equal(sum(w2, na.rm = TRUE), 1058080) - expect_equal(round(mf_h(20),6), 0.16484) - expect_equal(round(mf_k(20),6), 0.189974) - expect_equal(round(mf_m(20),6), 0.2) - expect_equal(round(imf_h(.15),2), 17.83) - expect_equal(round(imf_k(.15),2), 15.48) - expect_equal(round(imf_m(.15),2), 15) - expect_null(plot_compare_haplotypes(ploidy = 6, - hom.allele.p1 = maps.hexafake[[1]]$maps[[1]]$seq.ph$P[1:10], - hom.allele.q1 = maps.hexafake[[1]]$maps[[1]]$seq.ph$Q[1:10], - hom.allele.p2 = maps.hexafake[[1]]$maps[[1]]$seq.ph$P[1:10], - hom.allele.q2 = maps.hexafake[[1]]$maps[[1]]$seq.ph$Q[1:10])) - a <- print_mrk(input.data = hexafake, mrks = "M_1") - expect_true(is.list(a)) - expect_equal(sum(unlist(a)), expected = 304) - expect_equal(sum(perm_pars(1:5)), 900) - expect_equal(sum(perm_tot(1:5)), 1800) - a<-sample_data(tetra.solcap.geno.dist, n=30, type = "marker") - a<-sample_data(a, n=30) - w<-update_missing(a, prob.thres = .7) - expect_is(w, "mappoly.data") - w2 <- get_genomic_order(make_seq_mappoly(hexafake, "all")) - w3 <- as.numeric(crossprod(w2$ord$seq.pos)/10e17) - expect_is(plot(w2), "ggplot") - expect_equal(w3, 0.4952, tolerance = 1e-3) - ##test drop - w4<-get_submap(solcap.dose.map[[1]], 1:10, reestimate.rf = FALSE) - s4<-make_seq_mappoly(w4) - tpt2<-est_pairwise_rf(s4) - M2<-rf_list_to_matrix(tpt2, shared.alleles = TRUE) - a1 <- drop_marker(w4, 5) - a2 <- drop_marker(w4, "solcap_snp_c1_10930") - expect_equal(a1,a2, tolerance = 1e-3) - ##test add - a3 <- drop_marker(w4, 1) - a4 <- add_marker(a3, "solcap_snp_c2_51460", 0, M2, verbose = FALSE) - w4<-reest_rf(w4, tol = 10e-5) - a4<-reest_rf(a4, tol = 10e-5) - expect_equal(w4, a4, tolerance = 1e-3) - a3 <- drop_marker(w4, 5) - a4 <- add_marker(a3, "solcap_snp_c1_10930", 4, M2, verbose = FALSE) - w4<-reest_rf(w4, tol = 10e-5) - a4<-reest_rf(a4, tol = 10e-5) - expect_equal(w4, a4, tolerance = 1e-3) - a3 <- drop_marker(w4, 10) - a4 <- add_marker(a3, "solcap_snp_c2_36643", 9, M2, verbose = FALSE) - w4<-reest_rf(w4, tol = 10e-5) - a4<-reest_rf(a4, tol = 10e-5) - expect_equal(w4, a4, tolerance = 1e-3) - ##merge data - b<-merge_datasets(hexafake, hexafake.geno.dist) - expect_equivalent(mean(as.matrix(b$geno.dose)), 1.078817, tol = 1e-3) - ##update map - map<-update_map(solcap.dose.map[[1]]) - expect_equal(length(map$maps[[1]]$seq.num) - length(solcap.dose.map[[1]]$maps[[1]]$seq.num), 20) -})