Commit 84ebb3e3 authored by Louis Duchemin's avatar Louis Duchemin
Browse files

Mutsel CpG simulator demo

parent 983286d5
......@@ -59,18 +59,18 @@ module Evolution_model = struct
context_CpG;
}
let flat_param () =
let flat_param ~rate_CpG =
let pi = Nucleotide.flat_profile () in
let rho = Linear_algebra.Lacaml.Vector.init 6 ~f:(fun _ -> 1. /. 6.) in
let nucleotide_rates = Rate_matrix.Nucleotide.gtr ~equilibrium_frequencies:pi ~transition_rates:rho in
make_param ~gBGC:0. ~scaled_fitness:(flat_fitness ()) ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG:1.
make_param ~gBGC:0. ~scaled_fitness:(flat_fitness ()) ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG
let random_param ~alpha_nuc ~alpha_fitness =
let random_param ~alpha_nuc ~alpha_fitness ~rate_CpG =
let pi = Nucleotide.random_profile alpha_nuc in
let rho = Utils.random_profile 6 in
let nucleotide_rates = Rate_matrix.Nucleotide.gtr ~equilibrium_frequencies:pi ~transition_rates:rho in
let scaled_fitness = Amino_acid.random_profile alpha_fitness |> fitness_of_profile in
make_param ~gBGC:0. ~scaled_fitness ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG:1.
make_param ~gBGC:0. ~scaled_fitness ~omega:1. ~nucleotide_rates ~nucleotide_stat_dist:pi ~rate_CpG
let rate_increase_CpG (pos, x_a, x_b) context =
match Nucleotide.(pos, inspect x_a, inspect x_b), context with
......@@ -158,3 +158,41 @@ module Evolution_model = struct
end
module Make(BI: Simulator.Branch_info) = Simulator.Make(NSCodon)(Evolution_model)(BI)
let demo seq_length ~rate_CpG ~branch_length =
let module Branch_info = struct
type t = Newick.branch_info
let length branch = Option.value_exn Newick.(branch.length)
end
in
let module M = Make(Branch_info) in
let root_sequence = Array.init seq_length ~f:(fun _ ->
let codon_seed = Owl.Stats.uniform_int_rvs ~a:0 ~b:(NSCodon.card - 1) in
match NSCodon.of_int codon_seed with
| Some codon -> codon
| None -> sprintf "Cannot create codon from integer : %d" codon_seed |> failwith
) in
let param = Evolution_model.random_param ~alpha_nuc:0.5 ~alpha_fitness:0.5 ~rate_CpG in
let tree =
sprintf "(leaf:%f);" branch_length
|>Newick.from_string
|>Newick.with_inner_tree ~f:(fun tree ->
M.sequence_gillespie_direct tree ~root:root_sequence ~param
)
in
let simulated_sequence = Tree.leaves tree |> List.hd_exn in
let counts = Evolution_model.dinucleotide_ratios simulated_sequence in
Array.init Nucleotide.card ~f:(fun i ->
let nuc = Nucleotide.of_int_exn i |> Nucleotide.to_char |> Char.to_string in
"\t" ^ nuc
)
|> String.concat_array ~sep:"\t"
|> print_endline ;
Array.iteri counts ~f:(fun i nuc_ratios ->
let nuc = Nucleotide.of_int_exn i |> Nucleotide.to_char |> Char.to_string in
print_string nuc;
Array.iteri nuc_ratios ~f:(fun _ ratio ->
print_string ("\t" ^ string_of_float ratio)
);
Out_channel.newline stdout;
)
......@@ -32,9 +32,9 @@ module Evolution_model : sig
rate_CpG:float ->
NSCodon.t array -> int -> 'a -> param
val flat_param: unit -> NSCodon.t array -> int -> 'a -> param
val flat_param: rate_CpG:float -> NSCodon.t array -> int -> 'a -> param
val random_param: alpha_nuc:float -> alpha_fitness:float -> NSCodon.t array -> int -> 'a -> param
val random_param: alpha_nuc:float -> alpha_fitness:float -> rate_CpG:float -> NSCodon.t array -> int -> 'a -> param
val dinucleotide_ratios: NSCodon.t array -> float array array
end
......@@ -42,3 +42,5 @@ end
module Make(BI: Simulator.Branch_info) : sig
include module type of Simulator.Make(NSCodon)(Evolution_model)(BI)
end
val demo : int -> rate_CpG:float -> branch_length:float -> unit
\ No newline at end of file
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