Commit b6a50ad7 authored by Philippe Veber's avatar Philippe Veber
Browse files

update wrt phylogenetics

parent b4a5fdfc
......@@ -46,10 +46,12 @@ let nucleotide_fasta_gc ?pos fa =
seq_gc ?pos seqs
let nucleotide_fasta_gc_ac ?pos tree fa =
let module BI = Phylogenetics_convergence.Simulator.Branch_info in
let tree = Reviewphiltrans_toolbox.Utils.tree_from_file tree in
let seqs = strings_from_fasta fa in
let root = { BI.length = 0. ; condition = `Ancestral } in
let leaf_state =
Phylogenetics.Tree.map_leaves tree ~root:(0., 0) ~f:(fun _ b -> snd b = 0)
Phylogenetics.Tree.map_leaves tree ~root ~f:Poly.(fun _ b -> BI.condition b = `Ancestral)
|> Array.of_list
in
let ancestral_seqs, convergent_seqs =
......
......@@ -21,11 +21,16 @@ let%pworkflow phenotype_of_tree nhx =
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 []
and branch (Tree.Branch b) acc =
node (Phylogenetics_convergence.Simulator.Branch_info.condition b.data) b.tip acc
in
node `Ancestral 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
let write_phenotypes leaves oc = List.iter leaves ~f:(fun cond ->
fprintf oc "%d\n" (match cond with `Ancestral -> 0 | `Convergent -> 1)
)
in
Out_channel.with_file [%dest] ~f:(write_phenotypes leaves)
let template_of_lmm lmm =
......
......@@ -43,7 +43,8 @@ let%pworkflow multinomial_ocaml_implementation ~meth ~(tree_sc:_ file) ~(faa:ami
|> Rresult.R.get_ok
in
let tree = Reviewphiltrans_toolbox.Utils.tree_from_file [%path tree_sc] in
let leaves = fold_leaves tree ~init:[] ~f:(fun acc (_, cond) ni ->
let leaves = fold_leaves tree ~init:[] ~f:(fun acc bi ni ->
let cond = Phylogenetics_convergence.Simulator.Branch_info.condition bi in
match ni.name with
| None -> failwith "Leaves of the tree should be named"
| Some n -> (
......@@ -54,9 +55,8 @@ let%pworkflow multinomial_ocaml_implementation ~meth ~(tree_sc:_ file) ~(faa:ami
)
in
let seqs0, seqs1 = List.partition_map leaves ~f:Either.(function
| (aa, 0) -> First aa
| (aa, 1) -> Second aa
| _ -> assert false
| (aa, `Ancestral) -> First aa
| (aa, `Convergent) -> Second aa
)
in
let counts seqs i =
......
......@@ -2,6 +2,7 @@ open Core_kernel
let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1) ~tree ~fitness_profiles () =
let open Phylogenetics in
let module Convsim = Phylogenetics_convergence.Simulator in
let () = Option.iter ~f:Random.init [%param seed] in
let n_h0 = [%param n_h0] in
let n_ha = [%param n_ha] in
......@@ -50,27 +51,27 @@ let%pworkflow simulator ?branch_factor ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~g
loop ()
)
in
let params i cond =
let params_of_cond i cond =
let p, q =
if i < n_h0 then h0_params.(i)
else ha_params.(i - n_h0)
in
match cond with
| 0 -> p
| 1 -> q
| _ -> assert false
| `Ancestral -> p
| `Convergent -> q
in
let root_condition =
Option.value_exn (Phylogenetics_convergence.Simulator.root_condition tree)
let params i bi =
params_of_cond i (Convsim.Branch_info.condition bi)
in
let root_condition = `Ancestral in
let root_dists = Array.init (n_h0 + n_ha) ~f:(fun i ->
Mutsel.stationary_distribution (params i root_condition)
|> Mutsel.NSCodon.Vector.to_array
|> Mutsel.NSCodon.Table.of_array_exn
Mutsel.stationary_distribution (params_of_cond i root_condition)
|> Mutsel.NSCodon.Vector.to_array
|> Mutsel.NSCodon.Table.of_array_exn
)
in
let root = Phylogenetics_convergence.Simulator.Mutsel.hmm0 ~len:(n_h0 + n_ha) ~dist:(Array.get root_dists) in
let ali = Phylogenetics_convergence.Simulator.Mutsel.alignment tree ~root params in
let root = Convsim.Mutsel.hmm0 ~len:(n_h0 + n_ha) ~dist:(Array.get root_dists) in
let ali = Convsim.Mutsel.alignment tree ~root params in
let species_name =
Phylogenetics.Tree.leaves tree
|> List.map ~f:(fun { name } -> Option.value_exn name) in
......
(library
(name reviewphiltrans_toolbox)
(libraries biotk biocaml.ez ocaml-r.graphics ocaml-r.grDevices phylogenetics)
(libraries biotk biocaml.ez ocaml-r.graphics ocaml-r.grDevices phylogenetics phylogenetics.convergence)
(inline_tests
(deps ../../tests/data/gemma_output.tsv))
(preprocess
......
open Core_kernel
let tree_from_file ?(alpha = 1.) fn =
match Phylogenetics.Newick.from_file fn with
| Branch (Branch {tip= t; _}) | Tree t ->
Phylogenetics.Tree.map t ~node:Fn.id ~leaf:Fn.id
~branch:
Phylogenetics.Newick.(
fun b ->
let length = Option.value_exn b.length *. alpha in
let condition =
match
List.Assoc.find ~equal:String.equal b.tags "Condition"
with
| Some s -> Int.of_string s
| None -> failwith "Missing Condition tag"
in
(length, condition))
let open Phylogenetics in
let module BI = Phylogenetics_convergence.Simulator.Branch_info in
Newick.from_file fn
|> Newick.with_inner_tree ~f:(fun t ->
Tree.map t ~node:Fn.id ~leaf:Fn.id ~branch:Phylogenetics.Newick.(fun b ->
let length = Option.value_exn b.length *. alpha in
let condition =
match List.Assoc.find ~equal:String.equal b.tags "Condition" with
| Some s -> (
match s with
| "0" -> `Ancestral
| "1" -> `Convergent
| _ -> failwithf "Invalid condition: %s" s ()
)
| None -> failwith "Missing Condition tag"
in
{ BI.length ; condition }
)
)
let range_iter a b ~f =
let rec loop i =
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment