-
Notifications
You must be signed in to change notification settings - Fork 3
Expand file tree
/
Copy pathREADME.html
More file actions
527 lines (450 loc) · 22.3 KB
/
Copy pathREADME.html
File metadata and controls
527 lines (450 loc) · 22.3 KB
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
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no">
<style>
body {
font-family: Helvetica, arial, sans-serif;
font-size: 14px;
line-height: 1.6;
padding-top: 10px;
padding-bottom: 10px;
background-color: white;
padding: 30px; }
body > *:first-child {
margin-top: 0 !important; }
body > *:last-child {
margin-bottom: 0 !important; }
a {
color: #4183C4; }
a.absent {
color: #cc0000; }
a.anchor {
display: block;
padding-left: 30px;
margin-left: -30px;
cursor: pointer;
position: absolute;
top: 0;
left: 0;
bottom: 0; }
h1, h2, h3, h4, h5, h6 {
margin: 20px 0 10px;
padding: 0;
font-weight: bold;
-webkit-font-smoothing: antialiased;
cursor: text;
position: relative; }
h1:hover a.anchor, h2:hover a.anchor, h3:hover a.anchor, h4:hover a.anchor, h5:hover a.anchor, h6:hover a.anchor {
background: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABAAAAAQCAYAAAAf8/9hAAAAGXRFWHRTb2Z0d2FyZQBBZG9iZSBJbWFnZVJlYWR5ccllPAAAA09pVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADw/eHBhY2tldCBiZWdpbj0i77u/IiBpZD0iVzVNME1wQ2VoaUh6cmVTek5UY3prYzlkIj8+IDx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IkFkb2JlIFhNUCBDb3JlIDUuMy1jMDExIDY2LjE0NTY2MSwgMjAxMi8wMi8wNi0xNDo1NjoyNyAgICAgICAgIj4gPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4gPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIgeG1sbnM6eG1wPSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvIiB4bWxuczp4bXBNTT0iaHR0cDovL25zLmFkb2JlLmNvbS94YXAvMS4wL21tLyIgeG1sbnM6c3RSZWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20veGFwLzEuMC9zVHlwZS9SZXNvdXJjZVJlZiMiIHhtcDpDcmVhdG9yVG9vbD0iQWRvYmUgUGhvdG9zaG9wIENTNiAoMTMuMCAyMDEyMDMwNS5tLjQxNSAyMDEyLzAzLzA1OjIxOjAwOjAwKSAgKE1hY2ludG9zaCkiIHhtcE1NOkluc3RhbmNlSUQ9InhtcC5paWQ6OUM2NjlDQjI4ODBGMTFFMTg1ODlEODNERDJBRjUwQTQiIHhtcE1NOkRvY3VtZW50SUQ9InhtcC5kaWQ6OUM2NjlDQjM4ODBGMTFFMTg1ODlEODNERDJBRjUwQTQiPiA8eG1wTU06RGVyaXZlZEZyb20gc3RSZWY6aW5zdGFuY2VJRD0ieG1wLmlpZDo5QzY2OUNCMDg4MEYxMUUxODU4OUQ4M0REMkFGNTBBNCIgc3RSZWY6ZG9jdW1lbnRJRD0ieG1wLmRpZDo5QzY2OUNCMTg4MEYxMUUxODU4OUQ4M0REMkFGNTBBNCIvPiA8L3JkZjpEZXNjcmlwdGlvbj4gPC9yZGY6UkRGPiA8L3g6eG1wbWV0YT4gPD94cGFja2V0IGVuZD0iciI/PsQhXeAAAABfSURBVHjaYvz//z8DJYCRUgMYQAbAMBQIAvEqkBQWXI6sHqwHiwG70TTBxGaiWwjCTGgOUgJiF1J8wMRAIUA34B4Q76HUBelAfJYSA0CuMIEaRP8wGIkGMA54bgQIMACAmkXJi0hKJQAAAABJRU5ErkJggg==) no-repeat 10px center;
text-decoration: none; }
h1 tt, h1 code {
font-size: inherit; }
h2 tt, h2 code {
font-size: inherit; }
h3 tt, h3 code {
font-size: inherit; }
h4 tt, h4 code {
font-size: inherit; }
h5 tt, h5 code {
font-size: inherit; }
h6 tt, h6 code {
font-size: inherit; }
h1 {
font-size: 28px;
color: black; }
h2 {
font-size: 24px;
border-bottom: 1px solid #cccccc;
color: black; }
h3 {
font-size: 18px; }
h4 {
font-size: 16px; }
h5 {
font-size: 14px; }
h6 {
color: #777777;
font-size: 14px; }
p, blockquote, ul, ol, dl, li, table, pre {
margin: 15px 0; }
hr {
background: transparent url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAAYAAAAECAYAAACtBE5DAAAAGXRFWHRTb2Z0d2FyZQBBZG9iZSBJbWFnZVJlYWR5ccllPAAAAyJpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADw/eHBhY2tldCBiZWdpbj0i77u/IiBpZD0iVzVNME1wQ2VoaUh6cmVTek5UY3prYzlkIj8+IDx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IkFkb2JlIFhNUCBDb3JlIDUuMC1jMDYwIDYxLjEzNDc3NywgMjAxMC8wMi8xMi0xNzozMjowMCAgICAgICAgIj4gPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4gPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIgeG1sbnM6eG1wPSJodHRwOi8vbnMuYWRvYmUuY29tL3hhcC8xLjAvIiB4bWxuczp4bXBNTT0iaHR0cDovL25zLmFkb2JlLmNvbS94YXAvMS4wL21tLyIgeG1sbnM6c3RSZWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20veGFwLzEuMC9zVHlwZS9SZXNvdXJjZVJlZiMiIHhtcDpDcmVhdG9yVG9vbD0iQWRvYmUgUGhvdG9zaG9wIENTNSBNYWNpbnRvc2giIHhtcE1NOkluc3RhbmNlSUQ9InhtcC5paWQ6OENDRjNBN0E2NTZBMTFFMEI3QjRBODM4NzJDMjlGNDgiIHhtcE1NOkRvY3VtZW50SUQ9InhtcC5kaWQ6OENDRjNBN0I2NTZBMTFFMEI3QjRBODM4NzJDMjlGNDgiPiA8eG1wTU06RGVyaXZlZEZyb20gc3RSZWY6aW5zdGFuY2VJRD0ieG1wLmlpZDo4Q0NGM0E3ODY1NkExMUUwQjdCNEE4Mzg3MkMyOUY0OCIgc3RSZWY6ZG9jdW1lbnRJRD0ieG1wLmRpZDo4Q0NGM0E3OTY1NkExMUUwQjdCNEE4Mzg3MkMyOUY0OCIvPiA8L3JkZjpEZXNjcmlwdGlvbj4gPC9yZGY6UkRGPiA8L3g6eG1wbWV0YT4gPD94cGFja2V0IGVuZD0iciI/PqqezsUAAAAfSURBVHjaYmRABcYwBiM2QSA4y4hNEKYDQxAEAAIMAHNGAzhkPOlYAAAAAElFTkSuQmCC) repeat-x 0 0;
border: 0 none;
color: #cccccc;
height: 4px;
padding: 0;
}
body > h2:first-child {
margin-top: 0;
padding-top: 0; }
body > h1:first-child {
margin-top: 0;
padding-top: 0; }
body > h1:first-child + h2 {
margin-top: 0;
padding-top: 0; }
body > h3:first-child, body > h4:first-child, body > h5:first-child, body > h6:first-child {
margin-top: 0;
padding-top: 0; }
a:first-child h1, a:first-child h2, a:first-child h3, a:first-child h4, a:first-child h5, a:first-child h6 {
margin-top: 0;
padding-top: 0; }
h1 p, h2 p, h3 p, h4 p, h5 p, h6 p {
margin-top: 0; }
li p.first {
display: inline-block; }
li {
margin: 0; }
ul, ol {
padding-left: 30px; }
ul :first-child, ol :first-child {
margin-top: 0; }
dl {
padding: 0; }
dl dt {
font-size: 14px;
font-weight: bold;
font-style: italic;
padding: 0;
margin: 15px 0 5px; }
dl dt:first-child {
padding: 0; }
dl dt > :first-child {
margin-top: 0; }
dl dt > :last-child {
margin-bottom: 0; }
dl dd {
margin: 0 0 15px;
padding: 0 15px; }
dl dd > :first-child {
margin-top: 0; }
dl dd > :last-child {
margin-bottom: 0; }
blockquote {
border-left: 4px solid #dddddd;
padding: 0 15px;
color: #777777; }
blockquote > :first-child {
margin-top: 0; }
blockquote > :last-child {
margin-bottom: 0; }
table {
padding: 0;border-collapse: collapse; }
table tr {
border-top: 1px solid #cccccc;
background-color: white;
margin: 0;
padding: 0; }
table tr:nth-child(2n) {
background-color: #f8f8f8; }
table tr th {
font-weight: bold;
border: 1px solid #cccccc;
text-align: left;
margin: 0;
padding: 6px 13px; }
table tr td {
border: 1px solid #cccccc;
text-align: left;
margin: 0;
padding: 6px 13px; }
table tr th :first-child, table tr td :first-child {
margin-top: 0; }
table tr th :last-child, table tr td :last-child {
margin-bottom: 0; }
img {
max-width: 100%; }
span.frame {
display: block;
overflow: hidden; }
span.frame > span {
border: 1px solid #dddddd;
display: block;
float: left;
overflow: hidden;
margin: 13px 0 0;
padding: 7px;
width: auto; }
span.frame span img {
display: block;
float: left; }
span.frame span span {
clear: both;
color: #333333;
display: block;
padding: 5px 0 0; }
span.align-center {
display: block;
overflow: hidden;
clear: both; }
span.align-center > span {
display: block;
overflow: hidden;
margin: 13px auto 0;
text-align: center; }
span.align-center span img {
margin: 0 auto;
text-align: center; }
span.align-right {
display: block;
overflow: hidden;
clear: both; }
span.align-right > span {
display: block;
overflow: hidden;
margin: 13px 0 0;
text-align: right; }
span.align-right span img {
margin: 0;
text-align: right; }
span.float-left {
display: block;
margin-right: 13px;
overflow: hidden;
float: left; }
span.float-left span {
margin: 13px 0 0; }
span.float-right {
display: block;
margin-left: 13px;
overflow: hidden;
float: right; }
span.float-right > span {
display: block;
overflow: hidden;
margin: 13px auto 0;
text-align: right; }
code, tt {
margin: 0 2px;
padding: 0 5px;
white-space: nowrap;
border: 1px solid #eaeaea;
background-color: #f8f8f8;
border-radius: 3px; }
pre code {
margin: 0;
padding: 0;
white-space: pre;
border: none;
background: transparent; }
.highlight pre {
background-color: #f8f8f8;
border: 1px solid #cccccc;
font-size: 13px;
line-height: 19px;
overflow: auto;
padding: 6px 10px;
border-radius: 3px; }
pre {
background-color: #f8f8f8;
border: 1px solid #cccccc;
font-size: 13px;
line-height: 19px;
overflow: auto;
padding: 6px 10px;
border-radius: 3px; }
pre code, pre tt {
background-color: transparent;
border: none; }
@media screen and (min-width: 914px) {
body {
width: 854px;
margin:0 auto;
}
}
@media print {
table, pre {
page-break-inside: avoid;
}
}
</style>
<title>TdG09 program: identifying changes in selective constraints</title>
</head>
<body>
<h1>TdG09 program: identifying changes in selective constraints</h1>
<p>This program is an implementation of the model described in:</p>
<p>Tamuri AU, dos Reis M, Hay AJ, Goldstein RA (2009) Identifying Changes in
Selective Constraints: Host Shifts in Influenza. <em>PLoS Comput Biol</em> 5(11): e1000564. <a href="http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1000564">doi:10.1371/journal.pcbi.1000564</a></p>
<p>This phylogenetic model uses site-specific amino acid frequencies to distinguish patterns of substitution between two (or more) lineages or groups of taxa. We first estimate the site-specific amino acid frequencies at a given location in a protein alignment assuming that the pattern of substitution is homogeneous across all branches. We then estimate multiple sets of site-specific frequencies, allowing them to differ among different branches, producing a non-homogeneous model of evolutionary change. Using statistical tests, we then see whether the non-homogeneous model provides a significantly better fit to the data than the homogenous model.</p>
<h2>Tutorial</h2>
<ol>
<li><p><strong>Install Java:</strong> The program requires a recent version of the Java Runtime (JRE), which, if not already installed, can be downloaded from <a href="http://www.oracle.com/technetwork/java/javase/downloads/index.html">Oracle</a>. Linux packages are usually available in the distribution's repository (e.g. <code>sudo apt-get install openjdk-7-jre</code> for Debian, Ubuntu etc. distributions).</p></li>
<li><p><strong>Download the program:</strong> The latest version of the program is available for download from the <a href="https://github.com/tamuri/tdg09/tags">repository</a>. The download includes a compiled binary as well as a <code>build.xml</code> to compile from sources using <a href="http://ant.apache.org/">ant</a>. Unzip the download and check that the program works:</p>
<pre><code> $ ls -F
README.md build.xml dist/ etc/ lib/ src/
$ java -cp dist/tdg09.jar tdg09.Analyse
Usage: java -cp tdg09.jar tdg09.Analyse [options]
Options:
* -alignment
Alignment in PHYLIP format
* -groups
Group labels to partition tree e.g. Av Hu
-threads
Number of threads to use
Default: 1
* -tree
Tree in NEWICK format
</code></pre></li>
<li><p><strong>Preparing your data:</strong> The program requires a protein sequence alignment in <a href="http://www.bioperl.org/wiki/PHYLIP_multiple_alignment_format">PHYLIP format</a> and the corresponding tree in <a href="http://en.wikipedia.org/wiki/Newick_format">Newick format</a>. The tree branch lengths should be optimised for an amino acid model, such as WAG, using a program such as <a href="http://sco.h-its.org/exelixis/software.html">RAxML</a> or <a href="http://abacus.gene.ucl.ac.uk/software/paml.html">PAML</a>. Most importantly, sequence names must be prefixed by a two-letter identifier that is used to indicate its lineage/grouping. For example, to identify changes in selective constraints between avian and human flu viral proteins, we can prefix every sequence name with 'Av' or 'Hu' to indicate its lineage:</p>
<pre><code> $ cat etc/H1.faa
434 566
Hu_HA_AAX56530_H1N1 MKAKLLVLLCAFTATYADTI...
Hu_HA_AAY78939_H1N2 MKVKLLILLCTFTATYADTI...
Av_HA_ABB19507_H1N6 MEAKLFVLFCTFTVLKADTI...
Av_HA_ABB19518_H1N1 MEAKLFVLFCTFTALKADTI...
...
</code></pre>
<p>Both the sequence alignment and the tree must follow this convention.</p>
<pre><code> $ cat etc/H1.tree
((Av_HA_ABB19607_H1N1:0.0852880,Av_HA_ABG88212_H1N1:0.1036700):
0.0248320,((Av_HA_ABB19618_H1N1:...
</code></pre>
<p>This is the only way by which the tdg09 program determines which of the different non-homogeneous models a particular tree branch should use. Example data sets of flu viral proteins (used in Tamuri <em>et al.</em> 2009) are included in the <code>etc/</code> directory.</p></li>
<li><p><strong>Running the program:</strong> The command-line options for the program are:</p>
<ul>
<li>-alignment : the sequence alignment file in PHYLIP format e.g. <code>etc/H1.faa</code></li>
<li>-tree : the tree file in Newick format e.g. <code>etc/H1.tree</code></li>
<li>-groups : the two-letter identifiers used to partition the sequences e.g. <code>Av Hu</code></li>
<li>-threads : specify the numbers of CPU cores/threads to utilise e.g. <code>2</code></li>
</ul>
<p> The program prints messages to standard out, so this should be captured using <code>></code> or <code>tee</code>. We are now ready to run the program:</p>
<pre><code> $ java -cp dist/tdg09.jar tdg09.Analyse -alignment etc/H1.faa \
-tree etc/H1.tree -groups Av Hu -threads 2 > H1_out.txt
</code></pre>
<p> or if you have 'tee' installed:</p>
<pre><code> $ java -cp dist/tdg09.jar tdg09.Analyse -alignment etc/H1.faa \
-tree etc/H1.tree -groups Av Hu -threads 2 | tee H1_out.txt
</code></pre></li>
<li><p><strong>Inspect the results:</strong> In this example, program output is captured in <code>H1_out.txt</code>:</p>
<pre><code> $ cat H1_out.txt
StartTime: 2013-03-19 13:44:10.129
WorkingDirectory: /Users/Tester/Documents/tdg09
Options: -alignment etc/H1.faa -tree etc/H1.tree -groups Av Hu -threads 2
TreeFile: /Users/Tester/Documents/tdg09/etc/H1.tree
AlignmentFile: /Users/Tester/Documents/tdg09/etc/H1.faa
Alignment:
SequenceCount: 434
SiteCount: 566
Groups: [Av, Hu]
# The internal nodes of the tree are not labelled. Labelling...
# Node 432 from [Av, Hu] resolved
# Node 432 from [Av, Hu] resolved
# Assuming that root of tree is in group [Av]
# Switching from group [Av] to [Hu] at branch 432..431
LabelledTree: >
(((((((((((Av_HA_ABB19607_H1N1:0.0852880,Av_HA_ABG88212_H1N1:0.1036700)
Av:0.0248320,((Av_HA_ABB19618_H1N1:0.0910450,Av_HA_ABG88201_H1N1:0.0810480)
...
Hu:0.0618410)Hu:0.0978830)Hu:0.0436020)Hu:0.0482480)Hu:0.0255390)Hu:0.0362960)
Hu:0.0696670)Hu:0.0620790)Hu:0.0268645);
</code></pre>
<p>The output contains a "LabelledTree" that shows the inferred lineage for each ancestral node. This should be checked in a tree viewing program (such as <a href="http://www-ab.informatik.uni-tuebingen.de/software/dendroscope">Dendroscope</a>) to make sure that the lineages are correct. If not, they can be modified and the analysis can be rerun with the new, custom-labelled, tree. The output continues:</p>
<pre><code> # 2013-03-19 22:44:10.33 - site 1 complete.
# 2013-03-19 22:44:11.366 - site 3 complete.
# 2013-03-19 22:44:11.366 - site 4 complete.
...
LrtResults:
# Site, delta lnL, dof, LRT, FDR
- [ 204, 21.850690, 3, 0.0000000, 0.0000003 ]
- [ 169, 11.528011, 1, 0.0000016, 0.0001542 ]
- [ 289, 10.733096, 2, 0.0000218, 0.0008550 ]
- [ 252, 10.927903, 2, 0.0000180, 0.0008796 ]
- [ 9, 8.225516, 1, 0.0000499, 0.0008895 ]
- [ 300, 15.027774, 5, 0.0000144, 0.0009396 ]
- [ 62, 8.261257, 1, 0.0000481, 0.0009423 ]
- [ 303, 8.261836, 1, 0.0000480, 0.0010463 ]
- [ 239, 9.647869, 2, 0.0000646, 0.0010545 ]
- [ 315, 9.483944, 2, 0.0000761, 0.0011468 ]
- [ 253, 8.262439, 1, 0.0000480, 0.0011764 ]
...
FullResults:
# Site, WAG+ssF params, WAG+ssF lnL, WAG+lssF params, WAG+lssF params, delta lnL, dof, LRT, FDR
- [ 1, NA, NA, NA, NA, NA, NA, NA, NA ]
- [ 2, 3, -21.488149, 5, -13.378828, 8.109321, 2, 0.0003007, 0.0025627 ]
- [ 3, 2, -48.561153, 3, -47.763333, 0.797820, 1, 0.2065223, 0.3489515 ]
- [ 4, 2, -24.137478, 3, -23.407924, 0.729554, 1, 0.2270721, 0.3708845 ]
- [ 5, 2, -12.611075, 3, -12.456865, 0.154209, 1, 0.5786522, 0.6593943 ]
- [ 6, 3, -23.856049, 5, -22.761440, 1.094609, 2, 0.3346706, 0.5084918 ]
- [ 7, 2, -44.839152, 3, -43.761797, 1.077355, 1, 0.1421333, 0.2509740 ]
- [ 8, NA, NA, NA, NA, NA, NA, NA, NA ]
...
</code></pre>
<p> Of interest are the "LrtResults" and "FullResults" tables.</p>
<p> The LrtResults lists polymorphic sites (on which the non-homogeneous model was estimated) and orders them by the false discovery rate (a correction on the likelihood ratio test P-value required by multiple hypothesis testing). At a given FDR cutoff (e.g. 0.05), these sites are those at which the non-homogeneous model provides a statistically significant improvement over the homogenous model, indicating that the patterns of substitution are different between the different groups/lineages.</p>
<p> The FullResults table lists further results from all sites. This includes the log-likelihood for the WAG+ssF (site-specific frequencies or <em>homogeneous model</em>) and WAG+lssF (lineage and site-specific frequencies or <em>non-homogeneous model</em>). Conserved locations are not analysed, so their entries are 'NA'. Finally, after the FullResults table, the output file contains details of the per-location results, such as the amino acid frequencies estimated by the ssF and lssF models.</p></li>
<li><p><strong>Analysing the results</strong>: The output file is in <a href="http://www.yaml.org/">YAML</a> format, which means it can be read by any other programming language that has a YAML parsing library. Here we show an example of analysing the results using the programming language <a href="http://www.r-project.org/">R</a>. The code is available in the <code>src/R/example.R</code> file. Start the R console, install and load the <em>yaml</em> library, then load the output file using the <code>yaml.load_file</code> function:</p>
<pre><code> $ R
R version 2.15.3 (2013-03-01) -- "Security Blanket"
...
> install.packages("yaml")
Installing package(s) into ‘/Users/Tester/Library/R/2.15/library’
(as ‘lib’ is unspecified)
...
> library(yaml)
> out <- yaml.load_file(input='/Users/Tester/Documents/tdg09/H1_out.txt')
> summary(out)
Length Class Mode
StartTime 1 -none- character
WorkingDirectory 1 -none- character
Options 1 -none- character
TreeFile 1 -none- character
AlignmentFile 1 -none- character
Alignment 2 -none- list
Groups 2 -none- character
LabelledTree 1 -none- character
LrtResults 196 -none- list
FullResults 566 -none- list
ConservedPositions 2 -none- list
SiteResults 566 -none- list
EndTime 1 -none- character
</code></pre>
<p> To convert the R list objects into a type easier to work with, we convert the LrtResults and FullResults objects into data.frames:</p>
<pre><code> > lrt_results <- as.data.frame(matrix(unlist(out$LrtResults), ncol=5, byrow=T))
> names(lrt_results) <- c("site", "deltaLnL", "dof", "lrt", "fdr")
> head(lrt_results)
site deltaLnL dof lrt fdr
1 204 21.850690 3 0.00e+00 0.0000003
2 169 11.528011 1 1.60e-06 0.0001542
3 289 10.733096 2 2.18e-05 0.0008550
4 252 10.927903 2 1.80e-05 0.0008796
5 9 8.225516 1 4.99e-05 0.0008895
6 300 15.027774 5 1.44e-05 0.0009396
> sum(lrt_results$fdr <= 0.05) # how many sites identified with FDR < 0.05?
[1] 55
> full_results <- as.data.frame(matrix(unlist(out$FullResults), ncol=9, byrow=T))
> names(full_results) <- c("site", "ssfParams", "ssfLnL", "lssfParams", "lssfLnL", "deltaLnL", "dof", "lrt", "fdr")
> head(full_results)
site ssfParams ssfLnL lssfParams lssfLnL deltaLnL dof lrt fdr
1 1 NA NA NA NA NA NA NA NA
2 2 3 -21.488149 5 -13.378828 8.109321 2 0.0003007 0.0025627
3 3 2 -48.561153 3 -47.763333 0.79782 1 0.2065223 0.3489515
4 4 2 -24.137478 3 -23.407924 0.729554 1 0.2270721 0.3708845
5 5 2 -12.611075 3 -12.456865 0.154209 1 0.5786522 0.6593943
6 6 3 -23.856049 5 -22.76144 1.094609 2 0.3346706 0.5084918
</code></pre>
<p> To produce a plot of FDR values by site:</p>
<pre><code> > fdr <- as.numeric(levels(full_results$fdr)[full_results$fdr]) # FDR column should be numeric
Warning message:
NAs introduced by coercion
> fdr[is.na(fdr)] <- 1.0 # conserved locations implicitly have no evidence of non-homogeneity
> sites <- out$Alignment$SiteCount
> plot_ranges <- split(seq(1, sites), cut(seq(1, sites), 5)) # split sites in plot into 5 rows
> par(mfrow=c(5,1), mar=c(2.0,0.5,0.5,0.5))
for (p in 1:5) {
plot(1 - fdr,
xlim=c(plot_ranges[[p]][1], tail(plot_ranges[[p]], n=1)),
ty='h', lwd=1, main='', xlab='', ylab='', yaxt='n', col="#1B9E77")
lines(which(fdr <= 0.20), 1 - fdr[fdr <= 0.20], col="#D95F02", ty='h')
abline(h=0.95, lty='dashed')
points(which(fdr <= 0.05), 1 - fdr[fdr <= 0.05], pch=20, col="#DE2D26")
}
</code></pre>
<p> This produces the plot show below. Bars drawn in orange indicate locations with FDR < 0.20, and those locations with FDR < 0.05 have a red dot at their value (we drew the plot with "1 - fdr" values so that smaller FDRs are taller). We can see from this plot a cluster of identified sites at locations 200-210.</p>
<p><img src="./figure_fdr_plot.png" alt="image" /></p>
<p> The R code in <code>src/R/example.R</code> also provides an example of using the <a href="http://cran.r-project.org/web/packages/phangorn/index.html">phangorn</a> and <a href="http://cran.r-project.org/web/packages/ape/">ape</a> packages to simulate data under the homogeneous model for a particular site. We can then analyse the synthetic data with the TdG09 software, the parametric bootstrap providing significance of the non-homogeneous model on the original data. For example, here is the Monte Carlo distribution of Δ for 1000 parametric bootstrap replicates for flu protein HA (H1) site 204, showing that the homogenous model can be rejected in favour of the non-homogeneous model with P-value < 0.001.</p>
<p><img src="./figure_cox_plot.png" alt="image" /></p></li>
</ol>
</body>
</html>