change var count annotation to call annotation

polio-nanopore · Oct 24, 2023 · 826c743 · 826c743
1 parent d8eb66d
commit 826c743
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 8 deletions.
diff --git a/piranha/analysis/phylo_functions.py b/piranha/analysis/phylo_functions.py
@@ -53,10 +53,19 @@ def get_seqs_and_clusters(sample_seqs,supplementary_sequences,reference_sequence
                 seq_metadata[name][KEY_REFERENCE_GROUP] = ref_group
 
                 var_count = description_dict[KEY_VARIANT_COUNT]
-                if var_count == "NA":
-                    seq_metadata[name][KEY_VARIANT_COUNT] = 0
-                else:
-                    seq_metadata[name][KEY_VARIANT_COUNT] = int(var_count)
+
+                call = ref_group
+                if ref_group.startswith("Sabin"):
+                    # configured number of mutations in sabin for the call threshold of VDPV
+                    call_threshold = CALL_THRESHOLD_DICT[ref_group]
+                    if int(var_count) > call_threshold:
+                        call = "VDPV"
+                    elif var_count == 0:
+                        call = "Sabin"
+                    else:
+                        call = "Sabin-like"
+
+                seq_metadata[name][KEY_CALL] = call
 
 
     print(green("Reference groups for phylo pipeline:"))
@@ -73,7 +82,7 @@ def get_seqs_and_clusters(sample_seqs,supplementary_sequences,reference_sequence
                     seq_metadata[record.id][KEY_SAMPLE] = record.id
                     seq_metadata[record.id][KEY_SOURCE] = "Background"
                     seq_metadata[record.id][KEY_REFERENCE_GROUP] = ref_group
-                    seq_metadata[record.id][KEY_VARIANT_COUNT] = 0
+                    seq_metadata[record.id][KEY_CALL] = ref_group
 
     if supplementary_metadata:
         with open(supplementary_metadata, "r") as f:
@@ -98,7 +107,7 @@ def get_seqs_and_clusters(sample_seqs,supplementary_sequences,reference_sequence
                 seq_metadata[record.id][KEY_SAMPLE] = record.id
 
                 seq_metadata[record.id][KEY_REFERENCE_GROUP] = ref_group
-                seq_metadata[record.id][KEY_VARIANT_COUNT] = 0
+                seq_metadata[record.id][KEY_CALL] = ref_group
 
                 if "Sabin" in record.description:
                     seq_metadata[record.id][KEY_SOURCE] = "Sabin"

diff --git a/piranha/utils/config.py b/piranha/utils/config.py
@@ -39,6 +39,7 @@
 KEY_EPID="EPID"
 KEY_VARIANTS = "variants"
 KEY_VARIANT_COUNT = "variant_count"
+KEY_CALL = "call"
 KEY_WELL = "well"
 KEY_ALL_METADATA="all_metadata_to_header"
 KEY_REFERENCE = "reference"
@@ -65,7 +66,7 @@
 KEY_PHYLO_METADATA_COLUMNS = "phylo_metadata_columns"
 KEY_SUPPLEMENTARY_METADATA_COLUMNS = "supplementary_metadata_columns"
 
-VALUE_PHYLO_METADATA_COLUMNS = [KEY_VARIANT_COUNT,KEY_DATE,KEY_EPID]
+VALUE_PHYLO_METADATA_COLUMNS = [KEY_CALL,KEY_DATE,KEY_EPID]
 VALUE_SUPPLEMENTARY_METADATA_COLUMNS = ["location","lineage"]
 VALUE_SUPPLEMENTARY_METADATA_ID_COLUMN = "sequence_name"
 
@@ -123,7 +124,7 @@
 VALUE_OUTPUT_PREFIX = "analysis"
 VALUE_SUMMARY_HEADERS = ["taxon","sites","haplotype","num_reads","make_cns"]
 VALUE_REFERENCES_FOR_CNS = ["Sabin1-related","Sabin2-related","Sabin3-related","WPV1"]
-VALUE_PHYLO_HEADER = ["name","sample","barcode","source","reference_group","variant_count"]
+VALUE_PHYLO_HEADER = ["name","sample","barcode","source","reference_group","call"]
 
 VALUE_SAMPLE_TYPE = "stool"
 VALUE_ANALYSIS_MODE = "vp1"