Incorporate vibrio characterisation with srst2 into TheiaProk workflows

tasks/species_typing/task_srst2_vibrio.wdl

@@ -0,0 +1,159 @@
+version 1.0
+
+task srst2_vibrio {
+  meta {
+    description: "Computational method for finding spa types in Staphylococcus aureus"
+  }
+  input {
+    File reads1
+    File? reads2
+    String samplename
+    Int srst2_min_cov
+    Int srst2_max_divergence
+    Int srst2_min_depth
+    Int srst2_min_edge_depth
+    Int srst2_gene_max_mismatch
+    String docker = "quay.io/kapsakcj/srst2:0.2.0-vcholerae" # TODO: Update with container including vibrio db
+    Int disk_size = 100
+    Int cpu = 4
+  }
+  command <<<
+    if [ -z "~{reads2}" ] ; then
+      INPUT_READS="--input_se ~{reads1}"
+    else
+      # This task expects/requires that input FASTQ files end in "_1.clean.fastq.gz" and "_2.clean.fastq.gz"
+      # which is the syntax from TheiaProk read cleaning tasks
+      INPUT_READS="--input_pe ~{reads1} ~{reads2} --forward _1.clean --reverse _2.clean"
+    fi
+
+    srst2 --version 2>&1 | tee VERSION
+    srst2 \
+      ${INPUT_READS} \
+      --gene_db /vibrio-cholerae-db/vibrio_230224.fasta \
+      --output ~{samplename} \
+      --min_coverage ~{srst2_min_cov} \
+      --max_divergence ~{srst2_max_divergence} \
+      --min_depth ~{srst2_min_depth} \
+      --min_edge_depth ~{srst2_min_edge_depth} \
+      --gene_max_mismatch ~{srst2_gene_max_mismatch}
+
+    # capture output TSV
+    mv ~{samplename}__genes__*__results.txt ~{samplename}.tsv
+
+    # capture detailed output TSV - not available if no results are outputed
+    mv ~{samplename}__fullgenes__*__results.txt ~{samplename}.detailed.tsv || echo "No results" >  ~{samplename}.detailed.tsv
+
+    # parsing block to account for when output columns do not exist
+    python <<CODE
+    import csv
+    import re
+
+    # Converting TSV file into list of dictionaries
+    def tsv_to_dict(filename):
+      result_list=[]
+      with open(filename) as file_obj:
+          reader = csv.DictReader(file_obj, delimiter='\t')
+          for row in reader:
+              result_list.append(dict(row))
+      # only one sample is run, so there's only one row, flattening list
+      return result_list[0]
+
+    # Converting None to empty string
+    conv = lambda i : i or '-'
+
+    # Make characters human-readable 
+    def translate_chars(string):
+      translation = []
+      if '?' in string:
+        translation.append("low depth/uncertain")
+      if '-' in string:
+        translation.append("not detected")
+
+      # in case we want to retrieve the allele information
+      string = re.sub("\*|\?|-", "", string)
+
+      if len(translation) > 0:
+        return '(' + ';'.join(translation) + ')'
+      return ""
+
+    # load output TSV as dict 
+    row = tsv_to_dict('~{samplename}.tsv')
+
+    # presence or absence genes - ctxA, ompW and toxR
+    with open("ctxA", "wb") as ctxA_fh:
+      value = row.get("ctxA")
+      presence = translate_chars(conv(value))
+      if presence == "(not detected)":
+        ctxA_fh.write(presence)
+      else:
+        result = "present" + ' ' + presence
+        ctxA_fh.write(result.strip())
+
+    with open("ompW", "wb") as ompW_fh:
+      value = row.get("ompW")
+      presence = translate_chars(conv(value))
+      if presence == "(not detected)":
+        ompW_fh.write(presence)
+      else:
+        result = "present" + ' ' + presence
+        ompW_fh.write(result.strip())
+
+    with open("toxR", "wb") as toxR_fh:
+      value = row.get("toxR")
+      presence = translate_chars(conv(value))
+      if presence == "(not detected)":
+        toxR_fh.write(presence)
+      else:
+        result = "present" + ' ' + presence
+        toxR_fh.write(result.strip())
+
+    # biotype - tcpA classical or tcpA ElTor
+    with open("BIOTYPE", "wb") as biotype_fh:
+      value_ElTor = translate_chars(conv(row.get("tcpA_ElTor")))
+      value_classical = translate_chars(conv(row.get("tcpA_classical")))
+
+      if value_ElTor == "(not detected)" and value_classical == "(not detected)":
+        biotype_fh.write("(not detected)")
+      else:
+        if value_ElTor == "(not detected)":
+          result = "tcpA_Classical" + ' ' + value_classical
+          biotype_fh.write(result.strip())
+        else:
+          result = "tcpA_ElTor" + ' ' + value_ElTor
+          biotype_fh.write(result.strip())
+
+    # serotype - O1 or O139
+    with open("SEROTYPE", "wb") as serotype_fh:
+      value_O1 = translate_chars(conv(row.get("wbeN_O1")))
+      value_O139 = translate_chars(conv(row.get("wbfR_O139")))
+
+      if value_O1 == "(not detected)" and value_O139 == "(not detected)":
+        serotype_fh.write("(not detected)")
+      else:
+        if value_O1 == "(not detected)":
+          result = "O139" + ' ' + value_O139
+          serotype_fh.write(result.strip())
+        else:
+          result = "O1" + ' ' + value_O1
+          serotype_fh.write(result.strip())
+    CODE
+  >>>
+  output {
+      File srst2_detailed_tsv = "~{samplename}.detailed.tsv"
+      String srst2_version = read_string("VERSION")
+      String srst2_vibrio_ctxA = read_string("ctxA")
+      String srst2_vibrio_ompW = read_string("ompW")
+      String srst2_vibrio_toxR = read_string("toxR")
+      String srst2_vibrio_biotype = read_string("BIOTYPE")
+      String srst2_vibrio_serotype = read_string("SEROTYPE")
+  }
+  runtime {
+    docker: "~{docker}"
+    memory: "8 GB"
+    cpu: cpu
+    disks: "local-disk " + disk_size + " SSD"
+    disk: disk_size + " GB"
+    maxRetries: 3
+    preemptible: 0
+  }
+}

