Implement support for custom SNP sets

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mixder
 Type: Package
 Title: A workflow for performing SNP mixture deconvolution
-Version: 0.9.0
+Version: 0.10.0
 Author: Rebecca Mitchell <rebecca.mitchell@st.dhs.gov>
 Maintainer: Rebecca Mitchell <rebecca.mitchell@st.dhs.gov>
 Description: A workflow for performing SNP mixture deconvolution of ForenSeq Kintelligence SNPs. Mixture deconvolution is performed using EuroForMix (https://github.com/oyvble/euroformix/). After mixture deconvolution, the user can choose to calculate metrics such as genotype accuracy and heterozygosity for a range of allele probability thresholds (useful for validation work) or create GEDmatch PRO reports utilizing specified allele probability thresholds.

NAMESPACE

@@ -13,6 +13,7 @@ export(check_ratios)
 export(check_reads)
 export(checking_af)
 export(compile_metrics)
+export(convert_table_to_list)
 export(create_at)
 export(create_config)
 export(create_evid)
@@ -33,6 +34,7 @@ export(process_efm_files)
 export(process_kinreport)
 export(processing_evid_sample_reports)
 export(processing_ref_sample_reports)
+export(read_in_freq)
 export(read_in_table)
 export(reverse_comp)
 export(run_efm)
@@ -49,8 +51,10 @@ import(kgp)
 import(parallel)
 import(plotly)
 import(prompter)
+import(rrapply)
 import(shiny)
 import(shinyFiles)
+importFrom(data.table,fread)
 importFrom(grDevices,dev.off)
 importFrom(grDevices,png)
 importFrom(methods,show)

NEWS.md

@@ -1,3 +1,15 @@
+# mixder 0.10.0
+
+## Changes in version 0.10.0
+- Created more flexibility with file format for sample manifest
+- Implemented support for custom SNP panels:
+  - In GUI, Drop down menu to select either "Kintelligence" or "Custom" assay; will require SNP positions file if using custom
+  - Using R API, `assay` and `snp_pos` arguments used to specify if a custom panel is used
+  - Wrote functions to read in and create EFM-compatible files more efficiently
+  - Save large data files (i.e. reference SNP profiles, frequency files, etc.) as `.rda` files and will auto-detect these to load instead of reading in large file (significantly faster)
+- Load frequency data only once instead of with each sample when using Shiny app
+
+
 # mixder 0.9.0
 
 ## Changes in version 0.9.0

R/ancestry_prediction.R

@@ -28,12 +28,12 @@
 ancestry_prediction = function(report, path, id, analysis_type, contrib_status, testsnps, groups) {
   if (testsnps == "All Autosomal SNPs") {
     plotid="AllSNPs"
-    geno=mixder::ancestry_1000G_allsamples
+    geno=data.frame(mixder::ancestry_1000G_allsamples)
   } else {
     plotid="AncestrySNPsOnly"
     geno=mixder::ancestrysnps_1000G_allsamples
   }
-  ncols=ncol(geno)
+  ncols = ncol(geno)
   geno_filt=geno[,c(7:ncols)]
   snps = data.frame("snp_id"=colnames(geno_filt))
   snps = snps %>%
@@ -69,7 +69,6 @@ ancestry_prediction = function(report, path, id, analysis_type, contrib_status,
   newcol=ncols+1
   newcol2=ncols+4
   PCs_anc = cbind(geno_ancestry[,c(newcol:newcol2)], data.frame(PCs[,c(1:10)]))
-
   centroids(groups, PCs_anc, glue("{path}/PCA_plots"), glue("{id}_{contrib_status}_{analysis_type}_{plotid}"))
 
   dir.create(file.path(path, "PCA_plots"), showWarnings = FALSE, recursive=TRUE)

R/app.R

@@ -66,7 +66,8 @@ mixder = function() {
       textOutput("kin_inpath"),
       conditionalPanel(condition = "input.method == 'Calculate Metrics' | input.cond == 1", uiOutput("ref_GetFile"), uiOutput("ref_text")),
       conditionalPanel(condition = "input.cond == 1", uiOutput("ref_selector")),
-      conditionalPanel(condition = "input.skip_ancestry == 1", uiOutput("runmd"), uiOutput("twofreqs"), uiOutput("method")),
+      conditionalPanel(condition = "input.skip_ancestry == 1", uiOutput("runmd"), uiOutput("twofreqs"), uiOutput("method"), uiOutput("assay_choice")),
+      conditionalPanel(condition = "input.assay_choice == 'Custom'", uiOutput("assay_GetFile"), uiOutput("assay_text")),
       conditionalPanel(condition = "input.skip_ancestry == 1 & input.twofreqs == 0", uiOutput("freqselect")),
       conditionalPanel(condition = "input.twofreqs == 1", uiOutput("freqselect_major"), uiOutput("freqselect_minor")),
       conditionalPanel(condition = "input.skip_ancestry == 1 & input.twofreqs == 0 & input.uploadfreq == 'Upload Custom'", uiOutput("freq_GetFile"), uiOutput("freq_text")),
@@ -148,6 +149,29 @@ server = function(input, output, session) {
     add_prompt(message = "Optional- if selecting Calculate Metrics,\nmust include reference genotypes.", position = "right")
   ), c("", "Calculate Metrics", "Create GEDmatch PRO Report"))
   })
+  output$assay_choice = renderUI({
+    selectInput("assay_choice", tags$span("Select Assay", tags$span(icon(
+      name = "question-circle",
+    )
+    ) |>
+      add_prompt(message = "Select assay used to develop genotypes", position = "right")
+    ), c("Kintelligence", "Custom"))
+  })
+  output$assay_GetFile = renderUI({
+    fluidRow(column(10,shinyFilesButton("assay_GetFile", "Select a custom SNP positions file" ,
+                                        title = "Select a custom SNP positions file", multiple = FALSE,
+                                        buttonType = "default", class = NULL),
+                    tags$span(
+                      icon(
+                        name = "question-circle",
+                      )
+                    ) |>
+                      add_prompt(message = "File containing chromosomal locations for custom SNP panels.\nThis is not required if using the Kintelligence assay.\nColumns required: rsID, chromosome (e.g. 1, 2, 3,), and position (GRCh37).", position = "right")
+    ))
+  })
+  output$assay_text = renderUI({
+    textOutput("assay_file")
+  })
   output$twofreqs = renderUI({
     checkboxInput("twofreqs", tags$span("Use Different Allele Frequency Files For Each Contributor?", tags$span(icon(
     name = "question-circle",
@@ -371,52 +395,60 @@ server = function(input, output, session) {
       add_prompt(message = "When calculating metrics, a range of allele 2 probability thresholds\ncan be used to calculate the metrics at each combination of allele 1 and allele 2 probability thresholds.\nThis sets the maximum allele 2 probability threshold.\nThe threshold increases in increments of 0.01.", position = "right")
     ), value=1, min = 0, max = 1)
   })
-  volumes = getVolumes()
+  #volumes = getVolumes()
   ## sample file
-  shinyFileChoose(input, "sample_GetFile", roots=volumes, session=session)
-  samplefile = reactive({parseFilePaths(volumes, input$sample_GetFile)})
+  roots = c(home="~", wd=".")
+  shinyFileChoose(input, "sample_GetFile", roots=roots, session=session, defaultRoot="home")
+  samplefile = reactive({parseFilePaths(roots, input$sample_GetFile)})
   observe({
     if(!is.null(samplefile)){
       output$sample_file = renderText({if(input$Submit==0){as.character(samplefile()$datapath)} else {return()}})
     }
   })
-
   ## freq files
-  shinyFileChoose(input, "freq_GetFile", roots=volumes, session=session)
-  freq = reactive({parseFilePaths(volumes, input$freq_GetFile)})
+  shinyFileChoose(input, "freq_GetFile", roots=roots, session=session, defaultRoot="home")
+  freq = reactive({parseFilePaths(roots, input$freq_GetFile)})
   if (!is.null(freq)) {
     observe({
       output$freq_file = renderText({if(input$Submit==0){as.character(freq()$datapath)} else {return()}})
     })
   }
 
-  shinyFileChoose(input, "freq_GetFile_major", roots=volumes, session=session)
-  freq_major = reactive({parseFilePaths(volumes, input$freq_GetFile_major)})
+  shinyFileChoose(input, "freq_GetFile_major", roots=roots, session=session, defaultRoot="home")
+  freq_major = reactive({parseFilePaths(roots, input$freq_GetFile_major)})
   if (!is.null(freq_major)) {
     observe({
       output$freq_file_major = renderText({if(input$Submit==0){as.character(freq_major()$datapath)} else {return()}})
     })
   }
 
-  shinyFileChoose(input, "freq_GetFile_minor", roots=volumes, session=session)
-  freq_minor = reactive({parseFilePaths(volumes, input$freq_GetFile_minor)})
+  shinyFileChoose(input, "freq_GetFile_minor", roots=roots, session=session, defaultRoot="home")
+  freq_minor = reactive({parseFilePaths(roots, input$freq_GetFile_minor)})
   if (!is.null(freq_minor)) {
     observe({
       output$freq_file_minor = renderText({if(input$Submit==0){as.character(freq_minor()$datapath)} else {return()}})
     })
   }
 
+  shinyFileChoose(input, "assay_GetFile", roots=roots, session=session, defaultRoot="home")
+  assay = reactive({parseFilePaths(roots, input$assay_GetFile)})
+  if (!is.null(assay)) {
+    observe({
+      output$assay_file = renderText({if(input$Submit==0){as.character(assay()$datapath)} else {return()}})
+    })
+  }
+
   ## refs
-  shinyDirChoose(input, "ref_GetFile", roots=volumes(), session=session)
-  refs = reactive({parseDirPath(volumes, input$ref_GetFile)})
+  shinyDirChoose(input, "ref_GetFile", roots=roots, session=session, defaultRoot="home")
+  refs = reactive({parseDirPath(roots, input$ref_GetFile)})
   if (!is.null(refs)) {
     observe({
       output$refs_file = renderText({if(input$Submit==0){as.character(refs())} else {return()}})
     })
   }
 
-  shinyDirChoose(input, "kin_prefix", roots=volumes(), session=session)
-  kin_inpath = reactive({parseDirPath(volumes, input$kin_prefix)})
+  shinyDirChoose(input, "kin_prefix", roots=roots, session=session, defaultRoot="home")
+  kin_inpath = reactive({parseDirPath(roots, input$kin_prefix)})
   observe({
     if(!is.null(kin_inpath)){
       output$kin_inpath = renderText({if(input$Submit==0){as.character(kin_inpath())} else {return()}})
@@ -426,37 +458,78 @@ server = function(input, output, session) {
 
 ## Input the sample manifest and run the workflow on each line (sample)
   observeEvent(input$Submit, {
-    sample_list = read.table(samplefile()$datapath, sep="\t", header=T)
-    date = glue("{Sys.Date()}_{format(Sys.time(), '%H_%M_%S')}")
-    create_config(date, input$twofreqs, ifelse(!isTruthy(freq()$datapath),  input$uploadfreq, freq()$datapath), ifelse(!isTruthy(freq_major()$datapath), input$uploadfreq_major, freq_major()$datapath), ifelse(!isTruthy(freq_minor()$datapath), input$uploadfreq_minor, freq_minor()$datapath), refs(), samplefile()$datapath, input$output, input$run_mixdeconv, input$uncond, input$ref_selector, input$method, input$sets, kin_inpath(), input$dynamicAT, input$staticAT, input$minimum_snps, input$A1_threshold, input$A2_threshold, input$A1_threshmin_metrics, input$A1_threshmax_metrics, input$A2_threshmin_metrics, input$A2_threshmax_metrics, input$major_selector, input$minor_selector, input$filter_missing, input$skip_ancestry, input$ancestry_snps, input$pcagroups)
-    if (isTruthy(refs())) {
-      withProgress(message = "Loading References", value = 0, {
-        if (!file.exists(glue("{refs()}/EFM_references.csv"))) {
-          refData = euroformix::sample_tableToList(data.frame(processing_ref_sample_reports(refs())))
-        } else {
-          refData = euroformix::sample_tableToList(euroformix::tableReader(glue("{refs()}/EFM_references.csv")))
+    if (!isTruthy(samplefile()$datapath)) {
+      showModal(modalDialog(
+        title = "Input Error",
+        "Please provide a sample manifest before proceeding.",
+        easyClose = TRUE,
+        footer = modalButton("Dismiss")
+      ))
+    } else if (input$method == "Calculate Metrics" | isTruthy(input$cond)) {
+      refData = NULL
+      if (isTruthy(refs())) {
+        withProgress(message = "Loading References", value = 0, {
+          if (file.exists(glue("{refs()}/EFM_references.rda"))) {
+            load(glue("{refs()}/EFM_references.rda"))
+          } else if (!file.exists(glue("{refs()}/EFM_references.csv"))) {
+            refData = convert_table_to_list(data.frame(processing_ref_sample_reports(refs())))
+          } else {
+            refs = data.frame(fread(glue("{refs()}/EFM_references.csv")))
+            refData = convert_table_to_list(refs)
+            save(refData, file=glue("{refs()}/EFM_references.rda"))
+          }
+        })
+      } else {
+        showModal(modalDialog(
+          title = "Input Error",
+          "Please provide reference genotypes before proceeding.",
+          easyClose = TRUE,
+          footer = modalButton("Dismiss")
+        ))
+      }
+    } else if (!isTruthy(input$cond) & !isTruthy(input$uncond)) {
+      showModal(modalDialog(
+        title = "Input Error",
+        "Please select either an Unconditioned or Conditioned analysis (or both!) before proceeding.",
+        easyClose = TRUE,
+        footer = modalButton("Dismiss")
+      ))
+    } else {
+      sample_list = suppressWarnings(euroformix::tableReader(samplefile()$datapath))
+      date = glue("{Sys.Date()}_{format(Sys.time(), '%H_%M_%S')}")
+      out_path = glue("{kin_inpath()}/snp_sets/{input$output}/")
+      create_config(date, input$twofreqs, ifelse(!isTruthy(freq()$datapath),  input$uploadfreq, freq()$datapath), ifelse(!isTruthy(freq_major()$datapath), input$uploadfreq_major, freq_major()$datapath), ifelse(!isTruthy(freq_minor()$datapath), input$uploadfreq_minor, freq_minor()$datapath), refs(), samplefile()$datapath, NULL, NULL, out_path, input$run_mixdeconv, input$uncond, input$ref_selector, input$method, input$sets, kin_inpath(), input$dynamicAT, input$staticAT, input$minimum_snps, input$A1_threshold, input$A2_threshold, input$A1_threshmin_metrics, input$A1_threshmax_metrics, input$A2_threshmin_metrics, input$A2_threshmax_metrics, input$major_selector, input$minor_selector, input$filter_missing, input$skip_ancestry, input$ancestry_snps, input$pcagroups)
+      withProgress(message = "Loading Allele Frequency Data", value = 0, {
+        freq_both = ifelse(!isTruthy(freq()$datapath), input$uploadfreq, freq()$datapath)
+        freq_major = ifelse(!isTruthy(freq_major()$datapath), input$uploadfreq_major, freq_major()$datapath)
+        freq_minor = ifelse(!isTruthy(freq_minor()$datapath), input$uploadfreq_minor, freq_minor()$datapath)
+        popFreq = load_freq(input$twofreqs, freq_both, freq_major, freq_minor)
+      })
+      if (isTruthy(assay()$datapath)) {
+        snp_positions = read.table(assay()$datapath, sep="\t", header=T)
+      } else {
+        snp_positions = mixder::kintelligence_snp_positions
+      }
+      withProgress(message = "Running Samples", value = 0, {
+        n = nrow(sample_list)
+        for (row in 1:n) {
+          id = sample_list[row, 1]
+          if (ncol(sample_list)>1) {
+            replicate_id = ifelse(is.na(sample_list[row, 2]), "", sample_list[row, 2])
+          } else {
+            replicate_id = ""
+          }
+          incProgress((row-1)/n, detail = glue("On Sample {row} of {n}"))
+            withCallingHandlers({
+              shinyjs::html(id = "text", html = "")
+              run_workflow(date, id, replicate_id, input$twofreqs, popFreq, refData, refs(), out_path, input$run_mixdeconv, input$uncond, input$ref_selector, input$method, input$sets, kin_inpath(), input$dynamicAT, input$staticAT, input$minimum_snps, input$A1_threshold, input$A2_threshold, input$A1_threshmin_metrics, input$A1_threshmax_metrics, input$A2_threshmin_metrics, input$A2_threshmax_metrics, input$major_selector, input$minor_selector, input$min_cont_prob, input$keep_bins, input$filter_missing, input$skip_ancestry, input$ancestry_snps, input$pcagroups, input$assay_choice, snp_positions)
+            },
+            message = function(m) {
+              shinyjs::html(id = "text", html = m$message, add = TRUE)
+            })
         }
       })
-    } else if(input$method == "Calculate Metrics" | isTruthy(input$cond)) {
-      stop("No references provided but selected conditioned analyses or calculating metrics. Please re-run!")
-    } else {
-      refData = NULL
     }
-    withProgress(message = "Running Samples", value = 0, {
-      n = nrow(sample_list)
-      for (row in 1:n) {
-        id = sample_list[row, 1]
-        replicate_id = ifelse(is.na(sample_list[row, 2]), "", sample_list[row, 2])
-        incProgress((row-1)/n, detail = glue("On Sample {row} of {n}"))
-          withCallingHandlers({
-            shinyjs::html(id = "text", html = "")
-            run_workflow(date, id, replicate_id, input$twofreqs, ifelse(!isTruthy(freq()$datapath), input$uploadfreq, freq()$datapath),ifelse(!isTruthy(freq_major()$datapath), input$uploadfreq_major, freq_major()$datapath), ifelse(!isTruthy(freq_minor()$datapath), input$uploadfreq_minor, freq_minor()$datapath), refData, refs(), input$output, input$run_mixdeconv, input$uncond, input$ref_selector, input$method, input$sets, kin_inpath(), input$dynamicAT, input$staticAT, input$minimum_snps, input$A1_threshold, input$A2_threshold, input$A1_threshmin_metrics, input$A1_threshmax_metrics, input$A2_threshmin_metrics, input$A2_threshmax_metrics, input$major_selector, input$minor_selector, input$min_cont_prob, input$keep_bins, input$filter_missing, input$skip_ancestry, input$ancestry_snps, input$pcagroups)
-          },
-          message = function(m) {
-            shinyjs::html(id = "text", html = m$message, add = TRUE)
-          })
-      }
-    })
   })
 }
 

R/assigned_A2.R

@@ -12,15 +12,21 @@
 #'
 #' @param x Data frame of allele calls
 #' @param thresh Allele 2 probability threshold
+#' @param pos SNP positions data frame
 #'
 #' @return data frame
 #' @export
-assigned_A2 = function(x, thresh) {
+assigned_A2 = function(x, thresh, pos) {
   x$Allele2 = ifelse(x$A2_Prob>=thresh & x$A2!="99", x$A2, x$A1)
   names(x)[names(x) == "A1"] = "Allele1"
   x$Locus = ifelse(x$Locus=="RS201326893_Y152OCH", "rs201326893_Y152OCH", ifelse(x$Locus=="N29INSA", "N29insA", tolower(x$Locus)))
-  rsids = merge(mixder::kintelligence_snp_positions, x, by.x="rsid", by.y="Locus")
-  final_x = rsids[,c("rsid", "chromosome", "position", "Allele1", "Allele2")]
-  final_x_sorted = final_x[order(final_x$chromosome, final_x$position),]
+  colname = colnames(pos)
+  lcol = grep("rsid|marker|locus|snp", colname, ignore.case=TRUE, value=TRUE)
+  chrcol = grep("chr", colname, ignore.case=TRUE, value=TRUE)
+  poscol = grep("pos", colname, ignore.case=TRUE, value=TRUE)
+  rsids = merge(pos, x, by.x=lcol, by.y="Locus")
+  final_x = rsids[,c(lcol, chrcol, poscol, "Allele1", "Allele2")]
+  final_x_sorted = final_x %>%
+    arrange(as.numeric(!!sym(chrcol)), as.numeric(!!sym(poscol)))
   return(final_x_sorted)
 }

R/check_nas.R

@@ -16,13 +16,16 @@
 #' @export
 #'
 check_nas = function(df) {
-  if (sum(!is.na(df$Allele.4))==0) {
-    df[,c("Allele.4","Height.4")] = list(NULL)
+  if ("Allele.4" %in% names(df)) {
+    if (sum(!is.na(df$Allele.4))==0) {
+      df[,c("Allele.4","Height.4")] = list(NULL)
+    }
   }
-  if (sum(!is.na(df$Allele.3))==0) {
-    df[,c("Allele.3","Height.3")] = list(NULL)
+  if ("Allele.3" %in% names(df)) {
+    if (sum(!is.na(df$Allele.3))==0) {
+      df[,c("Allele.3","Height.3")] = list(NULL)
+    }
+    df[is.na(df)] = 0
   }
-  df[is.na(df)] = 0
-  #df %>% replace(is.na(.data), 0)
   return(df)
 }