Commit 04acf267 authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Mutsel_simulator_cpg: fixed context definition

context should not include information on a codon, since it is used to
compute a rate matrix between all codons.
parent b4e42135
......@@ -8,7 +8,6 @@ module Evolution_model = struct
type context_CpG = {
lh: Nucleotide.repr option;
rh: Nucleotide.repr option;
codon: (Nucleotide.repr * Nucleotide.repr * Nucleotide.repr)
}
type param = {
......@@ -25,12 +24,11 @@ module Evolution_model = struct
let x, y, z = NSCodon.nucleotides codon in
Nucleotide.(inspect x, inspect y, inspect z)
let context_CpG_of_codons p ?lh ?rh () =
let context_CpG_of_codons ?left_codon ?right_codon () =
let open Option.Monad_infix in
{
lh = lh >>| inspect_codon >>| trd3 ;
rh = rh >>| inspect_codon >>| fst3 ;
codon = inspect_codon p ;
lh = left_codon >>| inspect_codon >>| trd3 ;
rh = right_codon >>| inspect_codon >>| fst3 ;
}
let fitness_of_profile ?(beta = 1.) =
......@@ -43,11 +41,11 @@ module Evolution_model = struct
let make_param ~gBGC ~scaled_fitness ~omega ~nucleotide_rates ~nucleotide_stat_dist ~rate_CpG state pos _branch =
let context_CpG =
if pos = 0 then
context_CpG_of_codons state.(pos) ~rh:state.(pos+1) ()
context_CpG_of_codons ~right_codon:state.(pos+1) ()
else if pos = (Array.length state - 1 ) then
context_CpG_of_codons state.(pos) ~lh:state.(pos-1) ()
context_CpG_of_codons ~left_codon:state.(pos-1) ()
else
context_CpG_of_codons state.(pos) ~lh:state.(pos -1) ~rh:state.(pos + 1) ()
context_CpG_of_codons ~left_codon:state.(pos -1) ~right_codon:state.(pos + 1) ()
in
{
nucleotide_rates;
......@@ -72,14 +70,15 @@ module Evolution_model = struct
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
let rate_increase_CpG (pos, x_a, x_b) context =
match Nucleotide.(pos, inspect x_a, inspect x_b), context with
| (0, C, T), { codon = (_,G, _) ; _ }
| (1, C, T), { codon = (_,_, G) ; _ }
| (2, C, T), { rh = Some G ; _ }
| (0, G, A), { lh = Some C ; _ }
| (1, G, A), { codon = (C, _, _) ; _ }
| (2, G, A), { codon = (_,C,_) ; _ } -> true
let rate_increase_CpG (pos, x_a, x_b) context codon =
let open Nucleotide in
match (pos, inspect x_a, inspect x_b), context, codon with
| (0, C, T), _, (_, G, _)
| (1, C, T), _, (_, _, G)
| (2, C, T), { rh = Some G ; _ }, _
| (0, G, A), { lh = Some C ; _ }, _
| (1, G, A), _, (C, _, _)
| (2, G, A), _, (_,C,_) -> true
| _ -> false
let fixation_probability delta =
......@@ -94,7 +93,7 @@ module Evolution_model = struct
Mutsel.NSCodon_rate_matrix.make (fun p q ->
match NSCodon.neighbours p q with
| Some (pos, x_a, x_b) ->
let rate_CpG = if rate_increase_CpG (pos, x_a, x_b) context_CpG then rate_CpG else 1. in
let rate_CpG = if rate_increase_CpG (pos, x_a, x_b) context_CpG (inspect_codon p) then rate_CpG else 1. in
let _B_ = match Nucleotide.(inspect x_a, inspect x_b) with
| (A | T), (C | G) -> gBGC
| (C | G), (A | T) -> -. gBGC
......
......@@ -6,7 +6,6 @@ module Evolution_model : sig
type context_CpG = {
lh: Nucleotide.repr option;
rh: Nucleotide.repr option;
codon: (Nucleotide.repr * Nucleotide.repr * Nucleotide.repr)
}
type param = {
......
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