Gemma pipeline implementation : relatedness matrix

lib/gemma.ml

@@ -11,6 +11,29 @@ let%pworkflow genotype_of_fasta fasta =
   let module G = Reviewphiltrans_toolbox.Gemma in
   G.write_genotypes ~alignment:[%path fasta] ~output:[%dest]
 
+let%pworkflow phenotype_of_tree nhx =
+  let open Phylogenetics in
+  let module U = Reviewphiltrans_toolbox.Utils in 
+  let collect_leaves t =
+    let rec node condition t acc =
+      match t with 
+      | Tree.Node n -> 
+        List1.fold_right n.branches ~init:acc ~f:branch
+      | Leaf _ -> condition :: acc
+    and branch (Tree.Branch b) acc =
+      node (snd b.data) b.tip acc
+    in
+    node 0 t []
+  in
+  let leaves =
+    U.tree_from_file [%path nhx]  
+    |> collect_leaves
+  in
+  let write_phenotypes leaves oc =
+    List.iter leaves ~f:(fprintf oc "%d\n")
+  in
+  Out_channel.with_file [%dest] ~f:(write_phenotypes leaves)
+
 let template_of_mode m =
   int (
     match m with
@@ -18,12 +41,20 @@ let template_of_mode m =
     | `Standardized -> 2
   )
 
+let relatedness_filename = function
+  | `Centered -> "result.cXX.txt"
+  | `Standardized -> "result.sXX.txt"
+
 let calculate_relatedness_matrix ~mode ~genotype ~phenotype =
   Workflow.shell ~descr:"gemma.calculate_relatedness_matrix" [
     cmd "gemma" [
       opt "-g" dep genotype ;
       opt "-p" dep phenotype ;
       opt "-gk" template_of_mode mode ;
-      opt "-o" Fn.id dest ;
-    ]
+      opt "-outdir" Fn.id tmp ;
+    ] ;
+    cmd "cp" [
+      tmp // (relatedness_filename mode) ;
+      dest ;
+    ] ;
   ]
\ No newline at end of file

lib/gemma.mli

@@ -7,6 +7,7 @@ class type relatedness_matrix = object
 end
 
 val genotype_of_fasta : aminoacid_fasta file -> bimbam file
+val phenotype_of_tree : nhx file -> bimbam file
 
 val calculate_relatedness_matrix :
   mode:[`Centered | `Standardized] ->