Hotfix/update scripts to output stave

scripts/MultiLociBiallelicModel_wrapper/MultiLociBiallelicModel_wrapper.R

@@ -794,8 +794,7 @@ opts = list(
     "--aa_calls",
     type = "character", 
     help = str_c(
-      "TSV containing aminoacid calls per specimen, with the columns: specimen_id,
-      target_id, specimen_id, gene_id, aa_position, ref_aa, and aa."
+      "TSV containing aminoacid calls per specimen, with the columns: specimen_id, gene_id, aa_position, ref_aa, and aa."
     )
   ),
   make_option(
@@ -857,7 +856,7 @@ create_MultiLociBiallelicModel_input <- function(input_path, loci_group) {
   input_data <- read.csv(input_path, na.strings = "NA", sep = "\t")
 
   MLBM_data <- input_data %>%
-    mutate(identifier = paste(target_id, aa_position, sep = ":")) %>%
+    mutate(identifier = paste(gene_id, aa_position, sep = ":")) %>%
     group_by(specimen_id, identifier) %>%
     mutate(value = case_when(
       n_distinct(aa) > 2 ~ NA_real_,                      # More than two distinct amino acids
@@ -885,7 +884,7 @@ create_MultiLociBiallelicModel_input <- function(input_path, loci_group) {
 
   staves_data <- input_data %>%
     mutate(
-      staves = paste(target_id, aa_position, aa, sep = ":"),
+      staves = paste(gene_id, aa_position, aa, sep = ":"),
       state = if_else(
         substr(staves, nchar(staves), nchar(staves)) == ref_aa, 
         0, 
@@ -895,7 +894,7 @@ create_MultiLociBiallelicModel_input <- function(input_path, loci_group) {
     distinct(staves, .keep_all = TRUE) %>%
     select(staves, state) %>%
     mutate(prestave = sub(":([^:]+)$", "", staves))
-  
+
   MLBM_object = list(MLBM_data = MLBM_data, 
                      staves_data = staves_data)
 
@@ -930,8 +929,8 @@ create_MultiLociBiallelicModel_input <- function(input_path, loci_group) {
 
   #Matching groups to create lists of loci
   match_group = input_data %>% 
-    select(target_id, gene_id, aa_position) %>%
-    mutate(identifier = paste(target_id, aa_position, sep = ":")) 
+    select(gene_id, aa_position) %>%
+    mutate(identifier = paste(gene_id, aa_position, sep = ":")) 
 
   merged_data <- match_group %>%
     inner_join(loci_groups, by = c("gene_id", "aa_position"), relationship = "many-to-many")
@@ -948,7 +947,6 @@ create_MultiLociBiallelicModel_input <- function(input_path, loci_group) {
     subtable <- MLBM_data %>% select(all_of(columns))
     return(subtable)
   })
-
   # Name the list elements
   names(result_tables) <- names(result_list)
   MLBM_object$by_group_table <- result_tables
@@ -1055,7 +1053,7 @@ make_stave <- function(variant) {
     group_by(gene) %>%
     summarize(
       positions = str_c(position, collapse = "_"),
-      mutations = str_c(mutation, collapse = ""),
+      mutations = str_c(mutation, collapse = "_"),
       .groups = "drop"
     ) %>%
     mutate(formatted = str_c(gene, ":", positions, ":", mutations)) %>%

scripts/estimate_allele_frequency_naive/estimate_allele_frequency_naive.R

@@ -1,3 +1,8 @@
+# Get the script dir and make path to utils.R
+script_dir <- dirname(normalizePath(sub("--file=", "", commandArgs(trailingOnly = FALSE)[grep("--file=", commandArgs(trailingOnly = FALSE))])))
+utils_path <- file.path(script_dir, "..", "utils", "utils.R")
+
+source(utils_path)
 library(optparse)
 library(rlang)
 
@@ -122,5 +127,5 @@ out <- switch(
   read_count_prop = calculate_af_read_count_prop(aa_calls),
   presence_absence = calculate_af_presence_absence(aa_calls)
 )
-
-readr::write_tsv(out, args$output)
+freq_output <- convert_single_locus_table_to_stave(out, "freq")
+readr::write_tsv(freq_output, args$output)

scripts/slaf_from_stave_mlaf/slaf_from_stave_mlaf.R

@@ -1,3 +1,9 @@
+# Get the script dir and make path to utils.R
+script_dir <- dirname(normalizePath(sub("--file=", "", commandArgs(trailingOnly = FALSE)[grep("--file=", commandArgs(trailingOnly = FALSE))])))
+utils_path <- file.path(script_dir, "..", "utils", "utils.R")
+
+source(utils_path)
+
 library(variantstring)
 library(rlang)
 library(optparse)
@@ -52,10 +58,13 @@ convert_mlaf_to_slaf <- function(dat) {
     dplyr::group_by(.data$group_id, .data$gene, .data$pos, .data$aa) |>
     dplyr::summarize(freq = sum(.data$freq)) |>
     dplyr::ungroup() |>
-    dplyr::arrange("group_id", "gene", "pos", "aa")
+    dplyr::rename(gene_id = gene) |> 
+    dplyr::rename(aa_position = pos) |> 
+    dplyr::arrange("group_id", "gene_id", "aa_position", "aa")
   return(slaf)
 }
 
 mlaf <- load_mlaf(args$mlaf_input)
 slaf <- convert_mlaf_to_slaf(mlaf)
-readr::write_tsv(slaf, args$output)
+slaf_output <- convert_single_locus_table_to_stave(slaf, "freq")
+readr::write_tsv(slaf_output, args$output)

scripts/utils/utils.R

@@ -29,12 +29,11 @@ convert_single_locus_table_to_stave <- function(df, additional_columns=NULL) {
   df %>%
     ungroup %>%
     mutate(variant = paste(gene_id, aa_position, aa, sep = ":")) %>%
-    mutate(ref = paste(gene_id, aa_position, ref_aa, sep = ":")) %>%
     { 
       if (is.null(additional_columns)) {
-        select(., variant, ref)
+        select(., variant)
       } else {
-        select(., variant, ref, all_of(additional_columns))
+        select(., variant, all_of(additional_columns))
       }
     }