Condaify all the tools

src/bfx_tools/hlarp.ml

@@ -7,9 +7,8 @@ type hla_result = [
   ]
 
 let run ~(run_with:Machine.t)
-    ?(edges=[])
     ~(hla_result:hla_result)
-    ~output_path ~extract_alleles ()
+    ~output_path
   =
   let open KEDSL in
   let open Ketrew_pure.Target.Volume in
@@ -25,17 +24,10 @@ let run ~(run_with:Machine.t)
       ~requirements:[`Self_identification ["hlarp"]]
       Program.(
         Machine.Tool.init hlarp
-        && shf "hlarp %s %s > %s" subcommand hla_result_directory output_path
-        && sh
-          (if extract_alleles
-           then sprintf
-               "awk -F , '{ gsub(/^[ \t]+|[ \t]+$/,\"\", $2); print $2}' %s \
-                | tail -n +2 \
-                | sed \"s/'//\" > %s.tmp \
-                && mv %s.tmp %s"
-               output_path output_path output_path output_path
-           else "")) in
-  let edges = hla_result_dep :: edges @ [
+        && shf "hlarp %s %s > %s" subcommand hla_result_directory output_path)
+  in
+  let edges = [
+      hla_result_dep;
       depends_on (Machine.Tool.ensure hlarp);
       on_failure_activate
         (Workflow_utilities.Remove.file ~run_with output_path);

src/bfx_tools/isovar.ml

@@ -30,9 +30,7 @@ let run ~(run_with: Machine.t)
   ~configuration
   ~reference_build ~vcf ~bam ~output_file =
   let open KEDSL in
-  let isovar_tool =
-    Machine.get_tool run_with Machine.Tool.Definition.(create "isovar")
-  in
+  let isovar_tool = Machine.get_tool run_with Machine.Tool.Default.isovar in
   let genome = Machine.(get_reference_genome run_with reference_build)
     |> Reference_genome.name
   in

src/bfx_tools/optitype.ml

@@ -1,4 +1,3 @@
-
 open Biokepi_run_environment
 open Common
 
@@ -24,6 +23,24 @@ let move_optitype_product ?host ~path o =
     method path = path
   end
 
+let get_optitype_data_folder =
+  let tname, tversion =
+    let tool_def = Machine.Tool.Default.optitype in
+    Machine.Tool.Definition.(get_name tool_def, get_version tool_def)
+  in
+  sprintf "${CONDA_PREFIX}/share/%s-%s/"
+    tname (match tversion with None -> "*" | Some v -> v)
+
+(* copy sample config file and the required data over;
+   then, adjust the razers3 path in the config *)
+let prepare_optidata = 
+  let open KEDSL.Program in
+  let optidata_path = get_optitype_data_folder in
+  shf "cp -r %s/data ." optidata_path && (* HLA reference data *)
+  shf "cp -r %s/config.ini.example config.ini" optidata_path && 
+  sh "sed -i.bak \"s|\\/path\\/to\\/razers3|$(which razers3)|g\" config.ini"
+
+
 (**
    Run OptiType in [`RNA] or [`DNA] mode.
 
@@ -43,17 +60,13 @@ let hla_type ~work_dir ~run_with ~fastq ~run_name nt
         Machine.Tool.init tool
         && exec ["mkdir"; "-p"; work_dir]
         && exec ["cd"; work_dir]
-        && sh "cp -r ${OPTITYPE_DATA}/data ." (* HLA reference data *)
-        && (* config example *)
-        sh "cp -r ${OPTITYPE_DATA}/config.ini.example config.ini"
-        && (* adjust config razers3 path *)
-        sh "sed -i.bak \"s|\\/path\\/to\\/razers3|$(which razers3)|g\" config.ini"
-        &&
-        shf "OptiTypePipeline --verbose --input %s %s %s -o %s "
-          (Filename.quote r1_path)
-          (Option.value_map ~default:"" r2_path_opt ~f:Filename.quote)
-          (match nt with | `DNA -> "--dna" | `RNA -> "--rna")
-          run_name)
+        && prepare_optidata
+        && shf "OptiTypePipeline.py --verbose -c ./config.ini \
+                --input %s %s %s -o %s"
+            (Filename.quote r1_path)
+            (Option.value_map ~default:"" r2_path_opt ~f:Filename.quote)
+            (match nt with | `DNA -> "--dna" | `RNA -> "--rna")
+            run_name)
   in
   let product =
     let host = Machine.as_host run_with in
@@ -77,3 +90,121 @@ let hla_type ~work_dir ~run_with ~fastq ~run_name nt
           (Workflow_utilities.Remove.directory ~run_with work_dir);
       ]
     )
+
+(*
+  Optitype depends on alignment of reads onto the HLA-locus
+  as a preliminary filtering step, but the default aligner, razers3,
+  is so memory hungry that, the run fails on a ~50G memory machine
+  when the number of reads per FASTQ is approximately more than 200M.
+
+  The following variation of optitype run makes use of `bwa mem` instead
+  of `razers3` to do the initial filtering and some experimentation
+  with real patient data proved that this doesn't affect the results
+  despite their worrisome warning on the site.
+
+  This is only tested on the DNA arm of the pipeline, so is restricted
+  to that use case.
+*)
+let dna_hla_type_with_bwamem
+    ?(configuration = Bwa.Configuration.Mem.default)
+    ~work_dir ~run_with ~fastq ~run_name
+  =
+  let open KEDSL in
+  (* We need to pull in bwa mem and samtools to get some help *)
+  let optitype = Machine.get_tool run_with Machine.Tool.Default.optitype in
+  let bwa = Machine.get_tool run_with Machine.Tool.Default.bwa in
+  let samtools = Machine.get_tool run_with Machine.Tool.Default.samtools in
+  let bwa_wd = work_dir // "bwamem" in
+  let dna_hla_ref_path = 
+      get_optitype_data_folder // "data/hla_reference_dna.fasta"
+  in
+  (* Step 1: prepare hla reference indexes *)
+  let index_hla_wf =
+    let name = "Index OptiType's DNA-based HLA reference with bwa" in
+    let edges = [
+        depends_on (Machine.Tool.ensure bwa);
+        depends_on (Machine.Tool.ensure optitype);
+        on_failure_activate
+          (Workflow_utilities.Remove.directory ~run_with bwa_wd);]
+    in
+    let make = 
+      Machine.run_big_program run_with ~name
+        ~self_ids:["optitype"; "hla"; "dna"; "bwa index"]
+        Program.(
+          Machine.Tool.init optitype
+          && Machine.Tool.init bwa
+          && Machine.Tool.init optitype
+          && shf "mkdir -p %s" bwa_wd
+          && shf "bwa index %s" dna_hla_ref_path
+        )
+    in
+    let product = 
+      Workflow_utilities.Variable_tool_paths.single_file 
+        ~run_with ~tool:optitype (dna_hla_ref_path ^ ".bwt")
+    in
+    workflow_node product ~name ~make ~edges
+  in
+  (* Step 2: Map the input fastq to this pseudo reference genome via
+     `bwa mem` and turn the alignment back into fastq while keeping
+     only the reads that mapped to reduce OptiType's future memory
+     consumption.
+  *)
+  let filter_hla_reads_wf =
+    let name = sprintf "Filter out non-HLA-mapping reads: %s" run_name in
+    let edges = [
+        depends_on (Machine.Tool.ensure bwa);
+        depends_on (Machine.Tool.ensure samtools);
+        depends_on index_hla_wf; (* No indexing, no mapping *)
+        depends_on fastq;
+        on_failure_activate
+          (Workflow_utilities.Remove.directory ~run_with bwa_wd);]
+    in
+    let bwa2sam2fastq fqpath = (* the whole pipeline *)
+      let outfq_path = 
+        let fqbase = fqpath |> Filename.basename |> Filename.chop_extension in
+        bwa_wd // (sprintf "%s-hla_mapping.fastq" fqbase)
+      in
+      let processors = Machine.max_processors run_with in
+      let bwamem_part = Bwa.(
+        sprintf "bwa mem -t %d -O %d -E %d -B %d %s %s"
+          processors
+          configuration.Configuration.Mem.gap_open_penalty
+          configuration.Configuration.Mem.gap_extension_penalty
+          configuration.Configuration.Mem.mismatch_penalty
+          dna_hla_ref_path
+          (Filename.quote fqpath))
+      in
+      let samtools_part =
+        (* -F4 filters *out* all reads that do not map. 
+           See SAM/BAM flags for more information:
+              $ samtools flags
+        *)
+        sprintf "samtools fastq -F4 - -0 %s" outfq_path
+      in
+      String.concat ~sep:" | " [bwamem_part; samtools_part],
+      outfq_path
+    in
+    let in_r1, in_r2_opt = fastq#product#paths in
+    let filter_r1, out_r1 = bwa2sam2fastq in_r1 in
+    let filter_r2, out_r2 = 
+      match in_r2_opt with
+      | None -> "echo 'Second pair is missing'", None
+      | Some r2p -> let (f, o) = bwa2sam2fastq r2p in (f, Some o)
+    in
+    let make =
+      Machine.run_big_program run_with ~name
+        ~self_ids:["optitype"; "hla"; "dna"; "filtering"]
+        Program.(
+          Machine.Tool.init optitype
+          && Machine.Tool.init bwa
+          && Machine.Tool.init samtools
+          && shf "mkdir -p %s" bwa_wd
+          && sh filter_r1 && sh filter_r2
+        )
+    in
+    let product = transform_fastq_reads fastq#product out_r1 out_r2 in
+    workflow_node product ~name ~make ~edges
+  in
+  (* Step 3: Run OptiType as usual on the new filtered down FASTQ(s) *)
+  let owd = work_dir // "optitype" in
+  hla_type ~work_dir:owd ~run_with ~fastq:filter_hla_reads_wf ~run_name `DNA