tasks/taxon_id/task_gambit.wdl

@@ -175,6 +175,8 @@ task gambit {
       if [[ ${predicted_taxon} == *"Streptococcus pneumoniae"* ]]; then 
         merlin_tag="Streptococcus pneumoniae"
       fi
+    elif [[ ${predicted_taxon} == *"Vibrio"* ]]; then 
+      merlin_tag="Vibrio"
     else 
       merlin_tag="None"
     fi

tasks/utilities/task_broad_terra_tools.wdl

@@ -292,6 +292,13 @@ task export_taxon_tables {
     String? agrvate_agr_num_frameshifts
     String? agrvate_version
     String? agrvate_docker
+    File? srst2_vibrio_detailed_tsv
+    String? srst2_vibrio_version
+    String? srst2_vibrio_ctxA
+    String? srst2_vibrio_ompW
+    String? srst2_vibrio_toxR
+    String? srst2_vibrio_serotype
+    String? srst2_vibrio_biotype
   }
   command <<<
 
@@ -511,7 +518,7 @@ task export_taxon_tables {
       "tbprofiler_dr_type": "~{tbprofiler_dr_type}",
       "tbprofiler_resistance_genes": "~{tbprofiler_resistance_genes}",
       "tbprofiler_additional_outputs_csv": "~{tbprofiler_additional_outputs_csv}",
-      "tbprofiler_laboratorian_report_csv": "~{tbprofiler_laboratorian_report_csv}"
+      "tbprofiler_laboratorian_report_csv": "~{tbprofiler_laboratorian_report_csv}",
       "tbprofiler_gene_name": "~{tbprofiler_gene_name}",
       "tbprofiler_locus_tag": "~{tbprofiler_locus_tag}",
       "tbprofiler_variant_substitutions": "~{tbprofiler_variant_substitutions}",
@@ -609,7 +616,14 @@ task export_taxon_tables {
       "agrvate_agr_multiple": "~{agrvate_agr_multiple}",
       "agrvate_agr_num_frameshifts": "~{agrvate_agr_num_frameshifts}",
       "agrvate_version": "~{agrvate_version}",
-      "agrvate_docker": "~{agrvate_docker}"
+      "agrvate_docker": "~{agrvate_docker}",
+      "srst2_vibrio_detailed_tsv": "~{srst2_vibrio_detailed_tsv}",
+      "srst2_vibrio_version": "~{srst2_vibrio_version}",
+      "srst2_vibrio_ctxA": "~{srst2_vibrio_ctxA}",
+      "srst2_vibrio_ompW": "~{srst2_vibrio_ompW}",
+      "srst2_vibrio_toxR": "~{srst2_vibrio_toxR}",
+      "srst2_vibrio_serotype": "~{srst2_vibrio_serotype}",
+      "srst2_vibrio_biotype": "~{srst2_vibrio_biotype}",
     }
 
     with open("~{samplename}_terra_table.tsv", "w") as outfile:

tests/workflows/test_wf_theiaprok_illumina_pe.yml

@@ -107,7 +107,8 @@
     - path: miniwdl_run/call-clean_check_reads/work/_miniwdl_inputs/0/test_1.clean.fastq.gz
     - path: miniwdl_run/call-clean_check_reads/work/_miniwdl_inputs/0/test_2.clean.fastq.gz
     - path: miniwdl_run/call-gambit/command
-      md5sum: ffda45de2bad7a2206f507bf4485c930
+      md5sum: 7998373971cbc16b4654736280d8d4bd
+      md5sum: 1931093118b504186aa3776dbca5adb4
     - path: miniwdl_run/call-gambit/inputs.json
       contains: ["assembly", "fasta", "samplename", "test"]
     - path: miniwdl_run/call-gambit/outputs.json

workflows/wf_merlin_magic.wdl

@@ -23,6 +23,7 @@ import "../tasks/species_typing/task_pbptyper.wdl" as pbptyper
 import "../tasks/species_typing/task_poppunk_streppneumo.wdl" as poppunk_spneumo
 import "../tasks/species_typing/task_pasty.wdl" as pasty
 import "../tasks/gene_typing/task_abricate.wdl" as abricate_task
+import "../tasks/species_typing/task_srst2_vibrio.wdl" as srst2_vibrio_task
 
 workflow merlin_magic {
   meta {
@@ -46,6 +47,11 @@ workflow merlin_magic {
     Boolean call_shigeifinder_reads_input = false
     Boolean tbprofiler_additional_outputs = false
     String output_seq_method_type = "WGS"
+    Int srst2_min_cov = 80
+    Int srst2_max_divergence = 20
+    Int srst2_min_depth = 5
+    Int srst2_min_edge_depth = 2
+    Int srst2_gene_max_mismatch = 200
   }
     if (merlin_tag == "Acinetobacter baumannii") {
     call kaptive.kaptive {
@@ -226,6 +232,19 @@ workflow merlin_magic {
       }  
     }
   }
+  if (merlin_tag == "Vibrio") {
+    call srst2_vibrio_task.srst2_vibrio {
+      input:
+        reads1 = read1,
+        reads2 = read2,
+        samplename = samplename,
+        srst2_min_cov = srst2_min_cov,
+        srst2_max_divergence = srst2_max_divergence,
+        srst2_min_depth = srst2_min_depth,
+        srst2_min_edge_depth = srst2_min_edge_depth,
+        srst2_gene_max_mismatch = srst2_gene_max_mismatch
+    }
+  }
 
   output {
   # Ecoli Typing
@@ -419,5 +438,13 @@ workflow merlin_magic {
   String? seroba_ariba_serotype = seroba_task.seroba_ariba_serotype
   String? seroba_ariba_identity = seroba_task.seroba_ariba_identity
   File? seroba_details = seroba_task.seroba_details
- }
-}
+  # Vibrio
+  File? srst2_vibrio_detailed_tsv = srst2_vibrio.srst2_detailed_tsv
+  String? srst2_vibrio_version = srst2_vibrio.srst2_version
+  String? srst2_vibrio_ctxA = srst2_vibrio.srst2_vibrio_ctxA
+  String? srst2_vibrio_ompW = srst2_vibrio.srst2_vibrio_ompW
+  String? srst2_vibrio_toxR = srst2_vibrio.srst2_vibrio_toxR
+  String? srst2_vibrio_serotype = srst2_vibrio.srst2_vibrio_serotype
+  String? srst2_vibrio_biotype = srst2_vibrio.srst2_vibrio_biotype
+  }
+ }

workflows/wf_theiaprok_illumina_pe.wdl

@@ -484,7 +484,14 @@ workflow theiaprok_illumina_pe {
             pasty_docker = merlin_magic.pasty_docker,
             pasty_comment = merlin_magic.pasty_comment,
             qc_check = qc_check_task.qc_check,
-            qc_standard = qc_check_task.qc_standard
+            qc_standard = qc_check_task.qc_standard,
+            srst2_vibrio_detailed_tsv = merlin_magic.srst2_vibrio_detailed_tsv,
+            srst2_vibrio_version = merlin_magic.srst2_vibrio_version,
+            srst2_vibrio_ctxA = merlin_magic.srst2_vibrio_ctxA,
+            srst2_vibrio_ompW = merlin_magic.srst2_vibrio_ompW,
+            srst2_vibrio_toxR = merlin_magic.srst2_vibrio_toxR,
+            srst2_vibrio_serotype = merlin_magic.srst2_vibrio_serotype,
+            srst2_vibrio_biotype = merlin_magic.srst2_vibrio_biotype
         }
       }
     }
