Commit 10553ec5 authored by Louis Duchemin's avatar Louis Duchemin
Browse files

Adds dinucleotide ratio computation in Mutsel_simulator_cpg

parent 21f5f943
......@@ -115,6 +115,28 @@ module Evolution_model = struct
)
let stationary_distribution p = NSCodon_rate_matrix.stationary_distribution (rate_matrix p)
let dinucleotide_ratios (sequence:NSCodon.t array) =
let seq = Array.concat_map ~f:(fun codon ->
let x, y, z = NSCodon.nucleotides codon in
[|x;y;z|]
) sequence
in
let single_counts = Nucleotide.Table.init (fun nuc ->
Array.count seq ~f:(fun seq_nuc -> Nucleotide.equal seq_nuc nuc))
in
Nucleotide.Table.init (fun nuc1 ->
Nucleotide.Table.init (fun nuc2 ->
let expected = Nucleotide.Table.(get single_counts nuc1 * get single_counts nuc2)
in
let observed = Array.counti seq ~f:(fun pos nuc ->
if pos = Array.length seq then false
else Nucleotide.(equal nuc nuc1 && equal seq.(pos+1) nuc2)
)
in
float_of_int observed /. float_of_int expected
)
)
end
module Make(BI: Simulator.Branch_info) = Simulator.Make(NSCodon)(Evolution_model)(BI)
......@@ -35,6 +35,8 @@ module Evolution_model : sig
val flat_param: unit -> NSCodon.t array -> int -> 'a -> param
val random_param: alpha_nuc:float -> alpha_fitness:float -> NSCodon.t array -> int -> 'a -> param
val dinucleotide_ratios: NSCodon.t array -> float Nucleotide.table Nucleotide.table
end
module Make(BI: Simulator.Branch_info) : sig
......
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