mutsel_simulator.ml 2.27 KB
Newer Older
Philippe Veber's avatar
Philippe Veber committed
1
open Core_kernel
Philippe Veber's avatar
Philippe Veber committed
2
open Bistro
Philippe Veber's avatar
Philippe Veber committed
3

Philippe Veber's avatar
Philippe Veber committed
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
let simulation ?branch_scale ?seed ~n_h0 ~n_ha ~ne_s:(ne_s0, ne_s1) ~gBGC:(gBGC0, gBGC1) ~tree ~fitness_profiles () =
  let f = fun%workflow () ->
    let open Phylogenetics in
    let () = Option.iter ~f:Random.init [%param seed] in
    let n_h0 = [%param n_h0] in
    let n_ha = [%param n_ha] in
    let ne_s0 = [%param ne_s0] in
    let ne_s1 = [%param ne_s1] in
    let gBGC0 = [%param gBGC0] in
    let gBGC1 = [%param gBGC1] in
    let branch_scale = [%param branch_scale] in
    let tree =
      Codepitk.Convergence_tree.from_file [%path tree]
      |> Rresult.R.get_ok
    in
    let fitness_profiles =
20
      Codepitk.Profile_tsv.(
Philippe Veber's avatar
Philippe Veber committed
21 22 23 24
        read [%path fitness_profiles]
        |> to_fitness
        |> Array.map ~f: Amino_acid.Vector.of_array_exn
      ) in
Louis Duchemin's avatar
Louis Duchemin committed
25 26 27 28 29
    let tree = Option.value_map branch_scale ~default:tree 
        ~f:(fun scale -> Codepitk.Mutsel_simulator.tree_rescale_branch_length tree ~scale)
    in
    Codepitk.Mutsel_simulator.Site_independent.make
      ?seed ~n_h0 ~n_ha ~gBGC:(gBGC0, gBGC1)
Philippe Veber's avatar
Philippe Veber committed
30
      ~ne_s:(ne_s0, ne_s1) ~tree ~fitness_profiles ()
31
  in
Philippe Veber's avatar
Philippe Veber committed
32
  Workflow.plugin ~descr:"simulator.simulation" f
33

Philippe Veber's avatar
Philippe Veber committed
34 35
let alignment_of_simulation sim =
  let f = fun%workflow dest ->
Louis Duchemin's avatar
Louis Duchemin committed
36
    let sim : Codepitk.Mutsel_simulator.Site_independent.t =
Philippe Veber's avatar
Philippe Veber committed
37 38 39 40 41 42 43
      [%eval sim]
    in
    let species_name = Phylogenetics.Tree.leaves sim.tree in
    Out_channel.with_file dest ~f:(fun oc ->
        List.iter2_exn species_name sim.sequences ~f:(fun description sequence ->
            fprintf oc ">%s\n%s\n" description (sequence :> string)
          )
44 45
      )
  in
Philippe Veber's avatar
Philippe Veber committed
46
  Workflow.path_plugin ~descr:"simulator.alignment_of_simulation" f
Philippe Veber's avatar
Philippe Veber committed
47

Philippe Veber's avatar
Philippe Veber committed
48 49 50
let fitness_histogram sim =
  let f = fun%workflow dest ->
    let open Phylogenetics in
Louis Duchemin's avatar
Louis Duchemin committed
51
    let sim : Codepitk.Mutsel_simulator.Site_independent.t =
Philippe Veber's avatar
Philippe Veber committed
52
      [%eval sim]
53
    in
Philippe Veber's avatar
Philippe Veber committed
54 55 56 57 58 59 60 61 62 63 64
    let params = Array.append sim.h0_params sim.ha_params in
    let data =
      Array.fold params ~init:[] ~f:(fun acc (p,q) ->
          (Amino_acid.Vector.to_array p.scaled_fitness) ::
          (Amino_acid.Vector.to_array q.scaled_fitness) :: acc
        )
      |> Array.concat
    in
    OCamlR_grDevices.pdf dest ;
    ignore (OCamlR_graphics.hist data : OCamlR_graphics.hist) ;
    OCamlR_grDevices.dev_off ()
65
  in
Philippe Veber's avatar
Philippe Veber committed
66
  Workflow.path_plugin ~descr:"simulator.fitness_histogram" f