lmm_benchmark.ml 2.46 KB
Newer Older
Philippe Veber's avatar
Philippe Veber committed
1 2
open Core_kernel
open Bistro
3
open Codepi
Philippe Veber's avatar
Philippe Veber committed
4
open Codepi.File_formats
5

6 7
module Pipeline = Simulation_pipeline.Mutsel

Philippe Veber's avatar
Philippe Veber committed
8
type dataset = {
9
  label : string ;
Philippe Veber's avatar
Philippe Veber committed
10 11 12 13 14
  tree : nhx file ;
  rooted : bool ;
  branch_scale : float ;
  ne_s : float * float ;
}
15

Philippe Veber's avatar
Philippe Veber committed
16
let besnard2009 = {
17
  label = "besnard2009" ;
Philippe Veber's avatar
Philippe Veber committed
18 19
  tree = Bistro.Workflow.input "data/besnard2009/besnard2009.nhx" ;
  rooted = true ;
20
  branch_scale = 1. ;
21
  ne_s = 4., 4. ;
Philippe Veber's avatar
Philippe Veber committed
22 23 24
}

let oneline_rodent = {
25
  label = "online_rodent" ;
Philippe Veber's avatar
Philippe Veber committed
26 27
  tree = Bistro.Workflow.input "data/online_rodent/online_rodent.nhx"  ;
  rooted = true ;
28
  branch_scale = 1. ;
29
  ne_s = 4., 4. ;
Philippe Veber's avatar
Philippe Veber committed
30 31 32
}

let rubisco = {
33
  label = "rubisco" ;
Philippe Veber's avatar
Philippe Veber committed
34 35
  tree = Rubisco_dataset.(Path "data/rubisco" |> Query.tree ~branch_length_unit:`Amino_acid) ;
  rooted = false ;
36
  branch_scale = 1. ;
37
  ne_s = 4., 4. ;
Philippe Veber's avatar
Philippe Veber committed
38 39 40
}

let orthomam_echolocation = {
41
  label = "orthomam_echolocation" ;
Philippe Veber's avatar
Philippe Veber committed
42 43 44 45 46 47 48
  tree = (
    Orthomam.tree_of_db
      ~branch_length_unit:`Amino_acid
      ~convergent_species:Orthomam.species_with_echolocation
      (Codepitk.Orthomam_db.make "omm")
  ) ;
  rooted = false ;
49
  branch_scale = 1. ;
Philippe Veber's avatar
Philippe Veber committed
50 51 52 53
  ne_s = 4., 4. ;
}

type detection_method = {
54
  result : Pipeline.query -> cpt file ;
Philippe Veber's avatar
Philippe Veber committed
55 56 57 58
  label : string ;
  requires_rooted_tree : bool ;
}

59 60
let meth ?(requires_rooted_tree = false) result label =
  { result ; label ; requires_rooted_tree }
Philippe Veber's avatar
Philippe Veber committed
61

62
let methods = Pipeline.[
Philippe Veber's avatar
Philippe Veber committed
63
  meth tdg09 "tdg09" ;
64
  meth pcoc "pcoc" ;
65
  (* meth pcoc_v2 ~col:3 "pcoc v2" ; *)
Philippe Veber's avatar
Philippe Veber committed
66 67 68 69 70 71
  meth (gemma ~lmm_test:`Score ~relatedness_mode:`Standardized) "gemma" ;
  meth inhouse_lmm "LMM" ;
  meth multinomial_asymptotic_lrt "multinomial" ;
  meth topological "topological" ~requires_rooted_tree:true ;
]

72
let benchmark_rds ?(seed = 42) { tree = t ; rooted ; ne_s ; branch_scale ; _ } =
73 74
  let q =
    Pipeline.query ~seed ~tree:(NHX t) ~branch_scale ~ne_s
75 76 77
      ~profiles:"example/aa_fitness/263SelectedProfiles.tsv"
      ~n_h0:900 ~n_ha:100 ()
  in
78
  let simulation = Pipeline.simulation q in
79
  let results =
80 81
    List.filter_map methods ~f:(fun m ->
        if not m.requires_rooted_tree || rooted then
82
          Some (m.result q)
83 84 85
        else None
      )
  in
86
  Pipeline.benchmark_statistics simulation ~results
Philippe Veber's avatar
Philippe Veber committed
87 88

let () =
89 90 91
  let open Bistro_utils.Repo in
  let datasets = [ besnard2009 ; rubisco ; oneline_rodent ; orthomam_echolocation ] in
  let repo =
92
    List.map datasets ~f:(fun d ->
93 94 95
        item [d.label ^ ".rds"] (benchmark_rds d)
      )
  in
96
  let loggers = [ Bistro_utils.Console_logger.create () ] in
97
  build_main ~loggers ~np:4 ~mem:(`GB 4) repo ~outdir:"res"