/
fels81.eg
76 lines (65 loc) · 2.37 KB
/
fels81.eg
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
;; The grammar.
(grammar
(name Felsenstein81)
;; Transformation rules. These follow the pattern for a null model with rate matrix X.
;; See e.g. the Kimura two-parameter xgram model for more info.
(transform (from (S)) (to (X S*)))
(transform (from (S*)) (to (S)) (prob 1))
(transform (from (S*)) (to ()) (prob 1))
(update-rules 0)
;; Parameters for the Felsenstein '81 model
(params
((piA .25) (piC .25) (piG .25) (piT .25)) ;; base frequencies
((R 1)) ;; substitution rate
) ;; end params
;; Rate matrix
(chain
;; The state of this chain is a single nucleotide X.
(terminal (X))
;; Treat initial probabilities and mutation rates as functions, not variables.
(update-policy parametric)
;; initial probability distribution
(initial (state (a)) (prob piA))
(initial (state (c)) (prob piC))
(initial (state (g)) (prob piG))
(initial (state (t)) (prob piT))
;; mutation rates
(mutate (from (a)) (to (c)) (rate piC * R))
(mutate (from (a)) (to (g)) (rate piG * R))
(mutate (from (a)) (to (t)) (rate piT * R))
(mutate (from (c)) (to (a)) (rate piA * R))
(mutate (from (c)) (to (g)) (rate piG * R))
(mutate (from (c)) (to (t)) (rate piT * R))
(mutate (from (g)) (to (a)) (rate piA * R))
(mutate (from (g)) (to (c)) (rate piC * R))
(mutate (from (g)) (to (t)) (rate piT * R))
(mutate (from (t)) (to (a)) (rate piA * R))
(mutate (from (t)) (to (c)) (rate piC * R))
(mutate (from (t)) (to (g)) (rate piG * R))
;; "unobserved" mutations that are needed to accumulate probabilistic counts correctly
(mutate (from (a)) (to (a)) (rate piA * R))
(mutate (from (c)) (to (c)) (rate piC * R))
(mutate (from (g)) (to (g)) (rate piG * R))
(mutate (from (t)) (to (t)) (rate piT * R))
) ;; end chain X
) ;; end grammar
;; Standard DNA alphabet
(alphabet
(name DNA)
(token (a c g t))
(complement (t g c a))
(extend (to n) (from a) (from c) (from g) (from t))
(extend (to x) (from a) (from c) (from g) (from t))
(extend (to u) (from t))
(extend (to r) (from a) (from g))
(extend (to y) (from c) (from t))
(extend (to m) (from a) (from c))
(extend (to k) (from g) (from t))
(extend (to s) (from c) (from g))
(extend (to w) (from a) (from t))
(extend (to h) (from a) (from c) (from t))
(extend (to b) (from c) (from g) (from t))
(extend (to v) (from a) (from c) (from g))
(extend (to d) (from a) (from g) (from t))
(wildcard *)
) ;; end alphabet DNA