Fix missing reference record in afwdist results

Renames outgroup params for distance tree rules

Author: ahmig

workflow/rules/distances.smk

@@ -22,7 +22,7 @@ rule extract_afwdist_variants:
 rule afwdist_weighted_distances:
     conda: "../envs/afwdist.yaml"
     params:
-        extra_args = ""
+        extra_args = "",
     input:
         variants = OUTDIR/f"{OUTPUT_NAME}.variants.afwdist.csv",
         reference = OUTDIR/f"{OUTPUT_NAME}.ancestor.fasta",
@@ -41,7 +41,7 @@ rule afwdist_weighted_distances:
 rule format_afwdist_results:
     conda: "../envs/biopython.yaml"
     params:
-        samples = sorted(iter_samples()),
+        samples = sorted(iter_samples()) + [config["ALIGNMENT_REFERENCE"]],
     input:
         distances = OUTDIR/f"{OUTPUT_NAME}.distances.raw.csv",
     output:

workflow/rules/report.smk

@@ -260,7 +260,7 @@ rule allele_freq_tree_data:
     conda: "../envs/renv.yaml"
     params:
         use_bionj = config["USE_BIONJ"],
-        ref_name = config["ALIGNMENT_REFERENCE"],
+        outgroup_id = config["ALIGNMENT_REFERENCE"],
     input:
         dist = OUTDIR/f"{OUTPUT_NAME}.distances.csv",
     output:
@@ -275,7 +275,7 @@ rule allele_freq_tree_plot:
     conda: "../envs/renv.yaml"
     params:
         design = config["PLOTS"],
-        ref_name = config["ALIGNMENT_REFERENCE"],
+        outgroup_id = config["ALIGNMENT_REFERENCE"],
         plot_height_mm = 119.4,
         plot_width_mm = 159.2,
     input:
@@ -293,7 +293,7 @@ rule allele_freq_tree_plot:
 rule time_signal_data:
     conda: "../envs/renv.yaml"
     params:
-        ref_name = config["ALIGNMENT_REFERENCE"],
+        outgroup_id = config["ALIGNMENT_REFERENCE"],
     input:
         tree = report(REPORT_DIR_TABLES/"allele_freq_tree.nwk"),
         metadata = config["METADATA"],

workflow/scripts/format_afwdist_results.py

@@ -13,6 +13,7 @@
 
     logging.info("Initializing formatted output")
     output = pd.DataFrame(
+        data=0.0,
         columns=snakemake.params.samples,
         index=snakemake.params.samples,
         dtype="float64"

workflow/scripts/report/allele_freq_tree_data.R

@@ -46,7 +46,7 @@ tree <- matrix %>%
   as.matrix() %>%
   as.dist() %>%
   tree_method() %>%
-  root(snakemake@params$ref_name, resolve.root = TRUE)
+  root(snakemake@params$outgroup_id, resolve.root = TRUE)
 
 # Resolve possible negative edge lengths in tree
 log_info("Resolving negative edge lengths")

workflow/scripts/report/allele_freq_tree_plot.R

@@ -25,16 +25,16 @@ metadata <- read_csv(snakemake@input$metadata)
 
 # Extract names from records
 log_debug("Reading FASTA")
-study_names <- read.dna(
+study_records <- read.dna(
   snakemake@input$study_fasta,
   format = "fasta",
   as.matrix = FALSE,
 )
 log_debug("Extracting names")
-study_names <- study_names[
-  !startsWith(names(study_names), snakemake@params$ref_name)
+study_records <- study_records[
+  !startsWith(names(study_records), snakemake@params$outgroup_id)
 ]
-study_names <- names(study_names)
+study_names <- names(study_records)
 
 # Obtain sample names ordered by CollectionDate
 log_debug("Sorting names by collection date")
@@ -59,19 +59,21 @@ tree_tiplab <- data.frame(
   ) %>%
   ungroup() %>%
   add_row(
-    ID = snakemake@params$ref_name,
+    ID = snakemake@params$outgroup_id,
     order = 0,
-    tip_label = snakemake@params$ref_name
+    tip_label = snakemake@params$outgroup_id
   )
 
 # Read distance tree and root
 log_info("Reading tree")
 tree <- read.tree(snakemake@input$tree) %>%
-  root(snakemake@params$ref_name, resolve.root = TRUE)
+  root(snakemake@params$outgroup_id, resolve.root = TRUE)
 
 # Plot
 log_info("Plotting distance tree")
-max.tip.length <- max(node.depth.edgelength(tree)[1:length(tree$tip.label)])
+max.tip.length <- max(
+  node.depth.edgelength(tree)[seq_along(length(tree$tip.label))]
+)
 p <- ggtree(tree) %<+% tree_tiplab +
   geom_tiplab(aes(label = tip_label)) +
   geom_treescale(1.1 * max.tip.length) +

workflow/scripts/report/time_signal_data.R

@@ -18,7 +18,7 @@ log_threshold(INFO)
 # Read distance tree and root
 log_info("Reading tree")
 tree <- read.tree(snakemake@input$tree) %>%
-  root(snakemake@params$ref_name, resolve.root = TRUE)
+  root(snakemake@params$outgroup_id, resolve.root = TRUE)
 
 # Read metadata
 log_info("Reading metadata")
@@ -33,14 +33,15 @@ time.signal <- distRoot(
 ) %>%
   as.data.frame() %>%
   rownames_to_column(var = "ID") %>%
-  filter(ID != snakemake@params$ref_name) %>%
+  filter(ID != snakemake@params$outgroup_id) %>%
   rename(distance = ".") %>%
   left_join(
     select(
       metadata,
       ID,
       CollectionDate
-    )
+    ),
+    by = "ID"
   ) %>%
   mutate(
     date_interval = as.numeric(