Merge pull request #3 from mmollina/main

Update

DESCRIPTION

@@ -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,

R/homolog_probs.R

@@ -55,8 +55,23 @@ calc_homologprob <- function(input.genoprobs, verbose = TRUE){
     ploidy <- v[as.character(length(stt.names))]
     hom.prob <- array(NA, dim = c(ploidy*2, length(mrk.names), length(ind.names)))
     dimnames(hom.prob) <- list(letters[1:(2*ploidy)], mrk.names, ind.names)
-    for(i in letters[1:(2*ploidy)])
-      hom.prob[i,,] <- apply(input.genoprobs[[j]]$probs[grep(stt.names, pattern = i),,], c(2,3), function(x) round(sum(x, na.rm = TRUE),4))
+    ## for(i in letters[1:(2*ploidy)])
+    ##   hom.prob[i,,] <- apply(input.genoprobs[[j]]$probs[grep(stt.names, pattern = i),,], c(2,3), function(x) round(sum(x, na.rm = TRUE),4))
+    ## Incidence matrix for homologs
+    alleles = matrix(unlist(strsplit(stt.names, '')), ncol=(ploidy+1), byrow=TRUE)[,-c((ploidy/2)+1)]
+    ## Creating incidence matrix to transition from genotypes to homologs
+    alleles2 = matrix(0, length(stt.names), ploidy*2)
+    ## Creating dummy variables to associate genotypes with alleles
+    for(i in 1:nrow(alleles)){
+      for(k in 1:ncol(alleles)){
+        alleles2[i,match(alleles[i,k], letters)]=1.00
+      }
+    }
+    ## Getting homolog probabilities
+    for(m in 1:length(mrk.names)) {
+      hom.prob[,m,] <- round(t(alleles2)%*%input.genoprobs[[j]]$probs[,m,],4)
+      ## hom.prob[,m,] <- round(t(alleles2)%*%input.genoprobs[[j]]$probs[,m,], 4)
+    }
     df.hom <- reshape2::melt(hom.prob)
     map <- data.frame(map.position = input.genoprobs[[j]]$map, marker = names(input.genoprobs[[j]]$map))
     colnames(df.hom) <- c("homolog", "marker", "individual", "probability")

R/plot_progeny_dosage_change.R

@@ -1,31 +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 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)
-#'     plot_progeny_dosage_change(list(x), error=0.05)     
+#'       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){
+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")
@@ -37,7 +40,8 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T){
 
   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")
@@ -46,8 +50,8 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T){
 
   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)
 
 
 
@@ -91,7 +95,7 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T){
 
   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)
   })
@@ -168,15 +172,24 @@ plot_progeny_dosage_change <- function(map_list, error, verbose=T){
   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")
 
-  return(plot1) 
-
+  print(plot1)
+  if (output_corrected ==TRUE){
+  mrk_names=rownames(finished)
+  mrk_index=which(dat$mrk.names%in%mrk_names)
+  P1 = dat$dosage.p1[mrk_index]
+  P2 = dat$dosage.p2[mrk_index]
+  sequence = dat$chrom[mrk_index]
+  sequence_position = dat$genome.pos[mrk_index]
+
+  hmm_imputed = cbind(P1,P2,sequence,sequence_position, finished)
+  return(hmm_imputed)}
 }
 
 

src/Makevars.win

@@ -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" \