diff --git a/DESCRIPTION b/DESCRIPTION
index fc88253..872198d 100644
--- a/DESCRIPTION
+++ b/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
Maintainer: Rebecca Mitchell
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.
diff --git a/NAMESPACE b/NAMESPACE
index 3bc58b8..28ae4e2 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index fe63512..4608b40 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,17 @@
+# 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
+- Improved error handling in Shiny app (pop-up window will appear if required files or settings are missing)
+- Added option to change number of threads used by EFM; useful if running on a server or grid (default is 0)
+
+
# mixder 0.9.0
## Changes in version 0.9.0
diff --git a/R/ancestry_prediction.R b/R/ancestry_prediction.R
index 918c284..bd567ae 100644
--- a/R/ancestry_prediction.R
+++ b/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)
diff --git a/R/app.R b/R/app.R
index 8e723d0..58cbd61 100644
--- a/R/app.R
+++ b/R/app.R
@@ -62,11 +62,12 @@ mixder = function() {
name = "question-circle",
)
) |>
- add_prompt(message = "Select a folder containing the Mixture Sample Reports.", position = "right"),
+ add_prompt(message = "Select a folder containing the Mixture Sample Reports or a TSV file containing mixture genotypes.", position = "right"),
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,92 +395,154 @@ 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()
## sample file
- shinyFileChoose(input, "sample_GetFile", roots=volumes, session=session)
- samplefile = reactive({parseFilePaths(volumes, input$sample_GetFile)})
- observe({
- if(!is.null(samplefile)){
- output$sample_file = renderText({if(input$Submit==0){as.character(samplefile()$datapath)} else {return()}})
- }
- })
-
+ roots = c(home="~", wd=".")
+ shinyFileChoose(input, "sample_GetFile", roots=roots, session=session, defaultRoot="home")
+ samplefile = reactive({parseFilePaths(roots, input$sample_GetFile)})
+ if(!is.null(samplefile)){
+ observe({
+ output$sample_file = renderText({as.character(samplefile()$datapath)})
+ })
+ }
## 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()}})
+ output$freq_file = renderText({as.character(freq()$datapath)})
})
}
- 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()}})
+ output$freq_file_major = renderText({as.character(freq_major()$datapath)})
})
}
- 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()}})
+ output$freq_file_minor = renderText({as.character(freq_minor()$datapath)})
+ })
+ }
+
+ 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({as.character(assay()$datapath)})
})
}
## 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()}})
+ output$refs_file = renderText({as.character(refs())})
})
}
- 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()}})
+ output$kin_inpath = renderText({as.character(kin_inpath())})
}
})
## 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 = "Missing Input",
+ "Please provide a sample manifest before proceeding.",
+ easyClose = TRUE,
+ footer = modalButton("Dismiss")
+ ))
+ } else if (!isTruthy(kin_inpath())) {
+ showModal(modalDialog(
+ title = "Missing Input",
+ "Please provide a folder containing mixture data (Kintelligence Sample Reports or a .tsv file of genotypes).",
+ 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 if (input$assay_choice=="Custom" & !isTruthy(assay()$datapath)) {
+ showModal(modalDialog(
+ title = "Missing Input",
+ "When running a custom SNP panel, please provide the SNP positions file.",
+ easyClose = TRUE,
+ footer = modalButton("Dismiss")
+ ))
+ } else if ((input$method == "Calculate Metrics" | isTruthy(input$cond)) & !isTruthy(refs())) {
+ showModal(modalDialog(
+ title = "Missing Input",
+ "Please provide reference genotypes before proceeding.",
+ easyClose = TRUE,
+ footer = modalButton("Dismiss")
+ ))
+ } else {
+ 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 {
+ refData=NULL
+ }
+ 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, input$assay_choice, assay()$datapath)
+ 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)
- })
- }
- })
})
}
diff --git a/R/assigned_A2.R b/R/assigned_A2.R
index bd3b720..00e0a34 100644
--- a/R/assigned_A2.R
+++ b/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)
}
diff --git a/R/check_nas.R b/R/check_nas.R
index 40a796d..0457169 100644
--- a/R/check_nas.R
+++ b/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)
}
diff --git a/R/checking_af.R b/R/checking_af.R
index 7c72ec4..a306723 100644
--- a/R/checking_af.R
+++ b/R/checking_af.R
@@ -1,19 +1,36 @@
-
+# -------------------------------------------------------------------------------------------------
+# Copyright (c) 2024, DHS.
+#
+# This file is part of MixDeR and is licensed under the BSD license: see LICENSE.
+#
+# This software was prepared for the Department of Homeland Security (DHS) by the Battelle National
+# Biodefense Institute, LLC (BNBI) as part of contract HSHQDC-15-C-00064 to manage and operate the
+# National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and
+# Development Center.
+# -------------------------------------------------------------------------------------------------
#' Checking AF file for correct format; formatting if necessary
#'
#' @param affile Allele Frequency file
-#' @param outpath Path to write formatted AF file to if necessary
+#' @param contrib contributor assigned to AF file
#'
#' @return Frequency file loaded for EFM
#' @export
-checking_af = function(affile, outpath) {
- af=utils::read.csv(affile, header=T)
- if (isTruthy(af[c(1:4),1] == c("A", "C", "G", "T"))) {
- finalfreq = euroformix::freqImport(affile)[[1]]
+#'
+#' @importFrom data.table fread
+checking_af = function(affile) {
+ if (grepl(".rda", affile, ignore.case = TRUE)) {
+ load(affile)
} else {
- af_formatted = format_af(af)
- utils::write.csv(af_formatted, glue("{outpath}/custom_AF_formatted.csv"), row.names=F, quote=F)
- finalfreq = euroformix::freqImport(glue("{outpath}/custom_AF_formatted.csv"))[[1]]
+ af=data.frame(fread(affile, header=T, sep=","))
+ filename = gsub(".csv","", affile)
+ if (!isTruthy(af[c(1:4),1] == c("A", "C", "G", "T"))) {
+ af_formatted = format_af(af)
+ utils::write.csv(af_formatted, glue("{filename}_EFMformatted.csv"), row.names=F, quote=F)
+ } else {
+ af_formatted = af
+ }
+ popFreq=read_in_freq(af_formatted)
+ save(popFreq, file=glue("{filename}.rda"))
}
- return(finalfreq)
+ return(popFreq)
}
diff --git a/R/convert_table_to_list.R b/R/convert_table_to_list.R
new file mode 100644
index 0000000..077546b
--- /dev/null
+++ b/R/convert_table_to_list.R
@@ -0,0 +1,64 @@
+# -------------------------------------------------------------------------------------------------
+# Copyright (c) 2026, DHS.
+#
+# This file is part of MixDeR and is licensed under the BSD license: see LICENSE.
+#
+# This software was prepared for the Department of Homeland Security (DHS) by the Battelle National
+# Biodefense Institute, LLC (BNBI) as part of contract HSHQDC-15-C-00064 to manage and operate the
+# National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and
+# Development Center.
+# -------------------------------------------------------------------------------------------------
+#' Convert table to list
+#'
+#' @param input table
+#'
+#' @returns list of table
+#' @export
+#'
+#' @importFrom data.table fread
+#' @import rrapply
+convert_table_to_list = function(table) {
+ #table=data.frame(fread(input, header=T))
+ outL=rrapply(table[,c(1:4)], how="unmelt")
+ nestedlist = lapply(split(table, table$Sample.Name, drop = TRUE),
+ function(x) c(split(x, x[["Marker"]], drop = TRUE)))
+ colname = colnames(table) #colnames in file
+ lind = grep("marker",tolower(colname),fixed=TRUE) #locus col-ind
+ A_ind = grep("allele",tolower(colname),fixed=TRUE) #allele col-ind
+ H_ind = grep("height",tolower(colname),fixed=TRUE) #height col-ind
+ if(length(lind)==0) lind = grep("loc",tolower(colname),fixed=TRUE) #try another name
+ sind = grep("sample",tolower(colname),fixed=TRUE) #sample col-ind
+ if(length(sind)>1) sind = sind[grep("name",tolower(colname[sind]),fixed=TRUE)] #use only sample name
+ samplenames = unique(as.character(table[,sind])) #sample names
+ outL = list() #Init outList (insert non-empty characters):
+ for(samplename in samplenames) { #for each sample in matrix
+ filt_table=subset(table, table[,sind]==samplename)
+ locs = unique(toupper(filt_table[,lind])) #locus names: Use uniques and Convert to upper case
+ for(loc in locs) { #for each locus
+ loc_lower=ifelse(loc=="N29INSA", "N29insA", ifelse(loc=="RS201326893_Y152OCH", "rs201326893_Y152OCH", tolower(loc)))
+ if (length(A_ind) == 2) {
+ outL[[samplename]][[loc]]$adata = c(nestedlist[[samplename]][[loc_lower]][[A_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[2]]]])
+ if(length(H_ind)>0) {
+ outL[[samplename]][[loc]]$hdata = c(nestedlist[[samplename]][[loc_lower]][[H_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[2]]]])
+ }
+ } else if (length(A_ind) == 3) {
+ if (!is.na(nestedlist[[samplename]][[loc_lower]][[A_ind[[3]]]])) {
+ outL[[samplename]][[loc]]$adata = c(nestedlist[[samplename]][[loc_lower]][[A_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[3]]]])
+ outL[[samplename]][[loc]]$hdata = c(nestedlist[[samplename]][[loc_lower]][[H_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[3]]]])
+ } else {
+ outL[[samplename]][[loc]]$adata = c(nestedlist[[samplename]][[loc_lower]][[A_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[2]]]])
+ outL[[samplename]][[loc]]$hdata = c(nestedlist[[samplename]][[loc_lower]][[H_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[2]]]])
+ }
+ } else if (length(A_ind) == 4) {
+ if (!is.na(nestedlist[[samplename]][[loc_lower]][[A_ind[[4]]]])) {
+ outL[[samplename]][[loc]]$adata = c(nestedlist[[samplename]][[loc_lower]][[A_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[3]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[4]]]])
+ outL[[samplename]][[loc]]$hdata = c(nestedlist[[samplename]][[loc_lower]][[H_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[3]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[4]]]])
+ } else {
+ outL[[samplename]][[loc]]$adata = c(nestedlist[[samplename]][[loc_lower]][[A_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[A_ind[[3]]]])
+ outL[[samplename]][[loc]]$hdata = c(nestedlist[[samplename]][[loc_lower]][[H_ind[[1]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[2]]]], nestedlist[[samplename]][[loc_lower]][[H_ind[[3]]]])
+ }
+ }
+ }
+ }
+ return(outL)
+}
diff --git a/R/create_config.R b/R/create_config.R
index c2ae1a3..c68ad19 100644
--- a/R/create_config.R
+++ b/R/create_config.R
@@ -42,13 +42,19 @@
#' @param skipancestry TRUE/FALSE whether to skip ancestry prediction step
#' @param pcasnps SNPs used for PCA (ancestry prediction)
#' @param pcagroups Groups used for PCA (ancestry prediction), either Superpopulations or Subpopulations
+#' @param assay name of sequencing assay, either `kintelligence` or `custom`
+#' @param pos Path to file containing SNP positions of custom assay, required if assay==`custom`
#'
#' @export
#'
-create_config = function(date, twofreqs, freq_all, freq_major, freq_minor, refs, sample_manifest, sample, replicate, out_path, run_mixdeconv, unconditioned, cond, method, sets, kinpath, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, major, minor, filter_missing, skipancestry, pcasnps, pcagroups){
+create_config = function(date, twofreqs, freq_all, freq_major, freq_minor, refs, sample_manifest, sample, replicate, out_path, run_mixdeconv, unconditioned, cond, method, sets, kinpath, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, major, minor, filter_missing, skipancestry, pcasnps, pcagroups, assay, pos){
config = setNames(data.frame(matrix(ncol=2, nrow=0)), c("Setting", "Value"))
config = rbind(config, data.frame(Setting="MixDeR Version:", Value=getNamespaceVersion("mixder")[["version"]]))
config = rbind(config, data.frame(Setting="EuroForMix Version:", Value=getNamespaceVersion("euroformix")[["version"]]))
+ onfig = rbind(config, data.frame(Setting="Assay:", Value=assay))
+ if (assay == "custom") {
+ config = rbind(config, data.frame(Setting="SNP hg19 positions file for custom assay:", Value=pos))
+ }
if (isTruthy(sample_manifest)) {
config = rbind(config, data.frame(Setting="Path to sample manifest:", Value=sample_manifest))
} else if (isTruthy(sample)) {
@@ -73,7 +79,7 @@ create_config = function(date, twofreqs, freq_all, freq_major, freq_minor, refs,
config = rbind(config, data.frame(Setting="Frequency data for Minor Contributor:", Value=freq_all))
}
}
- config = rbind(config, data.frame(Setting="Output path:", Value=glue("{kinpath}/snp_sets/{out_path}/")))
+ config = rbind(config, data.frame(Setting="Output path:", Value=out_path))
config = rbind(config, data.frame(Setting="Number of SNP sets:", Value=sets))
config = rbind(config, data.frame(Setting="Minimum number of SNPs:", Value=minimum_snps))
config = rbind(config, data.frame(Setting="Static AT:", Value=staticAT))
@@ -102,6 +108,6 @@ create_config = function(date, twofreqs, freq_all, freq_major, freq_minor, refs,
config = rbind(config, data.frame(Setting="Assumed major contributor:", Value=major))
config = rbind(config, data.frame(Setting="Assumed minor contributor:", Value=minor))
}
- dir.create(file.path(kinpath, "snp_sets", out_path, "config_log_files", date), showWarnings = FALSE, recursive=TRUE)
- write.table(config, glue("{kinpath}/snp_sets/{out_path}/config_log_files/{date}/config_settings_run_{date}.txt"), row.names=F, quote=F, col.names=T, sep="\t")
+ dir.create(file.path(out_path, "config_log_files", date), showWarnings = FALSE, recursive=TRUE)
+ write.table(config, glue("{out_path}/config_log_files/{date}/config_settings_run_{date}.txt"), row.names=F, quote=F, col.names=T, sep="\t")
}
diff --git a/R/create_evid.R b/R/create_evid.R
index fa905fd..ae2aa9a 100644
--- a/R/create_evid.R
+++ b/R/create_evid.R
@@ -22,7 +22,8 @@ create_evid = function(sample, rep_sample, inpath) {
if (check_reads(evidfn)) {
evidData = vector()
} else {
- evidData = euroformix::sample_tableToList(read_in_table(evidfn))
+ evid_df_form = read_in_table(evidfn)
+ evidData = convert_table_to_list(evid_df_form)
}
} else {
if (check_reads(evidfn)) {
@@ -45,7 +46,7 @@ create_evid = function(sample, rep_sample, inpath) {
final_dfs = equalize_samples(corrected_tables[[1]], corrected_tables[[2]])
}
merge_table = rbind(final_dfs[[1]], final_dfs[[2]])
- evidData = euroformix::sample_tableToList(merge_table)
+ evidData = convert_table_to_list(merge_table)
}
return(evidData)
}
diff --git a/R/create_evid_all.R b/R/create_evid_all.R
index 831de98..7697ee1 100644
--- a/R/create_evid_all.R
+++ b/R/create_evid_all.R
@@ -16,11 +16,15 @@
#' @param keep_bins To use existing SNP bins or create new bins (and files)
#' @export
-create_evid_all = function(inpath, id, nsets, keep_bins){
+create_evid_all = function(inpath, id, nsets, keep_bins, assay){
evid_all = processing_evid_sample_reports(inpath, id)
evid_all$Total_Reads = rowSums(evid_all[,grepl("Height",colnames(evid_all))], na.rm = TRUE)
evid_sort = evid_all[order(evid_all$Total_Reads),]
- bin_size = round(10039/nsets)
+ if (assay == "Kintelligence") {
+ bin_size = round(10039/nsets)
+ } else {
+ bin_size = round(nrow(evid_sort)/nsets)
+ }
dir.create(file.path(inpath, "snp_sets"), showWarnings = FALSE, recursive=TRUE)
evid_final = cbind(Sample.Name = id, evid_sort)
write.table(evid_final, glue("{inpath}/snp_sets/{id}_snpsetscombined_evidence.tsv"), row.names=F, quote=F, col.names=T, sep="\t")
diff --git a/R/create_gedmatchpro_report.R b/R/create_gedmatchpro_report.R
index faaff87..82f7eef 100644
--- a/R/create_gedmatchpro_report.R
+++ b/R/create_gedmatchpro_report.R
@@ -23,12 +23,13 @@
#' @param A2max Maximum value for allele 2 probability threshold
#' @param minor_threshold If apply allele 1 probability to minor contributor
#' @param filter_missing TRUE/FALSE to filter SNPs with missing allele 2 values
+#' @param pos SNP position data frame
#'
#' @return list of two data frames, the GEDmatch PRO report and a data frame of metrics calculated for that report
#' @export
-create_gedmatchpro_report = function(path, x, contrib, contrib_status, minimum_snps, A1, A2, A1min, A1max, A2min, A2max, minor_threshold, filter_missing) {
+create_gedmatchpro_report = function(path, x, contrib, contrib_status, minimum_snps, A1, A2, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, pos) {
formatted_df = suppressWarnings(process_efm_files(x, contrib, NULL, minimum_snps, A1min, A1max, A2min, A2max, metrics=FALSE, filter_missing))
gedmatch_metrics = gedmatch_metrics(formatted_df, A1, A2, minimum_snps, path)
- report = filter_alleles(formatted_df, contrib_status, minimum_snps, A1, A2, minor_threshold, filter_missing)
+ report = filter_alleles(formatted_df, contrib_status, minimum_snps, A1, A2, minor_threshold, filter_missing, pos)
return(list(report, gedmatch_metrics[[1]], gedmatch_metrics[[2]]))
}
diff --git a/R/filter_alleles.R b/R/filter_alleles.R
index ba1a55c..abbf742 100644
--- a/R/filter_alleles.R
+++ b/R/filter_alleles.R
@@ -18,10 +18,11 @@
#' @param A2_threshold Allele 2 probability threshold
#' @param minor_threshold If apply allele 1 probability threshold to minor contributor
#' @param filter_missing TRUE/FALSE to filter SNPs with missing allele 2 values
+#' @param pos SNP position data frame
#'
#' @return Data frame with filtered allele calls
#' @export
-filter_alleles = function(all_files, contrib_status, minimum_snps, A1_threshold, A2_threshold, minor_threshold, filter_missing) {
+filter_alleles = function(all_files, contrib_status, minimum_snps, A1_threshold, A2_threshold, minor_threshold, filter_missing, pos) {
. = NULL
if (filter_missing) {
all_files = subset(all_files, !(all_files$A2 == 99 & all_files$A2_Prob%
.[c(1:minimum_snps),]
}
- final_df = assigned_A2(filt_df, A2_threshold)
+ final_df = assigned_A2(filt_df, A2_threshold, pos)
return(final_df)
}
diff --git a/R/format_ref.R b/R/format_ref.R
index ef9c34b..1042703 100644
--- a/R/format_ref.R
+++ b/R/format_ref.R
@@ -10,22 +10,17 @@
# -------------------------------------------------------------------------------------------------
#' @title Formatting reference file
#'
-#' @param refData Data frame of the reference genotypes
-#' @param refid Reference order in the reference genotypes file
+#' @param refid ID of the reference genotypes of interest
#' @param refs Path of the reference genotypes file
#'
#' @return Data frame containing reference of interest genotypes
#' @export
#'
#' @importFrom rlang .data
-format_ref = function(refData, refid, refs) {
- if (is.numeric(refid)) {
- ref_int = names(refData)[refid]
- } else {
- ref_int = refid
- }
- ref_profile = euroformix::tableReader(glue("{refs}/EFM_references.csv")) %>%
- filter(.data$Sample.Name == ref_int)
+#' @importFrom data.table fread
+format_ref = function(refid, refs) {
+ ref_profile = data.frame(fread(glue("{refs}/EFM_references.csv"))) %>%
+ filter(.data$Sample.Name == refid)
ref_profile$A1_order = ifelse(ref_profile$`Allele1`>ref_profile$`Allele2`, ref_profile$`Allele1`, ref_profile$`Allele2`)
ref_profile$A2_order = ifelse(ref_profile$`Allele1`>ref_profile$`Allele2`, ref_profile$`Allele2`, ref_profile$`Allele1`)
return(ref_profile)
diff --git a/R/get_ids.R b/R/get_ids.R
index d6a0de3..88d0d12 100644
--- a/R/get_ids.R
+++ b/R/get_ids.R
@@ -15,7 +15,10 @@
#' @return command to get list of files
#' @export
get_ids = function(inpath) {
- if (file.exists(glue("{inpath}/EFM_references.csv"))) {
+ if (file.exists(glue("{inpath}/EFM_references.rda"))) {
+ load(glue("{inpath}/EFM_references.rda"))
+ return(names(refData))
+ } else if (file.exists(glue("{inpath}/EFM_references.csv"))) {
all_samples = utils::read.csv(glue("{inpath}/EFM_references.csv"))
if ("Sample.Name" %in% colnames(all_samples)) {
return(unique(all_samples$Sample.Name))
diff --git a/R/load_freq.R b/R/load_freq.R
index ed5be2d..a971c62 100644
--- a/R/load_freq.R
+++ b/R/load_freq.R
@@ -8,9 +8,8 @@
# National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and
# Development Center.
# -------------------------------------------------------------------------------------------------
-#' Load in frequency data
+#' Load frequency data
#'
-#'@param out_path outpath directory
#' @param twofreqs If two different allele frequency files are to be used
#' @param freq_both_input Allele frequency file, if only using one
#' @param freq_major_input Allele frequency file for major contributor
@@ -18,11 +17,11 @@
#'
#' @return list of major AF data data frame and minor AF data data frame
#' @export
-load_freq = function(out_path, twofreqs, freq_both_input, freq_major_input, freq_minor_input) {
+load_freq = function(twofreqs, freq_both_input, freq_major_input, freq_minor_input) {
freqs_allowed = c("Global - 1000G", "Global - gnomAD", "AFR - 1000G", "AMR - 1000G", "EAS - 1000G", "EUR - 1000G", "SAS - 1000G")
if (!twofreqs) {
if (file.exists(freq_both_input)) {
- freq_minor = checking_af(freq_both_input, out_path)
+ freq_minor = checking_af(freq_both_input)
freq_major = freq_minor
} else {
freq_both = decode_freq_name(tolower(freq_both_input))
@@ -56,7 +55,7 @@ load_freq = function(out_path, twofreqs, freq_both_input, freq_major_input, freq
}
} else {
if (file.exists(freq_major_input)) {
- freq_major = checking_af(freq_major_input, out_path)
+ freq_major = checking_af(freq_major_input)
} else {
freq_major = decode_freq_name(tolower(freq_major_input))
if (!freq_major %in% freqs_allowed) {
@@ -81,7 +80,7 @@ load_freq = function(out_path, twofreqs, freq_both_input, freq_major_input, freq
}
}
if (file.exists(freq_minor_input)) {
- freq_minor = checking_af(freq_minor_input, out_path)
+ freq_minor = checking_af(freq_minor_input)
} else {
freq_minor = decode_freq_name(tolower(freq_minor_input))
if (!freq_minor %in% freqs_allowed) {
diff --git a/R/read_in_freq.R b/R/read_in_freq.R
new file mode 100644
index 0000000..2e133b9
--- /dev/null
+++ b/R/read_in_freq.R
@@ -0,0 +1,34 @@
+# -------------------------------------------------------------------------------------------------
+# Copyright (c) 2026, DHS.
+#
+# This file is part of MixDeR and is licensed under the BSD license: see LICENSE.
+#
+# This software was prepared for the Department of Homeland Security (DHS) by the Battelle National
+# Biodefense Institute, LLC (BNBI) as part of contract HSHQDC-15-C-00064 to manage and operate the
+# National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and
+# Development Center.
+# -------------------------------------------------------------------------------------------------
+#' Read in Frequency tables already formatted for EFM
+#'
+#' @param freq PATH to frequency table
+#'
+#' @returns list of frequency data in correct EFM list format
+#' @export
+#'
+#' @importFrom data.table fread
+#'
+read_in_freq = function(freq){
+ #efmaf=fread(freq, header=T, sep=",")
+ tab=data.frame(freq)
+ Anames = tab[,1] #first column is allele frequencies
+ tab = tab[,-1,drop=FALSE]
+ freqlist = vector("list", ncol(tab))
+ for(j in 1:ncol(tab)) { #for each locus
+ tmp = tab[,j]
+ tmp2 = tmp[!is.na(tmp) & as.numeric(tmp)>0] #require that allele is not NA and is>0
+ names(tmp2) = Anames[!is.na(tmp)]
+ freqlist[[j]] = tmp2
+ }
+ names(freqlist) = toupper(colnames(tab)) #LOCUS-names are assigned as Upper-case! This is important to do!
+ return(freqlist)
+}
diff --git a/R/read_in_table.R b/R/read_in_table.R
index aa4bd2e..56fab49 100644
--- a/R/read_in_table.R
+++ b/R/read_in_table.R
@@ -15,14 +15,15 @@
#' @return data frame of SNPs
#' @export
#'
+#' @importFrom data.table fread
read_in_table = function(df) {
- evid_df = read.table(df, header=T, sep="\t")
+ evid_df = data.frame(fread(df, header=T, sep="\t"))
myColClasses = sapply(evid_df, class)
needed_cols = c("Allele.1", "Allele.2", "Allele.3", "Allele.4")
myColClasses =
ifelse(names(myColClasses) %in% needed_cols,
"character",
myColClasses)
- evid_df_final = read.table(df, header=T, sep="\t", colClasses=myColClasses)
- return(evid_df_final)
+ evid_df_form = data.frame(fread(df, header=T, sep="\t", colClasses=myColClasses))
+ return(evid_df_form)
}
diff --git a/R/run_efm.R b/R/run_efm.R
index 7438ebc..a87856e 100644
--- a/R/run_efm.R
+++ b/R/run_efm.R
@@ -19,28 +19,30 @@
#' @param out_path Output path
#' @param attable AT table
#' @param nsets Number of SNP sets
-#' @param ancestry if not skipping ancestry prediction
+#' @param ancestry TRUE/FALSE if to skip ancestry
+#' @param assay assay used (kintelligence or custom)
#' @param cond Sample IDs to condition on
#' @param uncond TRUE/FALSE if performing unconditioned analysis
#' @param keep_bins To use existing SNP bins or create new bins (and files)
+#' @param threads number of threads used for EFM
#'
#' @export
#'
#' @import parallel
#'
-run_efm = function(date, popFreq, refData, id, replicate_id, inpath, out_path, attable, nsets, ancestry, cond = NULL, uncond=TRUE, keep_bins=TRUE) {
+run_efm = function(date, popFreq, refData, id, replicate_id, inpath, out_path, attable, nsets, ancestry, assay, cond = NULL, uncond=TRUE, keep_bins=TRUE, threads=0, parallel=FALSE) {
if (!ancestry) {
new_path = glue("{out_path}/ancestry_prediction/")
} else {
new_path = out_path
}
if (replicate_id == "") {
- create_evid_all(inpath, id, nsets, keep_bins)
+ create_evid_all(inpath, id, nsets, keep_bins, assay)
write_path = glue("{new_path}Single/{id}")
log_name = id
} else {
- create_evid_all(inpath, id, nsets, keep_bins)
- create_evid_all(inpath, replicate_id, nsets, keep_bins)
+ create_evid_all(inpath, id, nsets, keep_bins, assay)
+ create_evid_all(inpath, replicate_id, nsets, keep_bins, assay)
write_path = glue("{new_path}Replicates/{id}")
log_name = glue("{id}_{replicate_id}")
}
@@ -51,9 +53,17 @@ run_efm = function(date, popFreq, refData, id, replicate_id, inpath, out_path, a
} else {
ids = NULL
}
- results = list()
- for (i in 1:nsets) {
- results[[i]] = run_indiv_efm_set(i, ids, snps_input, popFreq, refData, id, replicate_id, write_path, attable, cond=cond, uncond=uncond)
+ if (parallel) {
+ numCores = ifelse(detectCores()"))
+ sample_at = create_at(evidData, sample, replicate, attable)
+ }
ratio_row = list()
ratio_row[glue("Set{i}_C1_Prob_uncond")] = NA
ratio_row[glue("Set{i}_C2_Prob_uncond")] = NA
if (uncond) {
+ print(glue("Running unconditioned analysis for set {i}"))
message("Running unconditioned mixture deconvolution ")
dir.create(file.path(write_path, "unconditioned"), showWarnings = FALSE, recursive=TRUE)
##unconditioned analysis
@@ -50,9 +63,9 @@ run_indiv_efm_set = function(i, ids, snps_input, popFreq, refData, id, replicate
repeat {
message(glue("Running unconditioned analysis for set {i}, attempt #{repeat_num+1} "))
if (substr(efm_v, 1,3)!="4.0" & substr(efm_v, 1,2) != "3.") {
- uncond_results = euroformix::calcMLE(2, evidData, popFreq, AT=sample_at, BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, resttol=0)
+ uncond_results = euroformix::calcMLE(2, evidData, popFreq, AT=sample_at, BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, resttol=0, maxThreads=threads)
} else {
- uncond_results = euroformix::calcMLE(2, evidData, popFreq, AT=sample_at, BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01)
+ uncond_results = euroformix::calcMLE(2, evidData, popFreq, AT=sample_at, BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, maxThreads=threads)
}
uncond_finaltable = euroformix::deconvolve(uncond_results)
if (check_allele_probabilities(data.frame(uncond_finaltable[["table4"]]), i)) break
@@ -63,14 +76,10 @@ run_indiv_efm_set = function(i, ids, snps_input, popFreq, refData, id, replicate
if (repeat_num < 10) {
ratio_row[glue("Set{i}_C1_Prob_uncond")] = uncond_results[["fit"]][["thetahat2"]][["Mix-prop. C1"]]
ratio_row[glue("Set{i}_C2_Prob_uncond")] = uncond_results[["fit"]][["thetahat2"]][["Mix-prop. C2"]]
- #ratio_row = list(Set_uncond=i, C1_Prob_uncond=uncond_results[["fit"]][["thetahat2"]][["Mix-prop. C1"]], C2_Prob_uncond=uncond_results[["fit"]][["thetahat2"]][["Mix-prop. C2"]])
- #uncond_ratios = rbind(uncond_ratios, ratio_row)
write.table(uncond_finaltable[["table4"]], glue("{write_path}/unconditioned/{id}_set{i}_uncond.tsv"), quote=F, row.names=F, sep="\t")
- #uncond_finaltable_all = rbind(uncond_finaltable_all, uncond_finaltable[["table4"]])
- #set_results = uncond_finaltable[["table4"]]
+ write.table(uncond_finaltable[["table3"]], glue("{write_path}/unconditioned/{id}_set{i}_uncond_table3.tsv"), quote=F, row.names=F, sep="\t")
} else {
message(glue("Repeated unconditioned analysis 10 times unsuccessfully. Will skip set {i}! "))
- #set_results = data.frame()
}
}
final_list = c(ratio_row)
@@ -80,20 +89,21 @@ run_indiv_efm_set = function(i, ids, snps_input, popFreq, refData, id, replicate
cond_vector = rep.int(0, total_refs)
for (cond_on in cond) {
ratio_row = list()
+ print(glue("Running conditioned analysis on {cond_on} for set {i}"))
message(glue("Running mixture deconvolution conditioned on {cond_on} "))
- #nam = glue("df_{cond_on}")
list_num = match(cond_on, names(refData))
dir.create(file.path(write_path, glue("/conditioned/cond_on_{cond_on}")), showWarnings = FALSE, recursive=TRUE)
message(glue("Running conditioned analysis on {cond_on} for set {i}. "))
if (substr(efm_v, 1,3)!="4.0" & substr(efm_v, 1,2) != "3.") {
- condresults = euroformix::calcMLE(2, evidData, popFreq, refData, AT=sample_at, condOrder=replace(cond_vector, list_num, 1), BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, resttol=0)
+ condresults = euroformix::calcMLE(2, evidData, popFreq, refData, AT=sample_at, condOrder=replace(cond_vector, list_num, 1), BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, resttol=0, maxThreads=threads)
} else {
- condresults = euroformix::calcMLE(2, evidData, popFreq, refData, AT=sample_at, condOrder=replace(cond_vector, list_num, 1), BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01)
+ condresults = euroformix::calcMLE(2, evidData, popFreq, refData, AT=sample_at, condOrder=replace(cond_vector, list_num, 1), BWS=FALSE, FWS=FALSE, DEG=FALSE, steptol=0.001, pC=0.01, lambda=0.05, fst=0.01, maxThreads=threads)
}
final_condresults = euroformix::deconvolve(condresults)
ratio_row[glue("Set{i}_C1_Prob_cond_on_{cond_on}")] = condresults[["fit"]][["thetahat2"]][["Mix-prop. C1"]]
ratio_row[glue("Set{i}_C2_Prob_cond_on_{cond_on}")] = condresults[["fit"]][["thetahat2"]][["Mix-prop. C2"]]
write.table(final_condresults[["table4"]], glue("{write_path}/conditioned/cond_on_{cond_on}/unknown_cond_on_{cond_on}_set{i}.tsv"), quote=F, row.names=F, sep="\t")
+ write.table(final_condresults[["table3"]], glue("{write_path}/conditioned/cond_on_{cond_on}/unknown_cond_on_{cond_on}_set{i}_table3.tsv"), quote=F, row.names=F, sep="\t")
final_list = c(final_list, ratio_row)
}
}
diff --git a/R/run_mixder_ancestry.R b/R/run_mixder_ancestry.R
index 798ab2d..9e11c8e 100644
--- a/R/run_mixder_ancestry.R
+++ b/R/run_mixder_ancestry.R
@@ -28,23 +28,31 @@
#' @param A2_threshold Allele 2 probability threshold (default=0.60)
#' @param minor_contrib_threshold Whether to apply the allele 1 probability threshold to the minor contributor, regardless of the minimum number of SNPs (default=FALSE)
#' @param keep_bins Use existing binned SNP data, if exists (default=TRUE)
-#' @param snps SNPs to use for ancestry prediction (either ancestry only or all SNPs)
-#' @param pcagroups How to color PCA plots (superpopulations and/or subpopulations)
+#' @param snps SNPs to use for ancestry prediction (either ancestry only or all SNPs; default="ancestry")
+#' @param pcagroups How to color PCA plots (superpopulations and/or subpopulations, default="superpopulations")
+#' @param threads threads for EFM (default=0)
+#' @param parallel parallelize EFM sets, will detect number of cores on machine and use all available (default=FALSE)
#'
#' @export
#'
-run_mixder_ancestry = function(sample_manifest=NULL, sample=NULL, replicate=NULL, sample_reports = getwd(), output = "output", refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1_threshold=0.99, A2_threshold=0.6, minor_contrib_threshold=FALSE, keep_bins=TRUE, snps, pcagroups) {
+run_mixder_ancestry = function(sample_manifest=NULL, sample=NULL, replicate="", sample_reports = getwd(), output = "output", refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1_threshold=0.99, A2_threshold=0.6, minor_contrib_threshold=FALSE, keep_bins=TRUE, snps="ancestry", pcagroups="superpopulations", threads=0, parallel=FALSE) {
date = glue("{Sys.Date()}_{format(Sys.time(), '%H_%M_%S')}")
+ popFreq = list(mixder::popFreq_1000G, mixder::popFreq_1000G)
+ out_path = glue("{sample_reports}/snp_sets/{output}/")
## load in references
if (isTruthy(refpath)) {
print("loading references")
- if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
- refData = euroformix::sample_tableToList(data.frame(processing_ref_sample_reports(refpath)))
+ if (file.exists(glue("{refpath}/EFM_references.rda"))) {
+ load(glue("{refpath}/EFM_references.rda"))
} else {
- refData = euroformix::sample_tableToList(euroformix::tableReader(glue("{refpath}/EFM_references.csv")))
+ if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
+ refData = convert_table_to_list(data.frame(processing_ref_sample_reports(refpath)))
+ } else {
+ refsdf = data.frame(fread(glue("{refpath}/EFM_references.csv")))
+ refData = convert_table_to_list(refsdf)
+ }
+ save(refData, file=glue("{refpath}/EFM_references.rda"))
}
- } else if (cond) {
- stop("No references provided but selected conditioned analyses. Please re-run!")
} else {
refData = NULL
}
@@ -65,8 +73,9 @@ run_mixder_ancestry = function(sample_manifest=NULL, sample=NULL, replicate=NULL
if (!isTruthy(sample_manifest) & !isTruthy(sample)){
stop("Please provide sample manifest or sample ID!")
}
+ snp_positions = mixder::kintelligence_snp_positions
print("Running ancestry prediction")
- create_config(date, FALSE, "1000G_global", NA, NA, refpath, sample_manifest, sample, replicate, output, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, FALSE, FALSE, snpset, pcagroupcat)
+ create_config(date, FALSE, "1000G_global", NA, NA, refpath, sample_manifest, sample, replicate, out_path, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, FALSE, FALSE, snpset, pcagroupcat, "kintelligence", NULL)
if (isTruthy(sample_manifest)) {
## sample manifest; loop through each sample
manifest=suppressWarnings(euroformix::tableReader(sample_manifest))
@@ -75,12 +84,11 @@ run_mixder_ancestry = function(sample_manifest=NULL, sample=NULL, replicate=NULL
replicate_id = ifelse(is.na(manifest[i, 2]), "", manifest[i, 2])
message(glue("Sample ID: {id}"))
message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, id, replicate_id, FALSE, "1000G_global", NA, NA, refData, refpath, output, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, FALSE, FALSE, snpset, pcagroupcat)
+ run_workflow(date, id, replicate_id, FALSE, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, FALSE, FALSE, snpset, pcagroupcat, "kintelligence", snp_positions, threads, parallel)
}
} else if (isTruthy(sample)){
- replicate_id = ifelse(isTruthy(replicate), replicate, "")
message(glue("Sample ID: {sample}"))
- message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, sample, replicate_id, FALSE, "1000G_global", NA, NA, refData, refpath, output, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, FALSE, FALSE, snpset, pcagroupcat)
+ message(glue("Replicate ID: {replicate}"))
+ run_workflow(date, sample, replicate, FALSE, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, FALSE, FALSE, snpset, pcagroupcat, "kintelligence", snp_positions, threads, parallel)
}
}
diff --git a/R/run_mixder_metrics.R b/R/run_mixder_metrics.R
index aff5bce..9b4e6e0 100644
--- a/R/run_mixder_metrics.R
+++ b/R/run_mixder_metrics.R
@@ -37,21 +37,35 @@
#' @param filter_missing Whether to remove SNPs with a missing value for the allele 2 (default=FALSE)
#' @param minor_contrib_threshold Whether to apply the allele 1 probability threshold to the minor contributor, regardless of the minimum number of SNPs (default=FALSE)
#' @param keep_bins Use existing binned SNP data, if exists (default=TRUE)
+#' @param threads number of threads for EFM (default=0)
+#' @param parallel parallelize EFM sets, will detect number of cores on machine and use all available (default=FALSE)
#'
#' @export
#'
-run_mixder_metrics = function(sample_manifest=NULL, sample=NULL, replicate=NULL, sample_reports = getwd(), output = "output", twofreqs=FALSE, freq_both="global_1000g", freq_major=NULL, freq_minor=NULL, refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1min=0.95, A1max=0.99, A2min=0.5, A2max=0.65, major=NULL, minor=NULL, filter_missing=FALSE, minor_contrib_threshold=FALSE, keep_bins=TRUE) {
+run_mixder_metrics = function(sample_manifest=NULL, sample=NULL, replicate="", sample_reports = getwd(), output = "output", assay="kintelligence", twofreqs=FALSE, freq_both="global_1000g", freq_major=NULL, freq_minor=NULL, refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1min=0.95, A1max=0.99, A2min=0.5, A2max=0.65, major=NULL, minor=NULL, filter_missing=FALSE, minor_contrib_threshold=FALSE, keep_bins=TRUE, threads=0, parallel=FALSE) {
date = glue("{Sys.Date()}_{format(Sys.time(), '%H_%M_%S')}")
+ out_path = glue("{sample_reports}/snp_sets/{output}/")
if (!isTruthy(major) | !isTruthy(minor)) {
stop("Major and/or minor contributor not specified. Please rerun.")
}
+ popFreq = load_freq(twofreqs, freq_both, freq_major, freq_minor)
## load in references
+ print("loading references")
if (isTruthy(refpath)) {
- print("loading references")
- if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
- refData = euroformix::sample_tableToList(data.frame(processing_ref_sample_reports(refpath)))
+ if (cond) {
+ if (file.exists(glue("{refpath}/EFM_references.rda"))) {
+ load(glue("{refpath}/EFM_references.rda"))
+ } else {
+ if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
+ refData = convert_table_to_list(data.frame(processing_ref_sample_reports(refpath)))
+ } else {
+ refsdf = data.frame(fread(glue("{refpath}/EFM_references.csv")))
+ refData = convert_table_to_list(refsdf)
+ }
+ save(refData, file=glue("{refpath}/EFM_references.rda"))
+ }
} else {
- refData = euroformix::sample_tableToList(euroformix::tableReader(glue("{refpath}/EFM_references.csv")))
+ refData = NULL
}
} else {
stop("No references provided. Please re-run!")
@@ -60,7 +74,7 @@ run_mixder_metrics = function(sample_manifest=NULL, sample=NULL, replicate=NULL,
if (!isTruthy(sample_manifest) & !isTruthy(sample)){
stop("Please provide sample manifest or sample ID!")
}
- create_config(date, twofreqs, freq_both, freq_major, freq_minor, refpath, sample_manifest, sample, replicate, output, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, filter_missing, TRUE, NA, NA)
+ create_config(date, twofreqs, freq_both, freq_major, freq_minor, refpath, sample_manifest, sample, replicate, out_path, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, filter_missing, TRUE, NA, NA, assay, NULL)
if (isTruthy(sample_manifest)) {
manifest=suppressWarnings(euroformix::tableReader(sample_manifest))
for (i in nrow(manifest)) {
@@ -68,12 +82,11 @@ run_mixder_metrics = function(sample_manifest=NULL, sample=NULL, replicate=NULL,
replicate_id = ifelse(is.na(manifest[i, 2]), "", manifest[i, 2])
message(glue("Sample ID: {id}"))
message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, id, replicate_id, twofreqs, freq_both, freq_major, freq_minor, refData, refpath, output, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA)
+ run_workflow(date, id, replicate_id, twofreqs, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA, assay, NULL, threads, parallel)
}
} else if (isTruthy(sample)){
- replicate_id = ifelse(isTruthy(replicate), replicate, "")
message(glue("Sample ID: {sample}"))
- message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, sample, replicate_id, twofreqs, freq_both, freq_major, freq_minor, refData, refpath, output, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA)
+ message(glue("Replicate ID: {replicate}"))
+ run_workflow(date, sample, replicate, twofreqs, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "Calculate Metrics", sets, sample_reports, dynamicAT, staticAT, minimum_snps, NA, NA, A1min, A1max, A2min, A2max, major, minor, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA, assay, NULL, threads, parallel)
}
}
diff --git a/R/run_mixder_report.R b/R/run_mixder_report.R
index 3b5c0ab..80c5255 100644
--- a/R/run_mixder_report.R
+++ b/R/run_mixder_report.R
@@ -15,6 +15,8 @@
#' @param replicate If running individual sample along with a replicate, specify the replicate ID; requires sample_manifest=NULL (default=NULL)
#' @param sample_reports Path to directory of Sample Reports or CSV file of mixture genotypes (default=current directory)
#' @param output name of output directory, will be outputted into the sample reports directory (default="output")
+#' @param assay name of sequencing assay, either `kintelligence` or `custom` (default=kintelligence)
+#' @param snp_pos Path to file containing SNP positions of custom assay, required if assay==`custom` (default=NULL)
#' @param twofreqs TRUE if using different allele frequency data for each contributor (default=FALSE)
#' @param freq_both frequency data for both contributors; requires twofreqs=FALSE (default=1000G global)
#' @param freq_major frequency data for major contributor; requires twofreqs=TRUE (default=NULL)
@@ -33,18 +35,39 @@
#' @param filter_missing Whether to remove SNPs with a missing value for the allele 2 (default=FALSE)
#' @param minor_contrib_threshold Whether to apply the allele 1 probability threshold to the minor contributor, regardless of the minimum number of SNPs (default=FALSE)
#' @param keep_bins Use existing binned SNP data, if exists (default=TRUE)
+#' @param threads number of threads for EFM (default=0)
+#' @param parallel parallelize EFM sets, will detect number of cores on machine and use all available (default=FALSE)
#'
#' @export
#'
-run_mixder_report = function(sample_manifest=NULL, sample=NULL, replicate=NULL, sample_reports = getwd(), output = "output", twofreqs=FALSE, freq_both="global_1000g", freq_major=NULL, freq_minor=NULL, refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1_threshold=0.99, A2_threshold=0.6, filter_missing=FALSE, minor_contrib_threshold=FALSE, keep_bins=TRUE) {
+run_mixder_report = function(sample_manifest=NULL, sample=NULL, replicate="", sample_reports = getwd(), output = "output", assay="kintelligence", snp_pos=NULL, twofreqs=FALSE, freq_both="global_1000g", freq_major=NULL, freq_minor=NULL, refpath=NULL, refs=NULL, mixdeconv=TRUE, uncond=TRUE, cond=FALSE, sets=10, dynamicAT=0.015, staticAT=10, minimum_snps=6000, A1_threshold=0.99, A2_threshold=0.6, filter_missing=FALSE, minor_contrib_threshold=FALSE, keep_bins=TRUE, threads=0, parallel=FALSE) {
date = glue("{Sys.Date()}_{format(Sys.time(), '%H_%M_%S')}")
+ out_path = glue("{sample_reports}/snp_sets/{output}/")
+ popFreq = load_freq(twofreqs, freq_both, freq_major, freq_minor)
+ if (tolower(assay)=="kintelligence") {
+ snp_positions = mixder::kintelligence_snp_positions
+ } else if (tolower(assay)=="custom") {
+ if (file.exists(snp_pos)) {
+ snp_positions = euroformix::tableReader(snp_pos)
+ } else {
+ stop("Custom assay specified but no SNP positions file provided.")
+ }
+ } else {
+ stop("Problem with assay name. Please specify either 'kintelligence' or 'custom'.")
+ }
## load in references
if (isTruthy(refpath)) {
print("loading references")
- if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
- refData = euroformix::sample_tableToList(data.frame(processing_ref_sample_reports(refpath)))
+ if (file.exists(glue("{refpath}/EFM_references.rda"))) {
+ load(glue("{refpath}/EFM_references.rda"))
} else {
- refData = euroformix::sample_tableToList(euroformix::tableReader(glue("{refpath}/EFM_references.csv")))
+ if (!file.exists(glue("{refpath}/EFM_references.csv"))) {
+ refData = convert_table_to_list(data.frame(processing_ref_sample_reports(refpath)))
+ } else {
+ refsdf = data.frame(fread(glue("{refpath}/EFM_references.csv")))
+ refData = convert_table_to_list(refsdf)
+ }
+ save(refData, file=glue("{refpath}/EFM_references.rda"))
}
} else if (cond) {
stop("No references provided but selected conditioned analyses. Please re-run!")
@@ -54,7 +77,7 @@ run_mixder_report = function(sample_manifest=NULL, sample=NULL, replicate=NULL,
if (!isTruthy(sample_manifest) & !isTruthy(sample)){
stop("Please provide sample manifest or sample ID!")
}
- create_config(date, twofreqs, freq_both, freq_major, freq_minor, refpath, sample_manifest, sample, replicate, output, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, filter_missing, TRUE, NA, NA)
+ create_config(date, twofreqs, freq_both, freq_major, freq_minor, refpath, sample_manifest, sample, replicate, out_path, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, filter_missing, TRUE, NA, NA, assay, snp_pos)
if (isTruthy(sample_manifest)) {
manifest=suppressWarnings(euroformix::tableReader(sample_manifest))
for (i in nrow(manifest)) {
@@ -62,12 +85,11 @@ run_mixder_report = function(sample_manifest=NULL, sample=NULL, replicate=NULL,
replicate_id = ifelse(is.na(manifest[i, 2]), "", manifest[i, 2])
message(glue("Sample ID: {id}"))
message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, id, replicate_id, twofreqs, freq_both, freq_major, freq_minor, refData, refpath, output, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA)
+ run_workflow(date, id, replicate_id, twofreqs, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA, assay, snp_positions, threads, parallel)
}
} else if (isTruthy(sample)){
- replicate_id = ifelse(isTruthy(replicate), replicate, "")
message(glue("Sample ID: {sample}"))
- message(glue("Replicate ID: {replicate_id}"))
- run_workflow(date, sample, replicate_id, twofreqs, freq_both, freq_major, freq_minor, refData, refpath, output, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA)
+ message(glue("Replicate ID: {replicate}"))
+ run_workflow(date, sample, replicate, twofreqs, popFreq, refData, refpath, out_path, mixdeconv, uncond, refs, "Create GEDmatch PRO Report", sets, sample_reports, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, NA, NA, NA, NA, NA, NA, minor_contrib_threshold, keep_bins, filter_missing, TRUE, NA, NA, assay, snp_positions, threads, parallel)
}
}
diff --git a/R/run_workflow.R b/R/run_workflow.R
index 50e6a10..c4d515d 100644
--- a/R/run_workflow.R
+++ b/R/run_workflow.R
@@ -16,9 +16,7 @@
#' @param id Sample ID
#' @param replicate_id Sample ID of replicate, if specified
#' @param twofreqs TRUE if using separate AF data for major and minor contributors
-#' @param freq_both Path (or name) of allele frequency data if using same data for both
-#' @param freq_major Path (or name) of allele frequency data for major contributor
-#' @param freq_minor Path (or name) of allele frequency data for minor contributor
+#' @param popFreq List of properly formatted allele frequency data
#' @param refData Reference data (if available)
#' @param refs Path of reference genotype(s) file
#' @param output Name of output directory
@@ -45,6 +43,9 @@
#' @param skipancestry TRUE/FALSE whether to skip ancestry prediction
#' @param ancestrysnps SNPs to use for ancestry prediction (either ancestry only or all SNPs)
#' @param pcagroups How to color PCA plots (superpopulations and/or subpopulations)
+#' @param assay assay used (kintelligence or custom)
+#' @param positions SNP chromosomal positions
+#' @param threads number of threads to use for EFM
#'
#' @export
#'
@@ -53,8 +54,7 @@
#'@importFrom utils write.table write.csv read.table
#'@importFrom grDevices dev.off png
#'@importFrom methods show
-run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major, freq_minor, refData, refs, output, run_mixdeconv, unconditioned, cond, method, sets, kinpath, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, major, minor, minor_threshold, keep_bins, filter_missing, skipancestry, ancestrysnps, pcagroups) {
- out_path = glue("{kinpath}/snp_sets/{output}/")
+run_workflow = function(date, id, replicate_id, twofreqs, popFreq, refData, refs, out_path, run_mixdeconv, unconditioned, cond, method, sets, kinpath, dynamicAT, staticAT, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, major, minor, minor_threshold, keep_bins, filter_missing, skipancestry, ancestrysnps, pcagroups, assay, positions, threads=0, parallel=FALSE) {
if (replicate_id == "") {
logfile = file(glue("{out_path}config_log_files/{date}/run_log_{id}_{date}.txt"), open = "wt")
} else {
@@ -75,17 +75,20 @@ run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major,
message(glue("Replicate Sample: {replicate_id} "))
## run EFM
if (run_mixdeconv | !skipancestry) {
- message("Loading Frequency Data ")
- popFreq = load_freq(out_path, twofreqs, freq_both, freq_major, freq_minor)
- attable = process_kinreport(id, replicate_id, kinpath, dynamicAT, staticAT)
+ if (dynamicAT != 0) {
+ message("creating AT table ")
+ attable = process_kinreport(id, replicate_id, kinpath, dynamicAT, staticAT)
+ } else {
+ attable = staticAT
+ }
if (!skipancestry) {
- efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, cond, uncond=unconditioned, keep_bins)
+ efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, assay, cond, uncond=unconditioned, keep_bins, threads=threads, parallel=parallel)
efm_results_minor = efm_results_major
} else if (twofreqs) {
- efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, cond, uncond=unconditioned, keep_bins)
- efm_results_minor = run_efm(date, popFreq[[2]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, cond, uncond=unconditioned, keep_bins)
+ efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, assay, cond, uncond=unconditioned, keep_bins, threads=threads, parallel=parallel)
+ efm_results_minor = run_efm(date, popFreq[[2]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, assay, cond, uncond=unconditioned, keep_bins, threads=threads, parallel=parallel)
} else {
- efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, cond, uncond=unconditioned, keep_bins)
+ efm_results_major = run_efm(date, popFreq[[1]], refData, id, replicate_id, kinpath, out_path, attable, sets, skipancestry, assay, cond, uncond=unconditioned, keep_bins, threads=threads, parallel=parallel)
efm_results_minor = efm_results_major
}
}
@@ -126,21 +129,21 @@ run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major,
uncond_filename_major = glue("{write_path}/{id}/unconditioned/{id}_efm_output_unconditioned_major.tsv")
uncond_filename_minor = glue("{write_path}/{id}/unconditioned/{id}_efm_output_unconditioned_minor.tsv")
if (file.exists(uncond_filename_major)) {
- uncond_table_major = read.table(uncond_filename_major, header=T, sep="\t")
+ uncond_table_major = fread(uncond_filename_major, header=T, sep="\t")
} else {
stop(glue("{uncond_filename_major} does not exist. You may need to run EFM or check the correct SNP file input folder and Output folder are correct!"))
}
if (file.exists(uncond_filename_minor)) {
- uncond_table_minor = read.table(uncond_filename_minor, header=T, sep="\t")
+ uncond_table_minor = fread(uncond_filename_minor, header=T, sep="\t")
} else {
stop(glue("{uncond_filename_minor} does not exist. You may need to run EFM or check the correct SNP file input folder and Output folder are correct!"))
}
}
if (method == "Create GEDmatch PRO Report" | !skipancestry) {
message("Creating GEDmatch PRO report for major contributor in unconditioned analysis. ")
- major_report = create_gedmatchpro_report(write_path, uncond_table_major, major_c, "major", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing)
+ major_report = create_gedmatchpro_report(write_path, uncond_table_major, major_c, "major", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, positions)
message("Creating GEDmatch PRO report for minor contributor in unconditioned analysis. ")
- minor_report = create_gedmatchpro_report(write_path, uncond_table_minor, "C2", "minor", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing)
+ minor_report = create_gedmatchpro_report(write_path, uncond_table_minor, "C2", "minor", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, positions)
if (method == "Create GEDmatch PRO Report") {
write.table(major_report[[1]], glue("{write_path}/GEDMatchPROReports/{id}_uncond_major_{type}_GEDmatchPROReport.txt"), col.names=T, sep="\t", row.names=F, quote=F)
write.csv(major_report[[2]], glue("{write_path}/GEDMatchPROReports/Metrics/{id}_uncond_major_{type}_GEDmatchPROReport_Metrics.csv"), row.names=F, quote=F)
@@ -158,24 +161,24 @@ run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major,
write.table(minor_report[[1]], glue("{write_path}/{id}/unconditioned/{id}_uncond_minor_{type}_Inferred_Genotypes.txt"), col.names=T, sep="\t", row.names=F, quote=F)
ancestry_prediction(minor_report[[1]], glue("{write_path}/{id}/unconditioned/"), id, "unconditioned", "minor", ancestrysnps, pcagroups)
}
- major_report = create_gedmatchpro_report(write_path, uncond_table_major, major_c, "major", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing)
+ major_report = create_gedmatchpro_report(write_path, uncond_table_major, major_c, "major", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, positions)
write.table(major_report[[1]], glue("{write_path}/GEDMatchPROReports/{id}_uncond_major_{type}_GEDmatchPROReport.txt"), col.names=T, sep="\t", row.names=F, quote=F)
write.csv(major_report[[2]], glue("{write_path}/GEDMatchPROReports/Metrics/{id}_uncond_major_{type}_GEDmatchPROReport_Metrics.csv"), row.names=F, quote=F)
png(glue("{write_path}/GEDMatchPROReports/Metrics/{id}_uncond_major_{type}_GEDmatchPROReport_Allele1_Probabilities_Density_Plot.png"))
show(major_report[[3]])
dev.off()
message("Creating GEDmatch PRO report for minor contributor in unconditioned analysis. ")
- minor_report = create_gedmatchpro_report(write_path, uncond_table_minor, minor_c, "minor", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing)
+ minor_report = create_gedmatchpro_report(write_path, uncond_table_minor, minor_c, "minor", minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, positions)
write.table(minor_report[[1]], glue("{write_path}/GEDMatchPROReports/{id}_uncond_minor_{type}_GEDmatchPROReport.txt"), col.names=T, sep="\t", row.names=F, quote=F)
write.csv(minor_report[[2]], glue("{write_path}/GEDMatchPROReports/Metrics/{id}_uncond_minor_{type}_GEDmatchPROReport_metrics.csv"), row.names=F, quote=F)
png(glue("{write_path}/GEDMatchPROReports/Metrics/{id}_uncond_minor_{type}_GEDmatchPROReport_Allele1_Probabilities_Density_Plot.png"))
show(minor_report[[3]])
dev.off()
} else if (method == "Calculate Metrics") {
- major_ref = format_ref(refData, major, refs)
+ major_ref = format_ref(major, refs)
major_tables = suppressWarnings(process_efm_files(uncond_table_major, major_c, major_ref, minimum_snps, A1min, A1max, A2min, A2max, metrics=TRUE, filter_missing))
write_tables(major_tables, glue("{write_path}/{id}/unconditioned/{major}"), minimum_snps)
- minor_ref = format_ref(refData, minor, refs)
+ minor_ref = format_ref(minor, refs)
minor_tables = suppressWarnings(process_efm_files(uncond_table_minor, minor_c, minor_ref, minimum_snps, A1min, A1max, A2min, A2max, metrics=TRUE, filter_missing))
write_tables(minor_tables, glue("{write_path}/{id}/unconditioned/{minor}"), minimum_snps)
}
@@ -212,7 +215,7 @@ run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major,
}
if (method == "Create GEDmatch PRO Report" | !skipancestry) {
message(glue("Creating GEDmatch PRO report for {contrib_status} contributor conditioned on {cond_on} in conditioned analysis. "))
- cond_report = create_gedmatchpro_report(write_path, get(glue("efm_table_{contrib_status}")), "C2", contrib_status, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing)
+ cond_report = create_gedmatchpro_report(write_path, get(glue("efm_table_{contrib_status}")), "C2", contrib_status, minimum_snps, A1_threshold, A2_threshold, A1min, A1max, A2min, A2max, minor_threshold, filter_missing, positions)
if (method == "Create GEDmatch PRO Report") {
write.table(cond_report[[1]], glue("{write_path}/GEDMatchPROReports/{id}_{contrib_status}_contrib_conditioned_on_{cond_on}_{type}_GEDmatchPROReport.txt"), col.names=T, sep="\t", row.names=F, quote=F)
write.csv(cond_report[[2]], glue("{write_path}/GEDMatchPROReports/Metrics/{id}_{contrib_status}_contrib_conditioned_on_{cond_on}_{type}_GEDmatchPROReport_Metrics.csv"), row.names=F, quote=F)
@@ -225,7 +228,7 @@ run_workflow = function(date, id, replicate_id, twofreqs, freq_both, freq_major,
}
} else if (method == "Calculate Metrics") {
unk = ifelse(contrib_status == "major", major, minor)
- ref = format_ref(refData, unk, refs)
+ ref = format_ref(unk, refs)
geno_correct_tables = suppressWarnings(process_efm_files(get(glue("efm_table_{contrib_status}")), "C2", ref, minimum_snps, A1min, A1max, A2min, A2max, metrics=TRUE, filter_missing))
write_tables(geno_correct_tables, glue("{write_path}/{id}/conditioned/{unk}"), minimum_snps)
}
diff --git a/README.Rmd b/README.Rmd
index 8dd2021..04ecd6e 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -13,13 +13,13 @@ knitr::opts_chunk$set(
)
```
-# MixDeR - Current Version: 0.9.0
+# MixDeR - Current Version: 0.10.0
-MixDeR (**Mix**ture **De**convolution in **R**) is a workflow (with a Shiny app) for performing mixture deconvolution of ForenSeq Kintelligence SNP data for two-person mixtures using [EuroForMix](https://github.com/oyvble/euroformix/) and creating GEDmatch PRO reports for the individual contributor SNP profiles.
+MixDeR (**Mix**ture **De**convolution in **R**) is a workflow (with a Shiny app) for performing mixture deconvolution of genotypes developed from either the ForenSeq Kintelligence assay or a custom SNP panel. MixDeR performs the deconvolution of two-contributor mixtures using [EuroForMix](https://github.com/oyvble/euroformix/), ultimately developing single source GEDmatch PRO/FamilyTreeDNA reports of the individual contributor genotypes.
This method requires extensive validation of the settings. MixDeR provides the option of calculating various metrics for evaluating the accuracy of the deduced SNP genotypes. This is extremely useful when determining settings, specifically the allele 1 probability threshold, the allele 2 probability threshold, and the minimum number of SNPs.
@@ -44,6 +44,6 @@ library(mixder)
mixder()
```
-To use the command line, please see the CLI page in the MixDeR documentation.
+To use the command line, please see the R API page in the MixDeR documentation.
Please see [MixDeR documentation](bioforensics.github.io/mixder/) for further information.
diff --git a/README.md b/README.md
index c87ee24..9437c97 100644
--- a/README.md
+++ b/README.md
@@ -1,17 +1,19 @@
-# MixDeR - Current Version: 0.9.0
+# MixDeR - Current Version: 0.10.0
MixDeR (**Mix**ture **De**convolution in **R**) is a workflow (with a
-Shiny app) for performing mixture deconvolution of ForenSeq
-Kintelligence SNP data for two-person mixtures using
-[EuroForMix](https://github.com/oyvble/euroformix/) and creating
-GEDmatch PRO reports for the individual contributor SNP profiles.
+Shiny app) for performing mixture deconvolution of genotypes developed
+from either the ForenSeq Kintelligence assay or a custom SNP panel.
+MixDeR performs the deconvolution of two-contributor mixtures using
+[EuroForMix](https://github.com/oyvble/euroformix/), ultimately
+developing single source GEDmatch PRO/FamilyTreeDNA reports of the
+individual contributor genotypes.
This method requires extensive validation of the settings. MixDeR
provides the option of calculating various metrics for evaluating the
@@ -43,7 +45,7 @@ To launch the shiny app:
library(mixder)
mixder()
-To use the command line, please see the CLI page in the MixDeR
+To use the command line, please see the R API page in the MixDeR
documentation.
Please see [MixDeR documentation](bioforensics.github.io/mixder/) for
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 92d28ef..3ed0337 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -10,6 +10,7 @@ articles:
- Installation
- Required_Files
- Running_MixDecon
+ - Custom_SNPs
- Ancestry_Prediction_Tool
- Valmetrics
- gedmatch
diff --git a/docs/404.html b/docs/404.html
index 0bf73ce..7f05c11 100644
--- a/docs/404.html
+++ b/docs/404.html
@@ -20,7 +20,7 @@
mixder
- 0.9.0
+ 0.10.0
NOTE: There are several allele frequency data files
-stored within MixDeR for your use. 1000G global is the default when
-using the same frequency file for both contributors
-(twofreqs=FALSE). The following datasets are available for
-the freq_both, freq_major and
+stored within MixDeR for your use for Kintelligence data ONLY.
+1000G global is the default when using the same frequency file for both
+contributors (twofreqs=FALSE). The following datasets are
+available for the freq_both, freq_major and
freq_minor arguments and should be referenced as
indicated:
- 1000 Genomes global (1000g_global)
@@ -102,8 +103,9 @@
R API
- 1000 Genomes EAS (eas)
- 1000 Genomes EUR (eur)
- 1000 Genomes SAS (sas)
-Of course, a custom allele frequency file is allowed, requiring only the
-Path to the file.
+Of course, a custom allele frequency file is also an option for
+Kintelligence data but is a requirement for a custom SNP set. The Path
+to the file must be specified if using a custom frequency dataset.
NOTE x 2: Using the R API, there are two options to
specify samples:
@@ -123,6 +125,9 @@
Running
Full reference for this function, including default settings,
can be found
here.
+
The following examples are specifically for the Kintelligence
+data; please see the Custom SNP set page
+for required arguments if using a custom SNP panel.
The following function runs an unconditioned mixture deconvolution
and creates a report using the default settings (See Running Mixture Deconvolution for more
information about the defaults), using the 1000 Genomes Global allele
@@ -157,6 +162,9 @@
Runnin
Full reference for this function, including default settings,
can be found
here.
+
The following examples are specifically for the Kintelligence
+data; please see the Custom SNP set page
+for required arguments if using a custom SNP panel.
The options are similar as the run_mixder_report
function. The following function runs an unconditioned mixture
deconvolution and calculates metrics using the default settings (See Running Mixture Deconvolution for more
diff --git a/docs/articles/Required_Files.html b/docs/articles/Required_Files.html
index 50b62bd..ab41c94 100644
--- a/docs/articles/Required_Files.html
+++ b/docs/articles/Required_Files.html
@@ -20,7 +20,7 @@
mixder
- 0.9.0
+ 0.10.0
@@ -37,6 +37,7 @@
@Manual{,
title = {mixder: A workflow for performing SNP mixture deconvolution},
author = {Rebecca Mitchell},
year = {2026},
- note = {R package version 0.9.0},
+ note = {R package version 0.10.0},
url = {https://bioforensics.github.io/mixder},
}
diff --git a/docs/conf.py b/docs/conf.py
deleted file mode 100644
index 94f19d9..0000000
--- a/docs/conf.py
+++ /dev/null
@@ -1,69 +0,0 @@
-# Configuration file for the Sphinx documentation builder.
-#
-# This file only contains a selection of the most common options. For a full
-# list see the documentation:
-# https://www.sphinx-doc.org/en/master/usage/configuration.html
-
-# -- Path setup --------------------------------------------------------------
-
-# If extensions (or modules to document with autodoc) are in another directory,
-# add these directories to sys.path here. If the directory is relative to the
-# documentation root, use os.path.abspath to make it absolute, like shown here.
-#
-import os
-import sys
-sys.path.insert(0, os.path.abspath('.'))
-
-
-# -- Project information -----------------------------------------------------
-
-project = 'MixDeR'
-copyright = '2025, DHS'
-author = 'NBFAC'
-
-# The full version, including alpha/beta/rc tags
-release = '1.0.0'
-
-
-# -- General configuration ---------------------------------------------------
-
-# Add any Sphinx extension module names here, as strings. They can be
-# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
-# ones.
-extensions = [
- 'myst_parser',
- 'sphinx.ext.duration',
- 'sphinx.ext.doctest',
- 'sphinx.ext.autodoc',
- 'sphinx.ext.autosummary',
- 'sphinx.ext.intersphinx',
-]
-
-intersphinx_mapping = {
- 'python': ('https://docs.python.org/3/', None),
- 'sphinx': ('https://www.sphinx-doc.org/en/master/', None),
-}
-intersphinx_disabled_domains = ['std']
-
-# Add any paths that contain templates here, relative to this directory.
-templates_path = ['_templates']
-
-# List of patterns, relative to source directory, that match files and
-# directories to ignore when looking for source files.
-# This pattern also affects html_static_path and html_extra_path.
-exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
-
-
-# -- Options for HTML output -------------------------------------------------
-
-# The theme to use for HTML and HTML Help pages. See the documentation for
-# a list of builtin themes.
-#
-#html_theme = 'alabaster'
-html_theme = 'sphinx_rtd_theme'
-
-
-# Add any paths that contain custom static files (such as style sheets) here,
-# relative to this directory. They are copied after the builtin static files,
-# so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
diff --git a/docs/index.html b/docs/index.html
index 98ee94e..72b7f21 100644
--- a/docs/index.html
+++ b/docs/index.html
@@ -22,7 +22,7 @@
mixder
- 0.9.0
+ 0.10.0
@@ -39,6 +39,7 @@
MixDeR (Mixture Deconvolution in R) is a workflow (with a Shiny app) for performing mixture deconvolution of ForenSeq Kintelligence SNP data for two-person mixtures using EuroForMix and creating GEDmatch PRO reports for the individual contributor SNP profiles.
+
MixDeR (Mixture Deconvolution in R) is a workflow (with a Shiny app) for performing mixture deconvolution of genotypes developed from either the ForenSeq Kintelligence assay or a custom SNP panel. MixDeR performs the deconvolution of two-contributor mixtures using EuroForMix, ultimately developing single source GEDmatch PRO/FamilyTreeDNA reports of the individual contributor genotypes.
This method requires extensive validation of the settings. MixDeR provides the option of calculating various metrics for evaluating the accuracy of the deduced SNP genotypes. This is extremely useful when determining settings, specifically the allele 1 probability threshold, the allele 2 probability threshold, and the minimum number of SNPs.
Note: MixDeR (and EFM) assume the mixture samples are composed of two contributors. MixDeR is able to identify and alert the user to samples that may be potentially either single source or consist of a mixture with a large mixture ratio (i.e. > 1:100 ratio between contributors). In these scenarios, the user is warned to be cautious with the minor contributor inferred genotyping results.
MixDeR version 0.8 and later provides the option to perform ancestry prediction using principal component analysis (PCA). Additional information can be found below in the Ancestry Prediction section.
@@ -85,7 +86,7 @@
Usage
To launch the shiny app:
library(mixder)mixder()
-
To use the command line, please see the CLI page in the MixDeR documentation.
+
To use the command line, please see the R API page in the MixDeR documentation.
diff --git a/docs/index.rst b/docs/index.rst
deleted file mode 100644
index e094f71..0000000
--- a/docs/index.rst
+++ /dev/null
@@ -1,81 +0,0 @@
-.. documentation master file, created by sphinx-quickstart
- You can adapt this file completely to your liking, but it should at least
- contain the root `toctree` directive.
-
-reStructuredText
-================================
-
-.. raw:: html
-
-
-
-.. role:: red
-
-This main document is in `'reStructuredText' ("rst") format
-`_,
-which differs in many ways from standard markdown commonly used in R packages.
-``rst`` is richer and more powerful than markdown. The remainder of this main
-document demonstrates some of the features, with links to additional ``rst``
-documentation to help you get started. The definitive argument for the benefits
-of ``rst`` over markdown is the `official language format documentation
-`_, which starts with a very clear
-explanation of the `benefits
-`_.
-
-Examples
---------
-
-All of the following are defined within the ``docs/index.rst`` file. Here is
-some :red:`coloured` text which demonstrates how raw HTML commands can be
-incorporated. The following are examples of ``rst`` "admonitions":
-
-.. note::
-
- Here is a note
-
- .. warning::
-
- With a warning inside the note
-
-.. seealso::
-
- The full list of `'restructuredtext' directives `_ or a similar list of `admonitions `_.
-
-.. centered:: This is a line of :red:`centered text`
-
-.. hlist::
- :columns: 3
-
- * and here is
- * A list of
- * short items
- * that are
- * displayed
- * in 3 columns
-
-The remainder of this document shows three tables of contents for the main
-``README`` (under "Introduction"), and the vignettes and R directories of
-a package. These can be restructured any way you like by changing the main
-``docs/index.rst`` file. The contents of this file -- and indeed the contents
-of any `readthedocs `_ file -- can be viewed by
-clicking *View page source* at the top left of any page.
-
-.. toctree::
- :maxdepth: 1
- :caption: Introduction:
-
- mixder.md
-
-
-
-
-.. toctree::
- :maxdepth: 1
- :caption: Vignettes
-
-
-
-.. toctree::
- :maxdepth: 1
- :caption: Functions
-
diff --git a/docs/make.bat b/docs/make.bat
deleted file mode 100644
index 8084272..0000000
--- a/docs/make.bat
+++ /dev/null
@@ -1,35 +0,0 @@
-@ECHO OFF
-
-pushd %~dp0
-
-REM Command file for Sphinx documentation
-
-if "%SPHINXBUILD%" == "" (
- set SPHINXBUILD=sphinx-build
-)
-set SOURCEDIR=.
-set BUILDDIR=_build
-
-if "%1" == "" goto help
-
-%SPHINXBUILD% >NUL 2>NUL
-if errorlevel 9009 (
- echo.
- echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
- echo.installed, then set the SPHINXBUILD environment variable to point
- echo.to the full path of the 'sphinx-build' executable. Alternatively you
- echo.may add the Sphinx directory to PATH.
- echo.
- echo.If you don't have Sphinx installed, grab it from
- echo.https://www.sphinx-doc.org/
- exit /b 1
-)
-
-%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
-goto end
-
-:help
-%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O%
-
-:end
-popd
diff --git a/docs/mixder.md b/docs/mixder.md
deleted file mode 100644
index fdf0968..0000000
--- a/docs/mixder.md
+++ /dev/null
@@ -1,352 +0,0 @@
-
-
-# MixDeR - Current Version: 1.0
-
-
-
-
-MixDeR (**Mix**ture **De**convolution in **R**) is a workflow (with a
-Shiny app) for performing mixture deconvolution of ForenSeq
-Kintelligence SNP data for two-person mixtures using
-[EuroForMix](https://github.com/oyvble/euroformix/) and creating
-GEDmatch PRO reports for the individual contributor SNP profiles.
-
-This method requires extensive validation of the settings. MixDeR
-provides the option of calculating various metrics for evaluating the
-accuracy of the deduced SNP genotypes. This is extremely useful when
-determining settings, specifically the allele 1 probability threshold,
-the allele 2 probability threshold, and the minimum number of SNPs.
-
-*Note: MixDeR (and EFM) assume the mixture samples are composed of two
-contributors. MixDeR is able to identify and alert the user to samples
-that may be potentially either single source or consist of a mixture
-with a large mixture ratio (i.e. \> 1:100 ratio between contributors).
-In these scenarios, the user is warned to be cautious with the minor
-contributor inferred genotyping results.*
-
-Please cite the following paper if utilizing MixDeR:
-```
-Mitchell, R., Peck, M., Gorden, E., & Just, R. (2025). MixDeR: A SNP mixture
-deconvolution workflow for forensic genetic genealogy. Forensic Science
-International: Genetics, 76, 103224, doi: 10.1016/j.fsigen.2025.103224
-```
-
-
-## Installation
-
-For any installation, EuroForMix must be installed. Please follow the
-instructions from the [EuroForMix GitHub
-page](https://github.com/oyvble/euroformix/). EuroForMix version 4.2.4
-and earlier have been tested and are compatible with MixDeR. If using a
-newer version, please be aware it has not been tested and errors may
-occur!
-
-If installing from GitHub:
-The R package `devtools` is required to install from GitHub:
-
- install.packages("devtools")
- devtools::install_github("bioforensics/mixder")
-
-If installing from source, first install the following R packages:
-
- install.packages(c("dplyr", "ggplot2", "glue", "prompter", "readxl", "rlang", "shiny", "shinyFiles", "shinyjs", "tibble", "tidyr"))
-
-To install MixDeR from source (i.e. the `mixder_0.1.0.tar.gz` file):
-
- install.packages("/path/to/mixder_0.1.0.tar.gz", repos = NULL, type="source")
-
-For example, if the `mixder_0.1.0.tar.gz` file is located in your
-Documents folder:
-
- install.packages("~/Documents/mixder_0.1.0.tar.gz", repos = NULL, type="source")
-
-## Usage
-
-To launch the shiny app:
-
- library(mixder)
- mixder()
-
-## Required files
-
-1. **Mixture Kintelligence SNP profiles**
- This can be in the form of the UAS Sample Report or a TSV file with
- the below format:
-
-| Marker | Allele | Reads |
-|:----------:|:------:|:-----:|
-| rs12615742 | T | 0 |
-| rs12615742 | C | 134 |
-| rs16885694 | G | 43 |
-| rs16885694 | A | 63 |
-
-The files should be tab delimited and should be named as a `.tsv` file,
-such as: `SampleID.tsv`.
-
-If using the Shiny app, you must specify the folder containing these SNP
-files. Multiple samples (with multiple files each) can be in the same
-folder. Additional files may be present in the folder and will be
-ignored by MixDeR.
-
-MixDeR will divide the entire Kintelligence dataset into more manageable
-sets (organized by total SNP read depth) to run through EFM (ideal for
-best performance). The user may specify how many sets the program will
-use (see below); the default is 10 sets. The user must then specify how
-many sets are provided so MixDeR knows how many files to process per
-sample.
-
-The default is for MixDeR to use previously-created SNP sets, if present
-in the specified input folder. If this option is unselected, MixDeR will
-create new SNP sets files, overwriting any previously made files.
-
-2. **The sample manifest**
- This file lists the Sample IDs of the files to run. The `SampleID`
- is extracted from the `Sample Name` field in the Sample Report. The
- columns names do not matter, just the order and both columns
- **must** be present even if no replicates are included. If a single
- sample is run, only the single ID needs to be in the first column,
- the second column should be left blank. If a second sample is to be
- run in replicate, the replicate ID should be listed in the second
- column.
-
-| SampleID | ReplicateID |
-|:---------:|:-----------:|
-| Sample01a | |
-| Sample01a | Sample01b |
-
-## Other files which may be required
-
-**Allele frequency file**
-
-MixDeR provides general population allele frequencies for Kintelligence
-SNPs from either [1000 Genomes Phase 3
-dataset](https://www.internationalgenome.org/home) or the [gnomAD v4
-dataset](https://gnomad.broadinstitute.org/downloads#v4). However, it is
-ideal to use allele frequencies derived from the population that closely
-matches the contributor(s) of interest. Therefore, MixDeR provides the
-user the opportunity to upload a different allele frequency file.
-EuroForMix requires the below format for allele frequency files, for all
-SNPs with each SNP as its own column:
-
-| Allele | rs6690515 | rs424079 | rs2055204 |
-|:------:|:---------:|:--------:|:---------:|
-| A | 0.122837 | 0.64677 | 0.501441 |
-| C | NA | 0.35323 | NA |
-| G | 0.877163 | NA | 0.498559 |
-| T | NA | NA | NA |
-
-Given the difficulty of formatting the data as such, MixDeR will create
-this format for the user if the user provides the frequency data in a
-CSV file with the following format (NOTE: the column names MUST match
-below; the order of columns and additional columns will not affect it).
-
-| SNP | Ref | Alt | Ref_AF | Alt_AF |
-|:---------:|:---:|:---:|:--------:|:--------:|
-| rs6690515 | A | C | 0.35323 | 0.64677 |
-| rs424079 | T | C | 0.122837 | 0.877163 |
-| rs2055204 | G | A | 0.501441 | 0.498559 |
-
-MixDeR provides the option to use the same allele frequency file for
-both the major and minor contributor or select different allele
-frequency files for the major and minor contributor. If two different
-frequency files are selected, MixDeR will run EFM twice, once using the
-allele frequency file for the major contributor (and extracting the
-inferred genotypes for the major contributor from those results) and
-once using the allele frequency file for the minor contributor (and
-extracting the inferred genotypes for the minor contributor from those
-results).
-
-**Reference Genotypes**
-
-If calculating validation metrics or performing a conditioned
-deconvolution, the reference genotypes are required.
-
-There are two options for providing reference genotypes. MixDeR accepts
-the UAS Sample Report, stored in a separate folder. A second option is
-to provide a single CSV file containing all references named
-`EFM_references.csv` with the following format:
-
-| Sample.Name | Marker | Allele1 | Allele2 |
-|:-----------:|:----------:|:-------:|:-------:|
-| Sample01 | rs12615742 | T | C |
-| Sample01 | rs16885694 | G | G |
-
-The user must provide the folder containing either the Sample Reports or
-the CSV reference file.
-NOTE: MixDeR first searches for the CSV file in the provided folder. If
-there are additional wanted references not contained in this file,
-please remove the CSV file and run MixDeR again. MixDeR creates the CSV
-file containing genotypes from the Sample Reports within the provided
-folder.
-\_\_\_\_\_\_\_\_
-
-## MixDeR Details
-
-MixDeR has four modules:
-1. Ancestry Prediction
-1. EFM mixture deconvolution
-2. Calculate validation metrics
-3. Create GEDmatch PRO reports
-
-EFM mixture deconvolution must be run at least once. If it’s been run
-previously, the other modules can be run using the existing
-deconvolution data.
-
-MixDeR can either calculate validation metrics OR create the GEDmatch
-PRO reports during a single run, not both.
-
-NOTE about the allele probability thresholds:
-This workflow utilizes individual probabilities for each allele call
-from EFM. The reported allele 1 is the allele with the higher
-probability; allele 2 is the allele with the lower probability. When
-applying the allele 1 probability threshold, any SNP with a probability
-below the threshold will be removed completely from the dataset. When
-applying the allele 2 probability threshold to the remaining SNPs, if
-the probability is below the threshold, allele 2 is reported as the same
-allele as allele 1. If it is above, the allele will be reported as
-called. For example, if the genotype for SNP rs12615742 is C,T but the
-allele 2 probability is below the threshold, the SNP genotype will be
-reported as C,C. If it is above the threshold, the SNP genotype will be
-reported as C,T.
-
-## Ancestry Prediction - Background
-The mixture deconvolution algorithm requires allele frequency data. While we’ve found that using global allele frequencies generally performs quite well, using allele frequencies as closely matched to the population group of the contributor can result in increased accuracy.
-
-Given forensic samples are almost always from individuals of unknown ancestry, we’ve developed a method using the inferred genotypes from the mixture deconvolution method to predict the ancestry of each contributor using principal component analysis (PCA).
-
-PCA uses the genotypes from individuals of known ancestry, in this case the 2,504 1000 Genomes samples, to create a statistical framework for predicting the ancestry of the unknown sample. PCA transforms high dimensionality data (here, the genomic data) into a lower-dimensionality space in the form of principal components, or PCs. PC1 explains the highest amount of variance in the data, PC2 explains the second highest amount of variance, and so forth. Using genetic data, the biogeographical ancestry is driving the top PCs given it accounts for most of the variation in the data.
-
-Plotting the PCs against each other (and coloring samples by population) allows one to visually see the separation and clustering of populations, oftentimes down to the subpopulation level. Where the unknown sample falls within the plot (along with some distance calculations) allows the user to assign a population to the sample, assuming it falls within a known cluster.
-
-Generally, ancestry can be visualized by plotting PC1 vs. PC2. However, we found that adding the 3rd PC in a 3-dimensional space provides the best separation of the populations. MixDeR creates 3-D ancestry plots and saves them as a .html file. These interactive plots are more insightful than the standard 2-D PCA plots.
-
-MixDeR calculates the Euclidean distance from the centroid of each superpopulation to the unknown sample and rank the superpopulations based on the distance with the smallest distance on top. It should be noted here that the top ranked superpopulation does NOT always match the actual ancestry of the unknown but instead lists the closest population distance-wise. The user must consider the actual value of that distance and examine the accompanying PCA plot to determine the potential accuracy of the prediction. If the unknown sample is falling outside of the top-ranked population, it likely does not closely match that population. While this can be a reflection of the quality of the sample and/or PCA, it could also occur because the unknown ancestry does not match those in the 1000 Genomes database. For example, Ashkenazi Jews are not included in the 1000 Genomes dataset and given their genetic divergence from other European populations, would fall outside of the EUR superpopulation (at least that’s what we’ve found!). All of these factors must be considered when evaluating the ancestry prediction results.
-
-MixDeR does not make its own ancestry determination, but provides the data and plots for the user to make the determination themselves. Analyzing samples of known ancestry is pertinent to understanding the strengths and weakness of the method and the developers strongly encourage users to perform validation studies before making any ancestry predictions.
-
-## Running the Ancestry Prediction tool
-
-When first launching the Shiny app, the ancestry prediction tool first appears. This tool can be skipped by checking the ```Skip Ancestry Prediction Step``` box.
-
-The user must select whether they want to use the Superpopulation groups (AFR, AMR, EAS, EUR, and SAS) and/or the subpopulation groups (e.g. Toscani in Italy, Esan in Nigeria, etc.). If both are selected, PCA plots (coloring by either superpopulation or subpopulation) and centroid calculations are performed for both.
-The ancestry prediction tool will first perform mixture deconvolution using the 1000 Genomes Global allele frequency data and perform genotype filtering using the specified settings (allele 1/2 probability thresholds, # of SNP bins, the static/dynamic AT, minimum # of SNPs). A conditioned and/or unconditioned deconvolution can be performed. A sample manifest and folder containing the mixture Sample Reports (as well the reference Sample Reports, if performing a conditioned deconvolution) are required. Please see the section entitled “Running mixture deconvolution using EuroForMix†for further information and guidance on these settings.
-
-The user can select whether to use the only 94 ancestry SNPs or all autosomal SNPs for PCA (the deconvolution step is performed using all SNPs). While the 94 ancestry SNPs do an ample job with good quality samples and is quite fast, using all SNPs provides better clustering but is significantly slower (several minutes per contributor). It is important the user weights the benefits vs. drawbacks of each set of SNPs to determine the best option for their data and system.
-
-
-## Running mixture deconvolution using EuroForMix
-
-There are several options/settings to run EFM mixture deconvolution:
-The type of mixture deconvolution analysis to perform (one or both can
-be selected at once):
-- Unconditioned analysis
-- Conditioned analysis
-
-**Allele Frequency Data**: The user must choose which allele frequency
-data to use: 1000Genomes Phase 3 data, gnomAD v4 data, or upload a
-custom file. See above for more details about the format for uploading a
-custom AF file.
-**References to Condition on**: IF running a conditioned analysis, once
-the reference folder has been uploaded, this dropdown menu will
-auto-populate containing the reference sample IDs. The user must select
-which references to condition on. **As of now, MixDeR can only condition
-on a single reference sample**. However, multiple references can be
-selected for conditioning; the conditioned analyses will be run
-separately.
-**Number of SNP Bins**: The number of SNP sets for each sample. The
-default is 10.
-**Static Analytical Threshold**: The minimum number of reads required
-for a SNP to be included (default = 10).
-**Dynamic Analytical Threshold**: The percent of total SNP reads
-required for a SNP to. be included (default = 0.015 or 1.5%). (A quick
-note on the ATs: both the static and dynamic ATs are applied and the one
-producing the higher AT will be used in the EFM software).
-**Output Folder Name**: The name of the folder created within the folder
-containing the original SNP files to store the generated output for the
-entire workflow.
-**Minimum Number of SNPs**: The minimum number of SNPs
-
-## Calculating validation metrics
-
-There are several options/settings to calculate the validation
-metrics:
-**Minimum Allele 1 Probability Threshold** and **Maximum Allele 1
-Probability Threshold**: The range of allele 1 probability thresholds
-for calculating the validation metrics, increasing in increments of 0.01
-(i.e. if minimum is set to 0.5 and maximum is set to 1, will calculate
-metrics using a threshold of 0.5, 0.51, 0.52, 0.53, up to 1).
-**Minimum Allele 2 Probability Threshold** and **Maximum Allele 2
-Probability Threshold**: The range of allele 2 probability thresholds
-for calculating the validation metrics, increasing in increments of 0.01
-(i.e. if minimum is set to 0.5 and maximum is set to 1, will calculate
-metrics using a threshold of 0.5, 0.51, 0.52, 0.53, up to 1).
-**Major Contributor Sample ID**: The ID of the major contributor of the
-mixture. Once the reference folder is uploaded, this dropdown menu will
-auto-populate with the reference sample IDs.
-**Minor Contributor Sample ID**: The ID of the minor contributor of the
-mixture. Once the reference folder is uploaded, this dropdown menu will
-auto-populate with the reference sample IDs.
-**Remove SNPs If Missing Either Allele?**: If an allele 1 is inferred to
-be missing (reported as `99`), that SNP will be automatically dropped
-from the final dataset. By default, if an allele 2 is inferred to be
-missing and the allele 2 probability is above the allele 2 probability
-threshold, the SNP is reported as homozygous for allele 1. However,
-selecting this option will result in dropping the SNP if the allele 2
-probability of the missing allele 2 is above the allele 2 probability
-threshold, instead of reporting the SNP as homozygous for allele 1.
-
-If calculating validation metrics, reference genotypes are required to
-calculate genotype accuracy. MixDerR will calculate the metrics for the
-range of allele 1 probability thresholds and allele 2 probability
-thresholds specified by the user. The final output file looks as such:
-
-| A1 cutoff | A2 cutoff | Total SNPs | N No Ref | N SNPs tested | N Genotypes Correct | Genotype Accuracy | Heterozygosity |
-|:--:|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
-| 0.99 | 0.01 | 9735 | 8 | 9727 | 9549 | 0.9817 | 0.456 |
-| 0.99 | 0.02 | 9735 | 8 | 9727 | 9548 | 0.9815 | 0.456 |
-| 0.99 | 0.03 | 9735 | 8 | 9727 | 9548 | 0.9815 | 0.456 |
-
-## Creating GEDmatch PRO Reports
-
-**Allele 1 Probability Threshold to create GEDmatch PRO Report** and
-**Allele 2 Probability Threshold to create GEDmatch PRO Report**: The
-allele 1 and allele 2 probability thresholds. If the contributor is the
-major contributor to the mixture, MixDeR will apply the allele 1 and
-allele 2 probability thresholds UNLESS this results in a profile with
-fewer SNPs than the specified minimum (6000 is the default). If it does
-not meet this minimum, the top 6,000 SNPs (or whatever the user
-specifies as the minimum) with the highest allele 1 probabilities are
-used and then the allele 2 probability threshold is applied. For minor
-contributors, the default setting is that the minimum number of SNPs is
-automatically used and then the allele 2 probability threshold is
-applied. The option to apply the allele 1 probability threshold (similar
-in manner to the major contributor) can be applied.
-**Remove SNPs If Missing Either Allele?**: If an allele 1 is inferred to
-be missing (reported as `99`), that SNP will be automatically dropped
-from the final dataset. By default, if an allele 2 is inferred to be
-missing and the allele 2 probability is above the allele 2 probability
-threshold, the SNP is reported as homozygous for allele 1. However,
-selecting this option will result in dropping the SNP if the allele 2
-probability of the missing allele 2 is above the allele 2 probability
-threshold, instead of reporting the SNP as homozygous for allele 1.
-
-As a way to assist the analyst in evaluating the inferred genotypes of a
-mixture of unknown contributors, several metrics are calculated in this
-step for three different scenarios: (1) only the allele 2 probability
-threshold applied; (2) the allele 1 and allele 2 probability thresholds
-applied; and (3) the minimum number of SNPs used and the allele 2
-probability threshold applied. For each dataset, the follow metrics are
-calculated: number of SNPs, mean allele 1 probabilities, median allele 1
-probabilities, standard deviation (SD) of the allele 1 probabilities and
-heterozygosity. Below is an example of the table created by MixDeR in
-this step. A density plot of allele 1 probability thresholds is also
-created. In general, the more SNPs with higher allele 1 probabilities,
-the higher the accuracy of the inferred genotypes. However, determining
-exactly what qualifies as acceptable should be determined by individual
-labs.
-
-| Allele1 Threshold Applied | Allele2 Threshold Applied | Total SNPs | Mean A1 Prob | Median A1 Prob | SD A1 Prob | Heterozygosity |
-|:--:|:--:|:--:|:--:|:--:|:--:|:--:|
-| No | Yes | 10024 | 0.9984 | 1 | 0.0096 | 0.4626 |
-| Yes | Yes | 9718 | 0.9998 | 1 | 9.00E-04 | 0.4569 |
-| Minimum \# of SNPs | Yes | 6000 | 1 | 1 | 0 | 0.4495 |
diff --git a/docs/news/index.html b/docs/news/index.html
index 02c6c47..6df9953 100644
--- a/docs/news/index.html
+++ b/docs/news/index.html
@@ -7,7 +7,7 @@
mixder
- 0.9.0
+ 0.10.0
@@ -22,6 +22,7 @@
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
+
Improved error handling in Shiny app (pop-up window will appear if required files or settings are missing)
+
Added option to change number of threads used by EFM; useful if running on a server or grid (default is 0)