V1.1 features

DESCRIPTION

@@ -66,7 +66,8 @@ Imports:
     Matrix,
     matrixcalc,
     markdown,
-    Rsamtools
+    Rsamtools,
+    methods
 Remotes:
     github::jendelman/GWASpoly,
     github::Breeding-Insight/BIGr,

NAMESPACE

@@ -84,6 +84,7 @@ importFrom(httr,GET)
 importFrom(httr,content)
 importFrom(httr,status_code)
 importFrom(matrixcalc,is.positive.definite)
+importFrom(methods,new)
 importFrom(plotly,add_markers)
 importFrom(plotly,ggplotly)
 importFrom(plotly,layout)

R/DosageCall_functions.R

@@ -8,14 +8,20 @@
 #' @param backcross.gen ToDo
 #' @param intercross.gen ToDo
 #' @param selfing.gen ToDo
+#' @param min_ind_read ToDo
+#' @param min_ind_maf ToDo
 #' @param contamRate ToDo
+#' @param tol ToDo
 #' @param session ToDo
 #'
 #' @import polyRAD
 #' @importFrom vcfR read.vcfR is.biallelic write.vcf
 #' @importFrom R.utils gunzip
 #'
-polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL, backcross.gen = 0, intercross.gen = 0, selfing.gen = 0, contamRate = 0.001, session) {
+polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL,
+                                backcross.gen = 0, intercross.gen = 0, selfing.gen = 0,
+                                contamRate = 0.001, min_ind_read = 1, min_ind_maf = 0,
+                                tol = 1e-05, session) {
 
   # Variables
   vcf_path <- vcf
@@ -28,36 +34,36 @@ polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL, backcr
   # Having some issues formatting the output when multiallelic SNPs is input, so excluding for now
   temp_vcf <- vcfR::read.vcfR(vcf_path, verbose = FALSE)
   temp_vcf <- temp_vcf[is.biallelic(temp_vcf),]
-  
+
   # Function to randomly assign non-matching bases for polyRAD
   assign_random_bases <- function(vcf) {
     # Get the number of rows in the VCF
     num_rows <- nrow(vcf@fix)
-    
+
     # Generate random bases for REF and ALT
     bases <- c("A", "C", "G", "T")
     ref_bases <- sample(bases, num_rows, replace = TRUE)
     alt_bases <- character(num_rows) # Initialize an empty character vector
-    
+
     # Ensure REF and ALT are different for each row
     for (i in 1:num_rows) {
       alt_bases[i] <- sample(bases[bases != ref_bases[i]], 1)
     }
-    
+
     # Assign the new bases to the vcf object
     vcf@fix[, "REF"] <- ref_bases
     vcf@fix[, "ALT"] <- alt_bases
-    
+
     return(vcf)
   }
-  
+
   #Editing REF and ALT fields randomly temporarily if blank or AB
   ref_base <- temp_vcf@fix[1, "REF"]
   alt_base <- temp_vcf@fix[1, "ALT"]
   if (is.na(ref_base) || is.na(alt_base) || alt_base == "B") {
   temp_vcf <- assign_random_bases(temp_vcf)
   }
-  
+
   # Adding filtered VCF as a temp object
   temp_vcf_path <- tempfile(fileext = ".vcf.gz")
   vcfR::write.vcf(temp_vcf, file = temp_vcf_path)
@@ -73,10 +79,11 @@ polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL, backcr
   # Retaining all markers with at least 1 individual with reads (minimal filtering at this stage)
   # Need to determine the best contamRate to use; currently using the default
   polyRAD_obj <- polyRAD::VCF2RADdata(file = temp_unzipped_path,
-                                      min.ind.with.reads = 1,
-                                      min.ind.with.minor.allele = 0,
+                                      min.ind.with.reads = min_ind_read,
+                                      min.ind.with.minor.allele = min_ind_maf,
                                       taxaPloidy = ploidy,
                                       contamRate = contamRate,
+                                      tol = tol,
                                       phaseSNPs = FALSE
 
   )
@@ -106,13 +113,13 @@ polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL, backcr
   # Perform dosage calling
   # "If you expect that your species has high linkage disequilibrium, the functions IterateHWE_LD and IteratePopStructLD behave like IterateHWE and IteratePopStruct, respectively, but also update priors based on genotypes at linked loci."
   if (model == "IterateHWE") {
-    rad <- IterateHWE(polyRAD_obj, tol = 1e-5, overdispersion = OD)
+    rad <- IterateHWE(polyRAD_obj, tol = tol, overdispersion = OD)
   } else if (model == "IteratePopStruct") {
-    rad <- IteratePopStruct(polyRAD_obj, tol = 1e-5, overdispersion = OD)
+    rad <- IteratePopStruct(polyRAD_obj, tol = tol, overdispersion = OD)
   } else if (model == "IterateHWE_LD") {
-    rad <- IterateHWE_LD(polyRAD_obj, tol = 1e-5, overdispersion = OD)
+    rad <- IterateHWE_LD(polyRAD_obj, tol = tol, overdispersion = OD)
   } else if (model == "IteratePopStruct_LD") {
-    rad <- IteratePopStruct_LD(polyRAD_obj, tol = 1e-5, overdispersion = OD)
+    rad <- IteratePopStruct_LD(polyRAD_obj, tol = tol, overdispersion = OD)
   } else {
     if (!is.null(p1)) {
       # First check that p1 is a valid parent
@@ -135,7 +142,8 @@ polyRAD_dosage_call <- function(vcf, ploidy, model, p1 = NULL, p2 = NULL, backcr
                                    overdispersion = OD,
                                    n.gen.backcrossing = n.bx,
                                    n.gen.intermating = n.inter,
-                                   n.gen.selfing = n.self)
+                                   n.gen.selfing = n.self,
+                                   tol=tol)
   }
 
 

R/mod_DosageCall.R

@@ -91,7 +91,7 @@ mod_DosageCall_ui <- function(id){
                 )
               )
             ),
