Skip to content

Commit

Permalink
Updating documentation and eliminating C++11 requirement
Browse files Browse the repository at this point in the history
  • Loading branch information
mmollina committed Nov 24, 2023
1 parent 1ddbf87 commit 637c5a1
Show file tree
Hide file tree
Showing 5 changed files with 56 additions and 51 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ Imports:
plotly
LinkingTo: Rcpp, RcppParallel
RoxygenNote: 7.2.3
SystemRequirements: C++11, GNU make
SystemRequirements: GNU make
Encoding: UTF-8
Suggests:
testthat,
Expand Down
41 changes: 22 additions & 19 deletions R/plot_progeny_dosage_change.R
Original file line number Diff line number Diff line change
@@ -1,32 +1,34 @@
#' Look at genotypes that were imputed or changed by the HMM chain given a level of global genotypic error
#' Display genotypes imputed or changed by the HMM chain given a global genotypic error
#'
#' Outputs a graphical representation ggplot with the percent of data changed.
#'
#' Most recent update 8/29/2023:
#' -fixed issue where only worked on tetraploid to now working for diploid to octaploid.
#' -un-hardcoded linkage groups in maps. previously hard-coded for tetraploid rose.
#'
#'
#' @param map_list a list of multiple \code{mappoly.map.list}
#'
#' @param error error rate used in global error in the `calc_genoprob_error()`
#' @param verbose T or F for `calc_genoprob_error()` and `calc_homologprob()`
#' @param output_corrected TRUE or FALSE, if F only the ggplot of the changed dosage is printed, if T then a new corrected dosage matrix is output.
#' @param verbose if TRUE (default), current progress is shown; if FALSE,
#' no output is produced
#' @param output_corrected logical. if FALSE only the ggplot of the changed
#' dosage is printed, if TRUE then a new corrected dosage matrix is output.
#'
#' @return A ggplot of the changed and imputed genotypic dosages
#'
#' @examples
#' x <- get_submap(solcap.err.map[[1]], 1:30, reestimate.rf = FALSE) # subset to run fast in tests
#' plot_progeny_dosage_change(list(x), error=0.05, output_corrected=F) # no new corrected dosages
#' corrected_matrix <- plot_progeny_dosage_change(list(x), error=0.05, output_corrected=F) #output corrected
#' x <- get_submap(solcap.err.map[[1]], 1:30, reestimate.rf = FALSE)
#' plot_progeny_dosage_change(list(x), error=0.05, output_corrected=FALSE)
#' corrected_matrix <- plot_progeny_dosage_change(list(x), error=0.05,
#' output_corrected=FALSE) #output corrected
#'
#' @author Jeekin Lau, \email{[email protected]}, with optimization by Cristiane Taniguti, \email{[email protected]}
#' @author Jeekin Lau, \email{[email protected]}, with optimization by Cristiane
#' Taniguti, \email{[email protected]}
#'
#' @import ggplot2
#' @import reshape2
#' @export plot_progeny_dosage_change

plot_progeny_dosage_change <- function(map_list, error, verbose=T, output_corrected=F){
plot_progeny_dosage_change <- function(map_list,
error,
verbose = TRUE,
output_corrected = FALSE){
Var1 <- Var2 <- value <- NULL
map=map_list
if(!exists(map[[1]]$info$data.name)) stop("mappoly.data object not here")
Expand All @@ -38,7 +40,8 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T, output_correc

genoprob <- vector("list", length(map))
for(i in 1:length(map)){
genoprob[[i]] <- calc_genoprob_error(input.map = map[[i]], error = error, verbose = verbose)
genoprob[[i]] <- calc_genoprob_error(input.map = map[[i]], error = error,
verbose = verbose)
}

print("calculating homologprob")
Expand All @@ -47,8 +50,8 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T, output_correc

print("comparing to orginal")

P=unlist(lapply(map, function(x) x$maps[[1]]$seq.ph$P), recursive=F)
Q=unlist(lapply(map, function(x) x$maps[[1]]$seq.ph$Q), recursive=F)
P=unlist(lapply(map, function(x) x$maps[[1]]$seq.ph$P), recursive=FALSE)
Q=unlist(lapply(map, function(x) x$maps[[1]]$seq.ph$Q), recursive=FALSE)



Expand Down Expand Up @@ -92,7 +95,7 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T, output_correc

test <- sapply(homoprob_ind, function(x) {
by_marker <- split.data.frame(x, x[,1])
final_matrix <- t(sapply(by_marker, function(y) y[order(y[,4], decreasing = T),][1:ploidy,2]))
final_matrix <- t(sapply(by_marker, function(y) y[order(y[,4], decreasing = TRUE),][1:ploidy,2]))
final_vector <- apply(final_matrix, 1, function(w) paste0(w, collapse = ""))
return(final_vector)
})
Expand Down Expand Up @@ -169,15 +172,15 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T, output_correc
empty_matrix[which(!original_geno==finished&!original_geno==ploidy+1)]="changed"
empty_matrix_melt=melt(empty_matrix)
plot1<-ggplot(empty_matrix_melt, aes(Var1, Var2, fill= factor(value, levels=c("changed","imputed","unchanged")))) +
geom_tile()+scale_fill_manual(values=colors, drop=F)+
geom_tile()+scale_fill_manual(values=colors, drop=FALSE)+
xlab("Markers")+
ylab("Individuals")+
ggtitle(paste0("changed = ",round(percent_changed, digits=3),"% ","imputed = ", round(percent_imputed, digits=3),"%"))+
guides(fill=guide_legend(title="Dosage change"))
print("done")

print(plot1)
if (output_corrected ==T){
if (output_corrected ==TRUE){
mrk_names=rownames(finished)
mrk_index=which(dat$mrk.names%in%mrk_names)
P1 = dat$dosage.p1[mrk_index]
Expand Down
30 changes: 17 additions & 13 deletions man/plot_progeny_dosage_change.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 0 additions & 2 deletions src/Makevars.win
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
PKG_LIBS = -lz

CXX_STD = CXX11

PKG_CXXFLAGS += -DRCPP_PARALLEL_USE_TBB=1

PKG_LIBS += $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" \
Expand Down
32 changes: 16 additions & 16 deletions src/RcppExports.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -59,22 +59,22 @@ BEGIN_RCPP
END_RCPP
}

RcppExport SEXP calc_genoprob(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
RcppExport SEXP calc_genoprob_prior(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
RcppExport SEXP calc_genprob_haplo(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
RcppExport SEXP calc_genprob_haplo_highprec(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 loglike_hmm(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, 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);
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 *, 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 *);

static const R_CallMethodDef CallEntries[] = {
{"_mappoly_vcf_get_probabilities", (DL_FUNC) &_mappoly_vcf_get_probabilities, 2},
Expand Down

0 comments on commit 637c5a1

Please sign in to comment.