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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
|
Allowing substitution model parameters to differ between branches
=================================================================
.. sectionauthor:: Gavin Huttley
A common task concerns assessing how substitution model exchangeability parameters differ between evolutionary lineages. This is most commonly of interest for the case of testing for natural selection. Here I'll demonstrate the different ways of scoping parameters across trees for the codon model case and how these can be used for evolutionary modelling.
We start with the standard imports, plus using a canned codon substitution model and then load the sample data set.
.. doctest::
>>> from cogent import LoadSeqs, LoadTree
>>> from cogent.evolve.models import MG94HKY
>>> aln = LoadSeqs("data/long_testseqs.fasta")
>>> tree = LoadTree("data/test.tree")
We construct the substitution model and likelihood function and set the alignment.
.. doctest::
>>> sm = MG94HKY()
>>> lf = sm.makeLikelihoodFunction(tree, digits=2, space=3)
>>> lf.setAlignment(aln)
At this point we have a likelihood function with two exchangeability parameters from the substitution model (``kappa`` the transition/transversion ratio; ``omega`` the nonsynonymous/synonymous ratio) plus branch lengths for all tree edges. To facilitate subsequent discussion I now display the tree
.. doctest::
>>> print tree.asciiArt()
/-Human
/edge.0--|
/edge.1--| \-HowlerMon
| |
| \-Mouse
-root----|
|--NineBande
|
\-DogFaced
In order to scope a parameter on a tree (meaning specifying a subset of edges for which the parameter is to be treated differently to the remainder of the tree) requires uniquely identifying the edges. We do this using the following arguments to the likelihood function ``setParamRule`` method:
- ``tip_names``: the name of two tips
- ``outgroup_name``: the name of a tip that is not part of the clade of interest
- ``is_clade``: if ``True``, all lineages descended from the tree node identified by the ``tip_names`` and ``outgroup_name`` argument are affected by the other arguments. If ``False``, then the ``is_stem`` argument must apply.
- ``is_stem``: Whether the edge leading to the node is included.
The next concepts include exactly what can be scoped and how. In the case of testing for distinctive periods of natural selection it is common to specify distinct values for ``omega`` for an edge. I'll first illustrate some possible uses for the arguments above by setting ``omega`` to be distinctive for specific edges. I will set a value for omega so that printing the likelihood function illustrates what edges have been effected, but I won't actually do any model fitting.
Specifying a clade
------------------
I'm going to cause ``omega`` to attain a different value for all branches aside from the primate clade and stem (``HowlerMon``, ``Human``, ``edge.0``).
.. doctest::
>>> lf.setParamRule('omega', tip_names=['DogFaced', 'Mouse'],
... outgroup_name='Human', init=2.0, is_clade=True)
>>> print lf
Likelihood Function Table
=====
kappa
-----
1.00
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.03 1.00
HowlerMon edge.0 0.04 1.00
edge.0 edge.1 0.04 1.00
Mouse edge.1 0.28 2.00
edge.1 root 0.02 2.00
NineBande root 0.09 2.00
DogFaced root 0.11 2.00
-----------------------------------
==============
motif mprobs
--------------
T 0.23
C 0.19
A 0.37
G 0.21
--------------
As you can see ``omega`` for the primate edges I listed above have the default parameter value (1.0), while the others have what I've assigned. In fact, you could omit the ``is_clade`` argument as this is the default, but I think for readability of scripts it's best to be explicit.
Specifying a stem
-----------------
This time I'll specify the stem leading to the primates as the edge of interest.
.. note:: I need to reset the ``lf`` so all edges have the default value again. I'll show this only for this example, but rest assured I'm doing it for all others too.
.. doctest::
>>> lf.setParamRule('omega', init=1.0)
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', init=2.0, is_stem=True, is_clade=False)
>>> print lf
Likelihood Function Table
=====
kappa
-----
1.00
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.03 1.00
HowlerMon edge.0 0.04 1.00
edge.0 edge.1 0.04 2.00
Mouse edge.1 0.28 1.00
edge.1 root 0.02 1.00
NineBande root 0.09 1.00
DogFaced root 0.11 1.00
-----------------------------------...
Specifying clade and stem
-------------------------
I'll specify that both the primates and their stem are to be considered.
.. doctest::
:hide:
>>> lf.setParamRule('omega', init=1.0)
.. doctest::
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', init=2.0, is_stem=True, is_clade=True)
>>> print lf
Likelihood Function Table
=====
kappa
-----
1.00
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.03 2.00
HowlerMon edge.0 0.04 2.00
edge.0 edge.1 0.04 2.00
Mouse edge.1 0.28 1.00
edge.1 root 0.02 1.00
NineBande root 0.09 1.00
DogFaced root 0.11 1.00
-----------------------------------...
Alternate arguments for specifying edges
----------------------------------------
The likelihood function ``setParamRule`` method also has the arguments of ``edge`` and ``edges``. These allow specific naming of the tree edge(s) to be affected by a rule. In general, however, the ``tip_names`` + ``outgroup_name`` combo is more robust.
Applications of scoped parameters
---------------------------------
The general use-cases for which a tree scope can be applied are:
1. constraining all edges identified by a rule to have a specific value which is constant and not modifiable
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True, is_const=True)
2. all edges identified by a rule have the same but different value to the rest of the tree
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True)
3. allowing all edges identified by a rule to have different values of the parameter with the remaining tree edges having the same value
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True, is_independent=True)
4. allowing all edges to have a different value
>>> lf.setParamRule('omega', is_independent=True)
I'll demonstrate these cases sequentially as they involve gradually increasing the degrees of freedom in the model. First we'll constrain ``omega`` to equal 1 on the primate edges. I'll then optimise the model.
.. note:: here I'm specifying a constant value for the parameter and so I **must** use the argument ``value`` to set it. This not to be confused with the argument ``init`` that is used for providing initial (starting) values for fitting.
.. doctest::
:hide:
>>> lf.setParamRule('omega', init=1.0)
.. doctest::
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True, value=1.0, is_const=True)
>>> lf.optimise(local=True, show_progress=False)
>>> print lf
Likelihood Function Table
=====
kappa
-----
3.87
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.09 1.00
HowlerMon edge.0 0.12 1.00
edge.0 edge.1 0.12 0.92
Mouse edge.1 0.84 0.92
edge.1 root 0.06 0.92
NineBande root 0.28 0.92
DogFaced root 0.34 0.92
-----------------------------------
==============
motif mprobs
--------------
T 0.23
C 0.19
A 0.37
G 0.21
--------------
>>> print lf.getLogLikelihood()
-8640.9...
>>> print lf.getNumFreeParams()
9
I'll now free up ``omega`` on the primate clade, but making it a single value shared by all primate lineages.
.. doctest::
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True, is_const=False)
>>> lf.optimise(local=True, show_progress=False)
>>> print lf
Likelihood Function Table
=====
kappa
-----
3.85
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.09 0.77
HowlerMon edge.0 0.12 0.77
edge.0 edge.1 0.12 0.92
Mouse edge.1 0.84 0.92
edge.1 root 0.06 0.92
NineBande root 0.28 0.92
DogFaced root 0.34 0.92
-----------------------------------
==============
motif mprobs
--------------
T 0.23
C 0.19
A 0.37
G 0.21
--------------
>>> print lf.getLogLikelihood()
-8639.7...
>>> print lf.getNumFreeParams()
10
Finally I'll allow all primate edges to have different values of ``omega``.
.. doctest::
>>> lf.setParamRule('omega', tip_names=['Human', 'HowlerMon'],
... outgroup_name='Mouse', is_clade=True, is_independent=True)
>>> lf.optimise(local=True, show_progress=False)
>>> print lf
Likelihood Function Table
=====
kappa
-----
3.85
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.09 0.59
HowlerMon edge.0 0.12 0.95
edge.0 edge.1 0.12 0.92
Mouse edge.1 0.84 0.92
edge.1 root 0.06 0.92
NineBande root 0.28 0.92
DogFaced root 0.34 0.92
-----------------------------------
==============
motif mprobs
--------------
T 0.23
C 0.19
A 0.37
G 0.21
--------------
>>> print lf.getLogLikelihood()
-8638.9...
>>> print lf.getNumFreeParams()
11
We now allow ``omega`` to be different on all edges.
.. doctest::
>>> lf.setParamRule('omega', is_independent=True)
>>> lf.optimise(local=True, show_progress=False)
>>> print lf
Likelihood Function Table
=====
kappa
-----
3.85
-----
===================================
edge parent length omega
-----------------------------------
Human edge.0 0.09 0.59
HowlerMon edge.0 0.12 0.95
edge.0 edge.1 0.12 1.13
Mouse edge.1 0.84 0.92
edge.1 root 0.06 0.38
NineBande root 0.28 1.27
DogFaced root 0.34 0.84
-----------------------------------
==============
motif mprobs
--------------
T 0.23
C 0.19
A 0.37
G 0.21
--------------
>>> print lf.getLogLikelihood()
-8636.1...
>>> print lf.getNumFreeParams()
15
|
