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

notes/Codepi_multinomial_note: added stuff on tdg09

parent 9aaf878a
......@@ -10,6 +10,13 @@ We investivate here how this simple approach behaves compared to more
sophisticated models that take the species phylogeny into account.
```ocaml
open Core;;
open Phylogenetics;;
open Codepi;;
open Codepitk;;
module Tk = Codepitk;;
module Pipeline = Simulation_pipeline.Mutsel
module Top = Bistro_utils.Toplevel_eval.Make(
struct
let np = 8
......@@ -38,11 +45,7 @@ module Show = struct
end
;;
open Core;;
open Top;;
open Codepi;;
module Tk = Codepitk;;
module Pipeline = Simulation_pipeline.Mutsel
```
For that purpose we use simulations under the mutsel model.
......@@ -53,11 +56,12 @@ let n_ha = 100
let tree_path = "data/besnard2009/besnard2009.nhx"
let tree = Simulation_pipeline.NHX (Bistro.Workflow.input tree_path)
let profiles = "example/aa_fitness/263SelectedProfiles.tsv"
let sim = Pipeline.query ~tree ~profiles ~n_h0 ~n_ha ()
let sim =
Pipeline.query ~tree ~profiles ~n_h0 ~n_ha ~ne_s:(2., 2.) ()
```
```ocaml
let benchmark1 =
let benchmark1 =
Pipeline.(benchmark sim [ multinomial_asymptotic_lrt ; tdg09 ])
;;
......@@ -213,3 +217,80 @@ Show.cached_png_rplot (fun%workflow () ->
## Simulating under tdg09's null distribution
```ocaml
let wag =
Bistro_unix.wget "https://www.ebi.ac.uk/goldman-srv/WAG/wag.dat"
|> Top.path
|> Wag.parse
;;
let tdg09_site_simulation ?(alpha = 1.) ?(scale = 1.) tree =
let open Tk.Tdg09.Pack in
let exchangeability_matrix = wag.rate_matrix in
let stationary_distribution = simulate_profile alpha in
let param = scale in
Model1.simulate_site ~exchangeability_matrix ~stationary_distribution ~param tree
;;
let tdg09_simulation ?alpha tree nsites =
Array.init nsites ~f:(fun _ ->
tdg09_site_simulation tree ?alpha
)
;;
let multinomial_on_tdg09_simulation site =
let leaves1, leaves2 =
Convergence_tree.leaves site
|> List.partition_map ~f:(function
| aa, `Ancestral -> Either.first aa
| aa, `Convergent -> Either.second aa
)
in
let counts x = (Amino_acid.counts (Sequence.of_list x) :> int array) in
let x1 = counts leaves1 in
let x2 = counts leaves2 in
let d = Multinomial_test.data ~x1 ~x2 in
let r = Multinomial_test.(
Permutation.test ~statistic:LRT.likelihood_log_ratio d
) in
r.pvalue
;;
let tdg09_on_tdg09_simulation tree sim =
let open Tk.Tdg09.Pack in
let site = Tree.leaves sim |> Array.of_list in
let _, _, lrt = Model3.lrt ~mode:`sparse wag tree site in
lrt.pvalue
;;
let tdg09_uniformity_test ?alpha tree f =
let sites = tdg09_simulation ?alpha tree 1_000 in
Array.map sites ~f
;;
let _ =
Show.cached_png_rplot (fun%workflow () ->
let tree =
Tk.Tdg09.Pack.pair_tree
~branch_length1:0.3 ~branch_length2:0.3
~npairs:10 in
tdg09_uniformity_test ~alpha:0.1 tree multinomial_on_tdg09_simulation
|> OCamlR_graphics.hist ~breaks:(`n 60)
|> fun _ -> ()
)
;;
```
```ocaml
let _ =
Show.cached_png_rplot (fun%workflow () ->
let tree =
Tk.Tdg09.Pack.pair_tree
~branch_length1:0.3 ~branch_length2:0.3
~npairs:10 in
tdg09_uniformity_test ~alpha:0.1 tree (tdg09_on_tdg09_simulation tree)
|> OCamlR_graphics.hist ~breaks:(`n 60)
|> fun _ -> ()
)
;;
```
......@@ -2,6 +2,7 @@
(name codepi_multinomial_note)
(libraries codepi scirep)
(modules Codepi_multinomial_note)
(modes native)
(preprocess (pps bistro.ppx)))
(rule
......
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