RERconverge : refactoring ; add relative rates plot

app/wip.ml

@@ -2,28 +2,35 @@ let db = Codepitk.Orthomam_db.make "/home/louis/Data/omm/"
 
 let echolocators = Bistro.Workflow.data Codepi.Orthomam.species_with_echolocation
 
-let subterranean = 
+let subterraneans = 
   Bistro.Workflow.input "/home/louis/dev/codepi/data/subterranean.tsv"
   |> Codepi.Orthomam.Binary_phenotype.load_convergent_species
 
-let aquatic = 
-  Bistro.Workflow.input "/home/louis/dev/codepi/data/orthomam_aquatic.tsv"
+let marines = 
+  Bistro.Workflow.input "/home/louis/dev/codepi/data/orthomam_marine.tsv"
   |> Codepi.Orthomam.Binary_phenotype.load_convergent_species
 
-let target =
-  Codepi.Orthomam.rer_converge ~db ~phenotypes:echolocators ()
+let pipeline ~phenotypes =
+  Codepi.Orthomam.rer_converge ~db ~phenotypes ()
 
 let () =
   try
     let loggers = [ Bistro_utils.Console_logger.create () ] in
     let allowed_containers = [ `Docker ] in
     let execute workflow = 
-    Bistro.Workflow.path workflow
-    |> Bistro_engine.Scheduler.simple_eval_exn ~allowed_containers ~np:10
-      ~mem:(`GB 8) ~loggers
-    |> print_endline 
-  in
-    execute (fst target);
-    execute (snd target);
+      Bistro.Workflow.path workflow
+      |> Bistro_engine.Scheduler.simple_eval_exn ~allowed_containers ~np:10
+        ~mem:(`GB 8) ~loggers
+      |> print_endline 
+    in
+    List.iter2 (fun label phenotypes -> 
+        let target = pipeline ~phenotypes in
+        print_endline label;
+        execute (target.rer_plot ~pat:"*LIM2*")
+      ) 
+      ["Subterraneans"]
+      (* ; "Marines "; "Echolocators"]  *)
+      [subterraneans]
+      (* ; marines; echolocators] *)
   with Failure msg -> print_endline msg
 

lib/R/plot_rers.R

+suppressPackageStartupMessages({
+  library(tidyverse)
+  library(RERconverge)
+})
+
+
+rds = readRDS("<<<dep rds>>>")
+gene_ids = <<<gene_ids>>>
+rers = rds$relative_rates
+
+pdf("<<<dest>>>")
+for (gene_id in gene_ids) {
+  plotRers(rers, gene_id, phenv=rds$traits_paths)
+}
+dev.off()
\ No newline at end of file

lib/R/rer_converge.R

@@ -4,7 +4,7 @@ suppressPackageStartupMessages({
 })
 
 master_tree <- read.tree(file="<<<dep master_tree>>>")
-gene_trees <- readTrees("<<<dep gene_tree_set>>>", max.read=<<<max_read>>>, master_tree=master_tree, minSpecs=<<<min_specs>>>)
+gene_trees <- readTrees("<<<dep gene_tree_set>>>", max.read=<<<max_read>>>, masterTree=master_tree, minSpecs=<<<min_specs>>>)
 
 relative_rates <- getAllResiduals(
   gene_trees,

lib/orthomam.ml

@@ -12,7 +12,7 @@ let alignments_of_db db =
   List.map (Orthomam_db.list_alignments db) ~f:(fun ali -> Alignment (db, ali))
 
 let search_alignments ~pat db =
-  let re = Re.(compile (Glob.glob pat)) in
+  let re = Re.(compile (Glob.glob ~expand_braces:true pat)) in
   alignments_of_db db
   |> List.filter ~f:(fun (Alignment (_, ali)) ->
       Re.execp re (Orthomam_db.id_of_alignment ali)
@@ -664,6 +664,12 @@ module Binary_phenotype = struct
     Workflow.plugin ~descr:"load_binary_phenotypes_tsv" f
 end
 
+type rer_converge_output = {
+  result_table : tsv file ;
+  best_candidate_summary : pdf file ;
+  rer_plot : pat:string -> pdf file ;
+}
+
 let rer_converge ?transform ?weighted ?scale ?continuous ?(summary_n_genes = 10) ~db
     ~phenotypes () =
   let open Bistro_bio in
@@ -701,9 +707,17 @@ let rer_converge ?transform ?weighted ?scale ?continuous ?(summary_n_genes = 10)
   let phenotypes = RER.phenotypes_of_convergent_species phenotypes in
   let rds = RER.rer_converge ?transform ?weighted ?scale ?continuous
     ~min_specs:10 ~master_tree:species_tree ~gene_tree_set ~phenotypes () in
-  let results_table = RER.result_table_of_rds rds in
-  let candidates_summary = RER.best_candidate_summary ~n_genes:summary_n_genes rds
+  let result_table = RER.result_table_of_rds rds in
+  let best_candidate_summary = RER.best_candidate_summary ~n_genes:summary_n_genes rds in
+  let rer_plot ~pat = 
+    search_alignments ~pat db
+    |> List.map ~f:(fun (Alignment (_, ali)) -> Orthomam_db.id_of_alignment ali)
+    |> RER.plot_relative_rates rds
   in
-  results_table, candidates_summary
+  {
+    result_table;
+    best_candidate_summary;
+    rer_plot;
+  }
 
 

lib/orthomam.mli

@@ -91,6 +91,12 @@ val site_ranking :
 
 val draw_site : query -> int -> pdf file
 
+type rer_converge_output = {
+  result_table : tsv file ;
+  best_candidate_summary : pdf file ;
+  rer_plot : pat:string -> pdf file ;
+}
+
 val rer_converge :
   ?transform:[ `none | `sqrt | `log ] ->
   ?weighted:bool ->
@@ -100,7 +106,7 @@ val rer_converge :
   db:Orthomam_db.t ->
   phenotypes:string list workflow ->
   unit ->
-  tsv file * pdf file
+  rer_converge_output
 
 module Binary_phenotype: sig
   type t = {

lib/rer_converge.ml

@@ -177,6 +177,12 @@ let match_species_tree_position ~gene_tree ~clipped_species_tree =
   in
   Workflow.path_plugin ~descr:"match_species_tree_position" f
 
+  let plot_relative_rates rds gene_ids  =
+    let module R = Bistro_utils.R_script in
+    let gene_ids = R.(make [strings gene_ids]) in
+    let script = [%include_script "lib/R/plot_rers.R"] in
+    R.workflow ~img ~descr:"RER_converge.plot_relative_rates" script
+
 
   let best_candidate_summary ?(n_genes = 10) rds=
     let n_genes = int n_genes in

lib/rer_converge.mli

@@ -29,6 +29,11 @@ val result_table_of_rds : rds file -> tsv file
 val match_species_tree_position :
   gene_tree:newick file -> clipped_species_tree:newick file -> newick file
 
+val plot_relative_rates : 
+  rds file ->
+  string list -> 
+  pdf file
+
 val best_candidate_summary : 
   ?n_genes:int ->
   rds file ->