-            numericInput(ns("cores"), "Number of CPU Cores*", min = 1, max = (future::availableCores() - 1), value = 1),
+            sliderInput(ns("cores"), "Number of CPU Cores*", min = 1, max = (availableCores() - 1), value = 1, step = 1),
           ),
           conditionalPanel(
             condition = "input.Rpackage == 'polyRAD'",
@@ -148,13 +148,21 @@ mod_DosageCall_ui <- function(id){
             status = "warning",
             icon = icon("info"), width = "500px",
             tooltip = tooltipOptions(title = "Click to see info!")
-          ))
+          )),
+          br(),
+          tags$hr(style="border-color: #d3d3d3; margin-top: 20px; margin-bottom: 20px;"),  # Lighter grey line
+          div(style="text-align: left; margin-top: 10px;",
+              actionButton(ns("advanced_options"),
+                           label = HTML(paste(icon("cog", style = "color: #007bff;"), "Advanced Options")),
+                           style = "background-color: transparent; border: none; color: #007bff; font-size: smaller; text-decoration: underline; padding: 0;"
+              )
+          )
         ),
         column(width=4,
-          valueBoxOutput(ns("MADCsnps"), width=12),
-          box(title = "Status", width = 12, collapsible = TRUE, status = "info",
-            progressBar(id = ns("pb_madc"), value = 0, status = "info", display_pct = TRUE, striped = TRUE, title = " ")
-          )
+               valueBoxOutput(ns("MADCsnps"), width=12),
+               box(title = "Status", width = 12, collapsible = TRUE, status = "info",
+                   progressBar(id = ns("pb_madc"), value = 0, status = "info", display_pct = TRUE, striped = TRUE, title = " ")
+               )
         )
       )
     )
@@ -249,6 +257,73 @@ mod_DosageCall_server <- function(input, output, session, parent_session){
     valueBox(snp_number(), "Markers in uploaded file", icon = icon("dna"), color = "info")
   })
 
