more fixes for rer converge, still not working though

app/wip.ml

@@ -7,9 +7,11 @@ let target =
     ()
 
 let () =
-  let loggers = [ Bistro_utils.Console_logger.create () ] in
-  let allowed_containers = [ `Docker ] in
-  Bistro.Workflow.path target
-  |> Bistro_engine.Scheduler.simple_eval_exn ~allowed_containers ~np:4
-       ~mem:(`GB 4) ~loggers
-  |> print_endline
+  try
+    let loggers = [ Bistro_utils.Console_logger.create () ] in
+    let allowed_containers = [ `Docker ] in
+    Bistro.Workflow.path target
+    |> Bistro_engine.Scheduler.simple_eval_exn ~allowed_containers ~np:4
+      ~mem:(`GB 4) ~loggers
+    |> print_endline
+  with Failure msg -> print_endline msg

lib/orthomam.ml

@@ -55,12 +55,27 @@ let remove_unobserved_sequences_from_alignment phylip : phylip file =
   in
   Workflow.path_plugin ~descr:"orthomam.remove_unobserved_sequences_from_alignment" f
 
-let clip_tree_on_alignment (tree:(#newick as 'a) file) (ali: phylip file) : 'a file =
+let clip_tree_on_alignment ?(handle_tags = true) (tree:(#newick as 'a) file) (ali: phylip file) : 'a file =
   let f = fun%workflow dest ->
     let open Phylogenetics in
+    let handle_tags = [%param handle_tags] in
     let tree = Newick.from_file [%path tree] in
     let ali = Phylip.read_exn ~strict:false [%path ali] in
     let ali_species = List.map ali.items ~f:(fun it -> it.name) in
+    let remove_nodes_with_single_child =
+      if handle_tags then
+        Convergence_tree.remove_nodes_with_single_child
+      else
+        Tree.simplify_node_with_single_child
+          ~merge_branch_data:(fun (branches : Newick.branch_info list) ->
+              let length =
+                List.fold branches ~init:0. ~f:(fun acc bi ->
+                    acc +. Option.value_exn bi.length
+                  )
+              in
+              { Newick.tags = [] ; length = Some length }
+            )
+    in
     let clipped_tree =
       Newick.map_inner_tree tree ~f:(fun tree ->
           match
@@ -68,7 +83,8 @@ let clip_tree_on_alignment (tree:(#newick as 'a) file) (ali: phylip file) : 'a f
               (fun bi -> bi.Newick.name) ali_species
           with
           | None -> failwith "Tree has no leaf in alignment"
-          | Some filtered_tree -> filtered_tree
+          | Some filtered_tree ->
+            remove_nodes_with_single_child filtered_tree
         )
     in
     Newick.to_file clipped_tree dest
@@ -596,16 +612,15 @@ let draw_site q pos =
   in
   Workflow.path_plugin ~descr:"orthomam." f
 
-
 let rer_converge ?max_read ?transform ?weighted ?scale ?continuous ~db
     ~phenotypes () =
   let open Bistro_bio in
   let module RER = Rer_converge in
-  let species_tree = estimated_codon_tree db in
+  let species_tree = estimated_nucleotide_tree db in
   let queries =
     alignments_of_db db
     |> List.map ~f:(query ~convergent_species:phenotypes)
-    |> Fn.flip List.take 10
+    |> Fn.flip List.take 100
   in
   let gene_ids = List.map ~f:family_id queries in
   let ali_files = List.map ~f:Q.alignment queries in
@@ -614,7 +629,9 @@ let rer_converge ?max_read ?transform ?weighted ?scale ?continuous ~db
     List.map
       ~f:(fun q ->
         let spp = iqtree_nexus_partition_file_of_alignment q in
-        let clipped_tree = clip_tree_on_alignment species_tree q in
+        let clipped_tree =
+          clip_tree_on_alignment ~handle_tags:false species_tree q
+        in
         Iqtree.iqtree ~te:clipped_tree (*~m:model*) ~st:`DNA ~spp
           (`phylip q)
         |> Iqtree.treefile)
@@ -624,4 +641,4 @@ let rer_converge ?max_read ?transform ?weighted ?scale ?continuous ~db
   let gene_tree_set = RER.gene_tree_file gene_trees in
   let phenotype = RER.phenotypes_of_convergent_species phenotypes in
   RER.rer_converge ?max_read ?transform ?weighted ?scale ?continuous
-    ~gene_tree_set ~phenotype
+    ~min_specs:10 ~gene_tree_set ~phenotype ~master_tree:species_tree

lib/rer_converge.ml

@@ -3,7 +3,7 @@ open Bistro
 open Bistro.Shell_dsl
 
 let img =
-  [ docker_image ~account:"lsdch" ~name:"rerconverge" ~tag:"latest" () ]
+  [ docker_image ~account:"lsdch" ~name:"rerconverge" ~tag:"0.1.0" () ]
 
 let phenotypes_of_convergent_species species_list =
   let f =
@@ -38,9 +38,13 @@ let template_transform t =
 
 let template_bool t = string (if t then "T" else "F")
 
-let script ?max_read ~transform ~weighted ~scale ~continuous ~gene_tree_set
-    ~phenotype =
-  let max_read = Option.value_map max_read ~default:(string "NA") ~f:int in
+let ropt f =
+  Option.value_map ~f ~default:(string "NA")
+
+let script ?max_read ?min_specs ~transform ~weighted ~scale ~continuous ~gene_tree_set
+    ~phenotype ~master_tree =
+  let max_read = ropt int  max_read in
+  let min_specs = ropt int min_specs in
   [%script
     {|
 
@@ -49,12 +53,13 @@ suppressPackageStartupMessages({
   library(RERconverge)
 })
 
-gene_trees <- readTrees({{dep gene_tree_set}}, max.read={{max_read}})
+masterTree <- read.tree(file="{{dep master_tree}}")
+gene_trees <- readTrees("{{dep gene_tree_set}}", max.read={{max_read}}, masterTree=masterTree, minSpecs={{min_specs}})
 
 relative_rates <- getAllResiduals(
   gene_trees,
   useSpecies = useSpecies,
-  transform = {{template_transform transform}},
+  transform = "{{template_transform transform}}",
   weighted = {{template_bool weighted}},
   scale = {{template_bool scale}}
 )
@@ -109,11 +114,11 @@ res <- if ({{template_bool continuous}}) {
   )
 }
 
-write_tsv({{dest}}, res)
+write_tsv("{{dest}}", res)
 |}]
 
-let rer_converge ?max_read ?(transform = `sqrt) ?(weighted = true)
-    ?(scale = true) ?(continuous = false) ~gene_tree_set ~phenotype =
+let rer_converge ?max_read ?min_specs ?(transform = `sqrt) ?(weighted = true)
+    ?(scale = true) ?(continuous = false) ~master_tree ~gene_tree_set ~phenotype =
   Bistro_utils.R_script.workflow ~img ~descr:"rer_converge"
     (script ?max_read ~transform ~weighted ~scale ~continuous
-       ~gene_tree_set ~phenotype)
+       ?min_specs ~gene_tree_set ~phenotype ~master_tree)

lib/rer_converge.mli

@@ -8,10 +8,12 @@ val gene_tree_file : (string * newick file) list -> tsv file
 
 val rer_converge :
   ?max_read:int ->
+  ?min_specs:int ->
   ?transform:[ `none | `sqrt | `log ] ->
   ?weighted:bool ->
   ?scale:bool ->
   ?continuous:bool ->
+  master_tree:newick file ->
   gene_tree_set:tsv file ->
   phenotype:tsv file ->
   tsv file