Commit 9e9674d1 authored by Louis Duchemin's avatar Louis Duchemin
Browse files

Mutsel CpG : hypermutability applies to convergent branches only

parent 583bb343
......@@ -17,7 +17,7 @@ let most_probable_aa (pref : Amino_acid.vector) =
| None -> assert false
| Some (_, i) -> i
let make_profiles n_h0 n_ha fitness_profiles =
let make_profiles n_h0 n_ha fitness_profiles =
let random_profile () = Random.int (Array.length fitness_profiles)
and most_probable_aa = Array.map fitness_profiles ~f:most_probable_aa in
let h0_profiles = Array.init n_h0 ~f:(fun _ -> random_profile ())
......@@ -32,11 +32,11 @@ let make_profiles n_h0 n_ha fitness_profiles =
) in
h0_profiles, ha_profiles
let make_params ~gBGC:(gBGC0, gBGC1) ~ne_s:(ne_s0, ne_s1) ~fitness_profiles ~h0_profiles ~ha_profiles =
let make_params ~gBGC:(gBGC0, gBGC1) ~ne_s:(ne_s0, ne_s1) ~fitness_profiles ~h0_profiles ~ha_profiles =
let base_param = Mutsel.random_param ~alpha_nucleotide:10. ~alpha_fitness:0.1 in
let make_param ~gBGC ~ne_s p =
let make_param ~gBGC ~ne_s p =
{
base_param with
base_param with
gBGC; scaled_fitness = rescale_fitness ne_s fitness_profiles.(p)
}
in
......@@ -44,18 +44,18 @@ let make_params ~gBGC:(gBGC0, gBGC1) ~ne_s:(ne_s0, ne_s1) ~fitness_profiles ~h0_
make_param ~gBGC:gBGC0 ~ne_s:ne_s0 p,
make_param ~gBGC:gBGC1 ~ne_s:ne_s1 q
in
let h0_params = Array.map h0_profiles ~f:(fun profile ->
let h0_params = Array.map h0_profiles ~f:(fun profile ->
make_param_pair profile profile
)
and ha_params = Array.map ha_profiles ~f:(fun (p,q) ->
make_param_pair p q
and ha_params = Array.map ha_profiles ~f:(fun (p,q) ->
make_param_pair p q
)
in
(base_param, h0_params, ha_params)
let params_of_cond h0_params ha_params =
let n_h0 = Array.length h0_params in
fun i cond ->
fun i cond ->
let p, q =
if i < n_h0 then h0_params.(i)
else ha_params.(i - n_h0)
......@@ -93,7 +93,7 @@ module Site_independent = struct
let make ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~tree ~fitness_profiles () =
let h0_profiles, ha_profiles = make_profiles n_h0 n_ha fitness_profiles in
let base_param, h0_params, ha_params =
let base_param, h0_params, ha_params =
make_params ~gBGC ~ne_s ~fitness_profiles ~h0_profiles ~ha_profiles in
let params_of_cond = params_of_cond h0_params ha_params in
let params i bi = params_of_cond i (Branch_info.condition bi) in
......@@ -110,7 +110,7 @@ module Site_independent = struct
gBGC ; ne_s ; h0_params ; ha_params ; base_param ; tree}
end
module With_CpG_hypermutability = struct
module With_CpG_hypermutability = struct
module Sim = Mutsel_simulator_cpg.Make(Branch_info)
......@@ -127,14 +127,18 @@ module With_CpG_hypermutability = struct
let make ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~rate_CpG ~(tree:Convergence_tree.t) ~fitness_profiles () =
let h0_profiles, ha_profiles = make_profiles n_h0 n_ha fitness_profiles in
let _, h0_params, ha_params =
let _, h0_params, ha_params =
make_params ~gBGC ~ne_s ~fitness_profiles ~h0_profiles ~ha_profiles in
let params_of_cond = params_of_cond h0_params ha_params in
let param_factory state pos branch =
let module EM = Mutsel_simulator_cpg.Evolution_model in
let left_codon, right_codon = EM.flanking_codons state pos in
let context_CpG = EM.context_CpG_of_codons ?left_codon ?right_codon () in
let mutsel = params_of_cond pos (Branch_info.condition branch) in
let mutsel = params_of_cond pos (Branch_info.condition branch)
and rate_CpG = match Branch_info.condition branch with
| `Ancestral -> 1.
| `Convergent -> rate_CpG
in
EM.{ mutsel ; rate_CpG ; context_CpG }
in
(* FIXME : this is stationary dist for Mutsel without CpG *)
......@@ -144,24 +148,24 @@ module With_CpG_hypermutability = struct
|> Mutsel.NSCodon.Table.of_array_exn
)
in
let root = Array.map root_dists ~f:(fun codon_dist ->
let root = Array.map root_dists ~f:(fun codon_dist ->
let dist = Gsl.Randist.discrete_preproc (codon_dist :> float array) in
Gsl.Randist.discrete rng dist
Gsl.Randist.discrete rng dist
|> Mutsel.NSCodon.of_int_exn
)
in
let sequences =
let sequences =
Sim.sequence_gillespie_direct rng tree ~root ~param:param_factory
|> Tree.leaves
|> List.map ~f:Dna.of_codons
|> Tree.leaves
|> List.map ~f:Dna.of_codons
in
{ sequences ; fitness_profiles ; h0_profiles ; ha_profiles ;
gBGC ; ne_s ; tree ; rate_CpG}
end
type t =
Site_independent of Site_independent.t
type t =
Site_independent of Site_independent.t
| With_CpG_hypermutability of With_CpG_hypermutability.t
let make_with_CpG_hypermutability ~rng ~n_h0 ~n_ha ~ne_s ~gBGC ~rate_CpG ~tree ~fitness_profiles () =
......
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