+  #Default model choices
+  advanced_options <- reactiveValues(
+    contamRate = 0.001,
+    min_ind_read = 1,
+    min_ind_maf = 0,
+    tol = 1e-05
+  )
+
+  #UI popup window for input
+  observeEvent(input$advanced_options, {
+    showModal(modalDialog(
+      title = "Advanced Options",
+      conditionalPanel(
+        condition = "input.Rpackage == 'polyRAD'", ns = ns,
+        div(
+          numericInput(
+            inputId = ns('contamRate'),
+            label = 'contamRate',
+            min = 0,
+            max = 1,
+            value = advanced_options$contamRate
+          ),
+          numericInput(
+            inputId = ns('min_ind_read'),
+            label = 'min.ind.with.reads',
+            min = 1,
+            value = advanced_options$min_ind_read
+          ),
+          numericInput(
+            inputId = ns('min_ind_maf'),
+            label = 'min.ind.with.minor.allele',
+            min = 0,
+            value = advanced_options$min_ind_maf
+          ),
+          numericInput(
+            inputId = ns('tol'),
+            label = 'tol',
+            min = 0,
+            max = 1,
+            value = advanced_options$tol
+          )
+        )
+      ),
+      conditionalPanel(
+        condition = "input.Rpackage == 'Updog'", ns = ns,
+        "Currently, there are no Advanced Options for Updog"
+      ),
+      footer = tagList(
+        modalButton("Close"),
+        actionButton(ns("save_advanced_options"), "Save")
+      )
+    ))
+  })
+
+
+
+  #Close popup window when user "saves options"
+  observeEvent(input$save_advanced_options, {
+    advanced_options$contamRate <- input$contamRate
+    advanced_options$min_ind_read <- input$min_ind_read
+    advanced_options$min_ind_maf <- input$min_ind_maf
+    advanced_options$tol <- input$tol
+    # Save other inputs as needed
+
+    removeModal()  # Close the modal after saving
+  })
+
   disable("download_updog_vcf")
 
   ##This is for performing Dosage Calling
@@ -461,7 +536,11 @@ mod_DosageCall_server <- function(input, output, session, parent_session){
                                            p2 = input$parent2,
                                            backcross.gen = input$bx.gen,
                                            intercross.gen = input$inter.gen,
-                                           selfing.gen = input$self.gen)
+                                           selfing.gen = input$self.gen,
+                                           contamRate = advanced_options$contamRate,
+                                           min_ind_read = advanced_options$min_ind_read,
+                                           min_ind_maf = advanced_options$min_ind_maf,
+                                           tol=advanced_options$tol)
       updateProgressBar(session = session, id = "pb_madc", value = 100, title = "Finished")
       polyrad_items
     }
@@ -476,12 +555,7 @@ mod_DosageCall_server <- function(input, output, session, parent_session){
   output$download_updog_vcf <- downloadHandler(
     filename = function() {
       output_name <- gsub("\\.vcf$", "", input$output_name)
-      if (input$Rpackage == "Updog") {
-        paste0(output_name, ".vcf.gz")
-      }else {
-        paste0(output_name, ".vcf")
-      }
-
+      paste0(output_name, ".vcf.gz")
     },
     content = function(file) {
 
@@ -492,24 +566,19 @@ mod_DosageCall_server <- function(input, output, session, parent_session){
           multidog.object = updog_out(),
           output.file = temp,
           updog_version = packageVersion("updog"),
-          compress = TRUE
+          compress = FALSE
         )
 
-        # Move the file to the path specified by 'file'
-        file.copy(paste0(temp, ".vcf.gz"), file)
-
       } else {
         polyRAD2vcf(updog_out()$Genos,
                     model = input$polyRAD_model,
                     vcf_path = input$madc_file$datapath,
                     hindhe.obj = updog_out()$RADHindHe,
                     ploidy = input$ploidy,
                     output.file = temp
-                    )
-
-        # Move the file to the path specified by 'file'
-        file.copy(paste0(temp, ".vcf"), file)
+        )
       }
+      bgzip_compress(paste0(temp, ".vcf"), file)
 
       # Delete the temporary file
       unlink(temp)