@@ -792,5 +799,13 @@ workflow theiaprok_illumina_pe {
     String? seroba_ariba_serotype = merlin_magic.seroba_ariba_serotype
     String? seroba_ariba_identity = merlin_magic.seroba_ariba_identity
     File? seroba_details = merlin_magic.seroba_details
+    # Vibrio Typing
+    File? srst2_vibrio_detailed_tsv = merlin_magic.srst2_vibrio_detailed_tsv
+    String? srst2_vibrio_version = merlin_magic.srst2_vibrio_version
+    String? srst2_vibrio_ctxA = merlin_magic.srst2_vibrio_ctxA
+    String? srst2_vibrio_ompW = merlin_magic.srst2_vibrio_ompW
+    String? srst2_vibrio_toxR = merlin_magic.srst2_vibrio_toxR
+    String? srst2_vibrio_biotype = merlin_magic.srst2_vibrio_biotype
+    String? srst2_vibrio_serotype = merlin_magic.srst2_vibrio_serotype
   }
 }

workflows/wf_theiaprok_illumina_se.wdl

@@ -424,7 +424,14 @@ workflow theiaprok_illumina_se {
             agrvate_agr_multiple = merlin_magic.agrvate_agr_multiple,
             agrvate_agr_num_frameshifts = merlin_magic.agrvate_agr_num_frameshifts,
             agrvate_version = merlin_magic.agrvate_version,
-            agrvate_docker = merlin_magic.agrvate_docker
+            agrvate_docker = merlin_magic.agrvate_docker,
+            srst2_vibrio_detailed_tsv = merlin_magic.srst2_vibrio_detailed_tsv,
+            srst2_vibrio_version = merlin_magic.srst2_vibrio_version,
+            srst2_vibrio_ctxA = merlin_magic.srst2_vibrio_ctxA,
+            srst2_vibrio_ompW = merlin_magic.srst2_vibrio_ompW,
+            srst2_vibrio_toxR = merlin_magic.srst2_vibrio_toxR,
+            srst2_vibrio_serotype = merlin_magic.srst2_vibrio_serotype,
+            srst2_vibrio_biotype = merlin_magic.srst2_vibrio_biotype
         }
       }
     }
@@ -712,5 +719,13 @@ workflow theiaprok_illumina_se {
     String? poppunk_GPS_db_version = merlin_magic.poppunk_GPS_db_version
     String? poppunk_version = merlin_magic.poppunk_version
     String? poppunk_docker = merlin_magic.poppunk_docker
+    # Vibrio Typing
+    File? srst2_vibrio_detailed_tsv = merlin_magic.srst2_vibrio_detailed_tsv
+    String? srst2_vibrio_version = merlin_magic.srst2_vibrio_version
+    String? srst2_vibrio_ctxA = merlin_magic.srst2_vibrio_ctxA
+    String? srst2_vibrio_ompW = merlin_magic.srst2_vibrio_ompW
+    String? srst2_vibrio_toxR = merlin_magic.srst2_vibrio_toxR
+    String? srst2_vibrio_biotype = merlin_magic.srst2_vibrio_biotype
+    String? srst2_vibrio_serotype = merlin_magic.srst2_vibrio_serotype
   }
 }