multinomial.ml 1.75 KB
Newer Older
1
open Core
Philippe Veber's avatar
Philippe Veber committed
2
open Bistro
3 4
open File_formats

5
let%pworkflow multinomial_ocaml_implementation ~meth ~(tree_sc:_ file) ~(faa:aminoacid_fasta file) (* : cpt file *) =
6
  let open Phylogenetics in
Philippe Veber's avatar
Philippe Veber committed
7 8
  let module MT = Phylogenetics_convergence.Multinomial_test in
  let open Codepitk in
Philippe Veber's avatar
Philippe Veber committed
9 10 11 12 13 14 15
  let meth = [%param meth] in
  let test = match meth with
    | `Asymptotic_LRT -> MT.LRT.asymptotic_test
    | `Simulation_LRT -> MT.LRT.simulation_test ~sample_size:10_000
    | `Asymptotic_sparse -> MT.Sparse.asymptotic_test
    | `Simulation_sparse -> MT.Sparse.simulation_test ~sample_size:10_000
  in
16 17 18 19
  let alignment =
    Alignment.from_fasta [%path faa]
    |> Rresult.R.get_ok
  in
Philippe Veber's avatar
Philippe Veber committed
20 21 22
  let tree =
    Convergence_tree.from_file [%path tree_sc]
    |> Rresult.R.failwith_error_msg
23
  in
Philippe Veber's avatar
Philippe Veber committed
24
  let site c0 c1 =
25 26
    let c0 = (c0 : int Amino_acid.table :> int array) in
    let c1 = (c1 : int Amino_acid.table  :> int array) in
27
    let d = MT.data ~x1:c0 ~x2:c1 in
Philippe Veber's avatar
Philippe Veber committed
28
    let r = test d in
29 30
    r._T_, r.pvalue
  in
31
  let res = Convergence_tree.alignment_counts_map tree alignment site in
32 33
  let res_lines =
    List.mapi res ~f:(fun i (lr, pval) ->
Philippe Veber's avatar
Philippe Veber committed
34
        sprintf "%d\t%f\t%f" (i + 1) (1. -. pval) lr
35 36
      )
  in
37
  let header = "Sites\tMultinomial_1mp\tMultinomial_LikelihoodRatio" in
38
  Out_channel.write_lines [%dest] (header :: res_lines)
Philippe Veber's avatar
Philippe Veber committed
39 40 41 42 43

let multinomial_asymptotic_lrt ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Asymptotic_LRT ~tree_sc ~faa
let multinomial_simulation_lrt ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Simulation_LRT ~tree_sc ~faa
let multinomial_asymptotic_sparse ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Asymptotic_sparse ~tree_sc ~faa
let multinomial_simulation_sparse ~tree_sc ~faa = multinomial_ocaml_implementation ~meth:`Simulation_sparse ~tree_sc ~faa