| line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
|
1
|
|
|
|
|
|
|
# |
|
2
|
|
|
|
|
|
|
# BioPerl module for Bio::Tools::Run::Phylo::SLR |
|
3
|
|
|
|
|
|
|
# |
|
4
|
|
|
|
|
|
|
# Please direct questions and support issues to |
|
5
|
|
|
|
|
|
|
# |
|
6
|
|
|
|
|
|
|
# Cared for by Albert Vilella |
|
7
|
|
|
|
|
|
|
# |
|
8
|
|
|
|
|
|
|
# Copyright Albert Vilella |
|
9
|
|
|
|
|
|
|
# |
|
10
|
|
|
|
|
|
|
# You may distribute this module under the same terms as perl itself |
|
11
|
|
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
# POD documentation - main docs before the code |
|
13
|
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
=head1 NAME |
|
15
|
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
Bio::Tools::Run::Phylo::SLR - Wrapper around the SLR program |
|
17
|
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
=head1 SYNOPSIS |
|
19
|
|
|
|
|
|
|
|
|
20
|
|
|
|
|
|
|
use Bio::Tools::Run::Phylo::SLR; |
|
21
|
|
|
|
|
|
|
use Bio::AlignIO; |
|
22
|
|
|
|
|
|
|
use Bio::TreeIO; |
|
23
|
|
|
|
|
|
|
use Bio::SimpleAlign; |
|
24
|
|
|
|
|
|
|
|
|
25
|
|
|
|
|
|
|
my $alignio = Bio::AlignIO->new |
|
26
|
|
|
|
|
|
|
(-format => 'fasta', |
|
27
|
|
|
|
|
|
|
-file => 't/data/219877.cdna.fasta'); |
|
28
|
|
|
|
|
|
|
|
|
29
|
|
|
|
|
|
|
my $aln = $alignio->next_aln; |
|
30
|
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
my $treeio = Bio::TreeIO->new |
|
32
|
|
|
|
|
|
|
(-format => 'newick', -file => 't/data/219877.tree'); |
|
33
|
|
|
|
|
|
|
|
|
34
|
|
|
|
|
|
|
my $tree = $treeio->next_tree; |
|
35
|
|
|
|
|
|
|
|
|
36
|
|
|
|
|
|
|
my $slr = Bio::Tools::Run::Phylo::SLR->new(); |
|
37
|
|
|
|
|
|
|
$slr->alignment($aln); |
|
38
|
|
|
|
|
|
|
$slr->tree($tree); |
|
39
|
|
|
|
|
|
|
# $rc = 1 for success, 0 for errors |
|
40
|
|
|
|
|
|
|
my ($rc,$results) = $slr->run(); |
|
41
|
|
|
|
|
|
|
|
|
42
|
|
|
|
|
|
|
my $positive_sites = $results->{'positive'}; |
|
43
|
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
print "# Site\tNeutral\tOptimal\tOmega\t", |
|
45
|
|
|
|
|
|
|
"lower\tupper\tLRT_Stat\tPval\tAdj.Pval\tResult\tNote\n"; |
|
46
|
|
|
|
|
|
|
foreach my $positive_site (@$positive_sites) { |
|
47
|
|
|
|
|
|
|
print |
|
48
|
|
|
|
|
|
|
$positive_site->[0], "\t", |
|
49
|
|
|
|
|
|
|
$positive_site->[1], "\t", |
|
50
|
|
|
|
|
|
|
$positive_site->[2], "\t", |
|
51
|
|
|
|
|
|
|
$positive_site->[3], "\t", |
|
52
|
|
|
|
|
|
|
$positive_site->[4], "\t", |
|
53
|
|
|
|
|
|
|
$positive_site->[5], "\t", |
|
54
|
|
|
|
|
|
|
$positive_site->[6], "\t", |
|
55
|
|
|
|
|
|
|
$positive_site->[7], "\t", |
|
56
|
|
|
|
|
|
|
$positive_site->[8], "\t", |
|
57
|
|
|
|
|
|
|
"positive\n"; |
|
58
|
|
|
|
|
|
|
} |
|
59
|
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
=head1 DESCRIPTION |
|
61
|
|
|
|
|
|
|
|
|
62
|
|
|
|
|
|
|
This is a wrapper around the SLR program. See |
|
63
|
|
|
|
|
|
|
http://www.ebi.ac.uk/goldman/SLR/ for more information. |
|
64
|
|
|
|
|
|
|
|
|
65
|
|
|
|
|
|
|
This module is more about generating the proper ctl file and |
|
66
|
|
|
|
|
|
|
will run the program in a separate temporary directory to avoid |
|
67
|
|
|
|
|
|
|
creating temp files all over the place. |
|
68
|
|
|
|
|
|
|
|
|
69
|
|
|
|
|
|
|
=head1 FEEDBACK |
|
70
|
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
=head2 Mailing Lists |
|
72
|
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
User feedback is an integral part of the evolution of this and other |
|
74
|
|
|
|
|
|
|
Bioperl modules. Send your comments and suggestions preferably to |
|
75
|
|
|
|
|
|
|
the Bioperl mailing list. Your participation is much appreciated. |
|
76
|
|
|
|
|
|
|
|
|
77
|
|
|
|
|
|
|
bioperl-l@bioperl.org - General discussion |
|
78
|
|
|
|
|
|
|
http://bioperl.org/wiki/Mailing_lists - About the mailing lists |
|
79
|
|
|
|
|
|
|
|
|
80
|
|
|
|
|
|
|
=head2 Support |
|
81
|
|
|
|
|
|
|
|
|
82
|
|
|
|
|
|
|
Please direct usage questions or support issues to the mailing list: |
|
83
|
|
|
|
|
|
|
|
|
84
|
|
|
|
|
|
|
I |
|
85
|
|
|
|
|
|
|
|
|
86
|
|
|
|
|
|
|
rather than to the module maintainer directly. Many experienced and |
|
87
|
|
|
|
|
|
|
reponsive experts will be able look at the problem and quickly |
|
88
|
|
|
|
|
|
|
address it. Please include a thorough description of the problem |
|
89
|
|
|
|
|
|
|
with code and data examples if at all possible. |
|
90
|
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
=head2 Reporting Bugs |
|
92
|
|
|
|
|
|
|
|
|
93
|
|
|
|
|
|
|
Report bugs to the Bioperl bug tracking system to help us keep track |
|
94
|
|
|
|
|
|
|
of the bugs and their resolution. Bug reports can be submitted via the |
|
95
|
|
|
|
|
|
|
web: |
|
96
|
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
http://redmine.open-bio.org/projects/bioperl/ |
|
98
|
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
=head1 AUTHOR - Albert Vilella |
|
100
|
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
Email avilella-at-gmail-dot-com |
|
102
|
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
=head1 CONTRIBUTORS |
|
104
|
|
|
|
|
|
|
|
|
105
|
|
|
|
|
|
|
Additional contributors names and emails here |
|
106
|
|
|
|
|
|
|
|
|
107
|
|
|
|
|
|
|
=head1 APPENDIX |
|
108
|
|
|
|
|
|
|
|
|
109
|
|
|
|
|
|
|
The rest of the documentation details each of the object methods. |
|
110
|
|
|
|
|
|
|
Internal methods are usually preceded with a _ |
|
111
|
|
|
|
|
|
|
|
|
112
|
|
|
|
|
|
|
=cut |
|
113
|
|
|
|
|
|
|
|
|
114
|
|
|
|
|
|
|
|
|
115
|
|
|
|
|
|
|
#' keep my emacs happy |
|
116
|
|
|
|
|
|
|
# Let the code begin... |
|
117
|
|
|
|
|
|
|
|
|
118
|
|
|
|
|
|
|
|
|
119
|
|
|
|
|
|
|
package Bio::Tools::Run::Phylo::SLR; |
|
120
|
1
|
|
|
1
|
|
136585
|
use vars qw(@ISA %VALIDVALUES $MINNAMELEN $PROGRAMNAME $PROGRAM); |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
62
|
|
|
121
|
1
|
|
|
1
|
|
5
|
use strict; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
16
|
|
|
122
|
1
|
|
|
1
|
|
4
|
use Bio::Root::Root; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
18
|
|
|
123
|
1
|
|
|
1
|
|
454
|
use Bio::AlignIO; |
|
|
1
|
|
|
|
|
71776
|
|
|
|
1
|
|
|
|
|
29
|
|
|
124
|
1
|
|
|
1
|
|
370
|
use Bio::TreeIO; |
|
|
1
|
|
|
|
|
13890
|
|
|
|
1
|
|
|
|
|
23
|
|
|
125
|
1
|
|
|
1
|
|
6
|
use Bio::SimpleAlign; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
16
|
|
|
126
|
1
|
|
|
1
|
|
395
|
use Bio::Tools::Run::WrapperBase; |
|
|
1
|
|
|
|
|
2
|
|
|
|
1
|
|
|
|
|
21
|
|
|
127
|
1
|
|
|
1
|
|
4
|
use Cwd; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
71
|
|
|
128
|
1
|
|
|
1
|
|
6
|
use File::Spec; |
|
|
1
|
|
|
|
|
2
|
|
|
|
1
|
|
|
|
|
209
|
|
|
129
|
|
|
|
|
|
|
|
|
130
|
|
|
|
|
|
|
@ISA = qw(Bio::Root::Root Bio::Tools::Run::WrapperBase); |
|
131
|
|
|
|
|
|
|
|
|
132
|
|
|
|
|
|
|
=head2 Default Values |
|
133
|
|
|
|
|
|
|
|
|
134
|
|
|
|
|
|
|
INCOMPLETE DOCUMENTATION OF ALL METHODS |
|
135
|
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
seqfile [incodon] |
|
137
|
|
|
|
|
|
|
File from which to read alignment of codon sequences. The file |
|
138
|
|
|
|
|
|
|
should be in PAML format. |
|
139
|
|
|
|
|
|
|
|
|
140
|
|
|
|
|
|
|
treefile [intree] |
|
141
|
|
|
|
|
|
|
File from which tree should be read. The tree should be in Nexus |
|
142
|
|
|
|
|
|
|
format |
|
143
|
|
|
|
|
|
|
|
|
144
|
|
|
|
|
|
|
outfile [slr.res] |
|
145
|
|
|
|
|
|
|
File to which results are written. If the file already exists, it will |
|
146
|
|
|
|
|
|
|
be overwritten. |
|
147
|
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
reoptimise [1] |
|
149
|
|
|
|
|
|
|
Should the branch lengths, omega and kappa be reoptimized? |
|
150
|
|
|
|
|
|
|
0 - no |
|
151
|
|
|
|
|
|
|
1 - yes. |
|
152
|
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
kappa [2.0] |
|
154
|
|
|
|
|
|
|
Value for kappa. If 'reoptimise' is specified, the value |
|
155
|
|
|
|
|
|
|
given will be used as am initial estimate, |
|
156
|
|
|
|
|
|
|
|
|
157
|
|
|
|
|
|
|
omega [0.1] |
|
158
|
|
|
|
|
|
|
Value for omega (dN/dS). If 'reoptimise' is specified, the value |
|
159
|
|
|
|
|
|
|
given will be used as an initial estimate. |
|
160
|
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
codonf [0] |
|
162
|
|
|
|
|
|
|
How codon frequencies are estimated: |
|
163
|
|
|
|
|
|
|
0: F61/F60 Estimates used are the empirical frequencies from the |
|
164
|
|
|
|
|
|
|
data. |
|
165
|
|
|
|
|
|
|
1: F3x4 The frequencies of nucleotides at each codon position |
|
166
|
|
|
|
|
|
|
are estimated from the data and then multiplied together to get the |
|
167
|
|
|
|
|
|
|
frequency of observing a given codon. The frequency of stop codons is |
|
168
|
|
|
|
|
|
|
set to zero, and all other frequencies scaled appropriately. |
|
169
|
|
|
|
|
|
|
2: F1x4 Nucleotide frequencies are estimated from the data |
|
170
|
|
|
|
|
|
|
(not taking into account at which position in the codon it occurs). |
|
171
|
|
|
|
|
|
|
The nucleotide frequencies are multiplied together to get the frequency |
|
172
|
|
|
|
|
|
|
of observing and then corrected for stop codons. |
|
173
|
|
|
|
|
|
|
|
|
174
|
|
|
|
|
|
|
freqtype [0] |
|
175
|
|
|
|
|
|
|
How codon frequencies are incorporated into the substitution matrix. |
|
176
|
|
|
|
|
|
|
0: q_{ij} = pi_{j} s_{ij} |
|
177
|
|
|
|
|
|
|
1: q_{ij} = \sqrt(pi_j/pi_i) s_{ij} |
|
178
|
|
|
|
|
|
|
2: q_{ij} = \pi_{n} s_{ij}, where n is the nucleotide that the |
|
179
|
|
|
|
|
|
|
subsitution is to. |
|
180
|
|
|
|
|
|
|
3: q_{ij} = s_{ij} / pi_i |
|
181
|
|
|
|
|
|
|
Option 0 is the tradition method of incorporating equilibrium frequencies |
|
182
|
|
|
|
|
|
|
into subsitution matrices (Felsenstein 1981; Goldman and Yang, 1994) |
|
183
|
|
|
|
|
|
|
Option 1 is described by Goldman and Whelan (2002), in this case with the |
|
184
|
|
|
|
|
|
|
additional parameter set to 0.5. |
|
185
|
|
|
|
|
|
|
Option 2 was suggested by Muse and Gaut (1994). |
|
186
|
|
|
|
|
|
|
Option 3 is included as an experiment, originally suggested by Bret Larget. |
|
187
|
|
|
|
|
|
|
it does not appear to describe evolution very successfully and should not |
|
188
|
|
|
|
|
|
|
be used for analyses. |
|
189
|
|
|
|
|
|
|
|
|
190
|
|
|
|
|
|
|
Kosakovsky-Pond has repeatedly stated that he finds incorporating codon |
|
191
|
|
|
|
|
|
|
frequencies in the manner of option 2 to be superior to option 0. We find |
|
192
|
|
|
|
|
|
|
that option 1 tends to perform better than either of these options. |
|
193
|
|
|
|
|
|
|
|
|
194
|
|
|
|
|
|
|
positive_only [0] |
|
195
|
|
|
|
|
|
|
If only positively selected sites are of interest, set this to "1". |
|
196
|
|
|
|
|
|
|
Calculation will be slightly faster, but information about sites under |
|
197
|
|
|
|
|
|
|
purifying selection is lost. |
|
198
|
|
|
|
|
|
|
|
|
199
|
|
|
|
|
|
|
gencode [universal] |
|
200
|
|
|
|
|
|
|
Which genetic code to use when determining whether a given mutation |
|
201
|
|
|
|
|
|
|
is synonymous or nonsynonymous. Currently only "universal" and |
|
202
|
|
|
|
|
|
|
"mammalian" mitochondrial are supported. |
|
203
|
|
|
|
|
|
|
|
|
204
|
|
|
|
|
|
|
nucleof [0] |
|
205
|
|
|
|
|
|
|
Allow for empirical exchangabilities for nucleotide substitution. |
|
206
|
|
|
|
|
|
|
0: No adjustment. All nucleotides treated the same, modulo |
|
207
|
|
|
|
|
|
|
transition / transversion. |
|
208
|
|
|
|
|
|
|
1: The rate at which a substitution caused a mutation from nucleotide |
|
209
|
|
|
|
|
|
|
a to nucleotide b is adjust by a constant N_{ab}. This adjustment is |
|
210
|
|
|
|
|
|
|
in addition to other adjustments (e.g. transition / transversion or |
|
211
|
|
|
|
|
|
|
base frequencies). |
|
212
|
|
|
|
|
|
|
|
|
213
|
|
|
|
|
|
|
aminof [0] |
|
214
|
|
|
|
|
|
|
Incorporate amino acid similarity parameters into substitution matrix, |
|
215
|
|
|
|
|
|
|
adjusting omega for a change between amino acid i and amino acid j. |
|
216
|
|
|
|
|
|
|
A_{ij} is a symmetric matrix of constants representing amino acid |
|
217
|
|
|
|
|
|
|
similarities. |
|
218
|
|
|
|
|
|
|
0: Constant omega for all amino acid changes |
|
219
|
|
|
|
|
|
|
1: omega_{ij} = omega^{A_{ij}} |
|
220
|
|
|
|
|
|
|
2: omega_{ij} = a_{ij} log(omega) / [ 1 - exp(-a_{ij} log(omega)) ] |
|
221
|
|
|
|
|
|
|
Option 1 has the same form as the original codon subsitution model |
|
222
|
|
|
|
|
|
|
proposed by Goldman and Yang (but with potentially different |
|
223
|
|
|
|
|
|
|
constants). |
|
224
|
|
|
|
|
|
|
Option 2 has a more population genetic derivtion, with omega being |
|
225
|
|
|
|
|
|
|
interpreted as the ratio of fixation probabilities. |
|
226
|
|
|
|
|
|
|
|
|
227
|
|
|
|
|
|
|
nucfile [nuc.dat] |
|
228
|
|
|
|
|
|
|
If nucleof is non-zero, read nucleotide substitution constants from |
|
229
|
|
|
|
|
|
|
nucfile. If this file does not exist, hard coded constants are used. |
|
230
|
|
|
|
|
|
|
|
|
231
|
|
|
|
|
|
|
aminofile [amino.dat] |
|
232
|
|
|
|
|
|
|
If aminof is non-zero, read amino acid similarity constants from |
|
233
|
|
|
|
|
|
|
aminofile. If this file does not exist, hard coded constants are used. |
|
234
|
|
|
|
|
|
|
|
|
235
|
|
|
|
|
|
|
timemem [0] |
|
236
|
|
|
|
|
|
|
Print summary of real time and CPU time used. Will eventually print |
|
237
|
|
|
|
|
|
|
summary of memory use as well. |
|
238
|
|
|
|
|
|
|
|
|
239
|
|
|
|
|
|
|
ldiff [3.841459] |
|
240
|
|
|
|
|
|
|
Twice log-likelihood difference used as a threshold for calculating |
|
241
|
|
|
|
|
|
|
support (confidence) intervals for sitewise omega estimates. This |
|
242
|
|
|
|
|
|
|
value should be the quantile from a chi-square distribution with one |
|
243
|
|
|
|
|
|
|
degree of freedom corresponding to the support required. |
|
244
|
|
|
|
|
|
|
E.g. qchisq(0.95,1) = 3.841459 |
|
245
|
|
|
|
|
|
|
0.4549364 = 50% support |
|
246
|
|
|
|
|
|
|
1.323304 = 75% support |
|
247
|
|
|
|
|
|
|
2.705543 = 90% support |
|
248
|
|
|
|
|
|
|
3.841459 = 95% support |
|
249
|
|
|
|
|
|
|
6.634897 = 99% support |
|
250
|
|
|
|
|
|
|
7.879439 = 99.5% support |
|
251
|
|
|
|
|
|
|
10.82757 = 99.9% support |
|
252
|
|
|
|
|
|
|
|
|
253
|
|
|
|
|
|
|
paramin [] |
|
254
|
|
|
|
|
|
|
If not blank, read in parameters from file given by the argument. |
|
255
|
|
|
|
|
|
|
|
|
256
|
|
|
|
|
|
|
paramout [] |
|
257
|
|
|
|
|
|
|
If not blank, write out parameter estimates to file given. |
|
258
|
|
|
|
|
|
|
|
|
259
|
|
|
|
|
|
|
skipsitewise [0] |
|
260
|
|
|
|
|
|
|
Skip sitewise estimation of omega. Depending on other options given, |
|
261
|
|
|
|
|
|
|
either calculate maximum likelihood or likelihood fixed at parameter |
|
262
|
|
|
|
|
|
|
values given. |
|
263
|
|
|
|
|
|
|
|
|
264
|
|
|
|
|
|
|
seed [0] |
|
265
|
|
|
|
|
|
|
Seed for random number generator. If seed is 0, then previously |
|
266
|
|
|
|
|
|
|
produced seed file (~/.rng64) is used. If this does not exist, the |
|
267
|
|
|
|
|
|
|
random number generator is initialised using the clock. |
|
268
|
|
|
|
|
|
|
|
|
269
|
|
|
|
|
|
|
saveseed [1] |
|
270
|
|
|
|
|
|
|
If non-zero, save finial seed in file (~/.rng64) to be used as initial |
|
271
|
|
|
|
|
|
|
seed in future runs of program. |
|
272
|
|
|
|
|
|
|
|
|
273
|
|
|
|
|
|
|
=head2 Results Format |
|
274
|
|
|
|
|
|
|
|
|
275
|
|
|
|
|
|
|
Results file (default: slr.res) |
|
276
|
|
|
|
|
|
|
------------ |
|
277
|
|
|
|
|
|
|
Results are presented in nine columns |
|
278
|
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
Site |
|
280
|
|
|
|
|
|
|
Number of sites in alignment |
|
281
|
|
|
|
|
|
|
|
|
282
|
|
|
|
|
|
|
Neutral |
|
283
|
|
|
|
|
|
|
(minus) Log-probability of observing site given that it was |
|
284
|
|
|
|
|
|
|
evolving neutrally (omega=1) |
|
285
|
|
|
|
|
|
|
|
|
286
|
|
|
|
|
|
|
Optimal |
|
287
|
|
|
|
|
|
|
(minus) Log-probability of observing site given that it was |
|
288
|
|
|
|
|
|
|
evolving at the optimal value of omega. |
|
289
|
|
|
|
|
|
|
|
|
290
|
|
|
|
|
|
|
Omega |
|
291
|
|
|
|
|
|
|
The value of omega which maximizes the log-probability of observing |
|
292
|
|
|
|
|
|
|
|
|
293
|
|
|
|
|
|
|
LRT_Stat |
|
294
|
|
|
|
|
|
|
Log-likelihood ratio statistic for non-neutral selection (or |
|
295
|
|
|
|
|
|
|
positive selection if the positive_only option is set to 1). |
|
296
|
|
|
|
|
|
|
LRT_Stat = 2 * (Neutral-Optimal) |
|
297
|
|
|
|
|
|
|
|
|
298
|
|
|
|
|
|
|
Pval |
|
299
|
|
|
|
|
|
|
P-value for non-neutral (or positive) selection at a site, |
|
300
|
|
|
|
|
|
|
unadjusted for multiple comparisons. |
|
301
|
|
|
|
|
|
|
|
|
302
|
|
|
|
|
|
|
Adj. Pval |
|
303
|
|
|
|
|
|
|
P-value for non-neutral (or positive) selection at a site, after |
|
304
|
|
|
|
|
|
|
adjusting for multiple comparisons using the Hochberg procedure |
|
305
|
|
|
|
|
|
|
(see the file "MultipleComparisons.txt" in the doc directory). |
|
306
|
|
|
|
|
|
|
|
|
307
|
|
|
|
|
|
|
Result |
|
308
|
|
|
|
|
|
|
A simple visual guide to the result. Sites detected as having been |
|
309
|
|
|
|
|
|
|
under positive selection are marked with a '+', sites under |
|
310
|
|
|
|
|
|
|
purifying selection are marked with '-'. The number of symbols |
|
311
|
|
|
|
|
|
|
Number symbols Threshold |
|
312
|
|
|
|
|
|
|
1 95% |
|
313
|
|
|
|
|
|
|
2 99% |
|
314
|
|
|
|
|
|
|
3 95% after adjustment |
|
315
|
|
|
|
|
|
|
4 99% after adjustment |
|
316
|
|
|
|
|
|
|
|
|
317
|
|
|
|
|
|
|
Occasionally the result may also contain an exclamation mark. This |
|
318
|
|
|
|
|
|
|
indicates that the observation at a site is not significantly |
|
319
|
|
|
|
|
|
|
different from random (equivalent to infinitely strong positive |
|
320
|
|
|
|
|
|
|
selection). This may indicate that the alignment at that site is bad |
|
321
|
|
|
|
|
|
|
|
|
322
|
|
|
|
|
|
|
Note |
|
323
|
|
|
|
|
|
|
|
|
324
|
|
|
|
|
|
|
The following events are flagged: |
|
325
|
|
|
|
|
|
|
Synonymous All codons at a site code for the same amino |
|
326
|
|
|
|
|
|
|
acid. |
|
327
|
|
|
|
|
|
|
Single character Only one sequence at the site is ungapped, |
|
328
|
|
|
|
|
|
|
the result of a recent insertion for example. |
|
329
|
|
|
|
|
|
|
All gaps All sequences at a site contain a gap |
|
330
|
|
|
|
|
|
|
character. |
|
331
|
|
|
|
|
|
|
|
|
332
|
|
|
|
|
|
|
Sites marked "Single character" or "All gaps" are not counted |
|
333
|
|
|
|
|
|
|
towards the number of sites for the purposes of correcting for |
|
334
|
|
|
|
|
|
|
multiple comparisons since it is not possible to detect selection |
|
335
|
|
|
|
|
|
|
from none or one observation under the assumptions made by the |
|
336
|
|
|
|
|
|
|
sitewise likelihood ratio test. |
|
337
|
|
|
|
|
|
|
|
|
338
|
|
|
|
|
|
|
=cut |
|
339
|
|
|
|
|
|
|
|
|
340
|
|
|
|
|
|
|
|
|
341
|
|
|
|
|
|
|
#' keep my emacs happy |
|
342
|
|
|
|
|
|
|
|
|
343
|
|
|
|
|
|
|
BEGIN { |
|
344
|
|
|
|
|
|
|
|
|
345
|
1
|
|
|
1
|
|
2
|
$MINNAMELEN = 25; |
|
346
|
1
|
|
|
|
|
1
|
$PROGRAMNAME = 'Slr_Linux_static'; |
|
347
|
1
|
50
|
|
|
|
8
|
if ($^O =~ /darwin/i) { |
|
|
|
50
|
|
|
|
|
|
|
348
|
0
|
|
|
|
|
0
|
$PROGRAMNAME = 'Slr_osx'; |
|
349
|
|
|
|
|
|
|
} elsif ($^O =~ /mswin/i) { |
|
350
|
0
|
|
|
|
|
0
|
$PROGRAMNAME = 'Slr_windows.exe'; |
|
351
|
|
|
|
|
|
|
} |
|
352
|
1
|
50
|
|
|
|
3
|
if( defined $ENV{'SLRDIR'} ) { |
|
353
|
0
|
0
|
|
|
|
0
|
$PROGRAM = Bio::Root::IO->catfile($ENV{'SLRDIR'},$PROGRAMNAME). ($^O =~ /mswin/i ?'_windows.exe':'');; |
|
354
|
|
|
|
|
|
|
} |
|
355
|
|
|
|
|
|
|
|
|
356
|
|
|
|
|
|
|
# valid values for parameters, the default one is always |
|
357
|
|
|
|
|
|
|
# the first one in the array |
|
358
|
|
|
|
|
|
|
# example file provided with the package |
|
359
|
|
|
|
|
|
|
%VALIDVALUES = ( |
|
360
|
1
|
|
|
|
|
1594
|
'outfile' => 'slr.res', |
|
361
|
|
|
|
|
|
|
'reoptimise' => [ 1,0], |
|
362
|
|
|
|
|
|
|
'kappa' => '2.0', |
|
363
|
|
|
|
|
|
|
'omega' => '0.1', |
|
364
|
|
|
|
|
|
|
'codonf' => [ 0, 1,2], |
|
365
|
|
|
|
|
|
|
'freqtype' => [ 0, 1,2,3], |
|
366
|
|
|
|
|
|
|
'positive_only' => [ 0, 1], |
|
367
|
|
|
|
|
|
|
'gencode' => [ "universal", "mammalian"], |
|
368
|
|
|
|
|
|
|
'nucleof' => [ 0, 1], |
|
369
|
|
|
|
|
|
|
'aminof' => [ 0, 1,2], |
|
370
|
|
|
|
|
|
|
'nucfile' => '', |
|
371
|
|
|
|
|
|
|
'aminofile' => '', |
|
372
|
|
|
|
|
|
|
'timemem' => [ 0, 1], |
|
373
|
|
|
|
|
|
|
'ldiff' => [ 3.841459, 0.4549364,1.323304,2.705543,6.634897,7.879439,10.82757], |
|
374
|
|
|
|
|
|
|
'paramin' => '', |
|
375
|
|
|
|
|
|
|
'paramout' => '', |
|
376
|
|
|
|
|
|
|
'skipsitewise' => [ 0, 1], |
|
377
|
|
|
|
|
|
|
'seed' => [0], |
|
378
|
|
|
|
|
|
|
'saveseed' => [ 1, 0] |
|
379
|
|
|
|
|
|
|
); |
|
380
|
|
|
|
|
|
|
} |
|
381
|
|
|
|
|
|
|
|
|
382
|
|
|
|
|
|
|
=head2 program_name |
|
383
|
|
|
|
|
|
|
|
|
384
|
|
|
|
|
|
|
Title : program_name |
|
385
|
|
|
|
|
|
|
Usage : $factory->program_name() |
|
386
|
|
|
|
|
|
|
Function: holds the program name |
|
387
|
|
|
|
|
|
|
Returns: string |
|
388
|
|
|
|
|
|
|
Args : None |
|
389
|
|
|
|
|
|
|
|
|
390
|
|
|
|
|
|
|
=cut |
|
391
|
|
|
|
|
|
|
|
|
392
|
|
|
|
|
|
|
sub program_name { |
|
393
|
6
|
|
|
6
|
1
|
26
|
return $PROGRAMNAME; |
|
394
|
|
|
|
|
|
|
} |
|
395
|
|
|
|
|
|
|
|
|
396
|
|
|
|
|
|
|
=head2 program_dir |
|
397
|
|
|
|
|
|
|
|
|
398
|
|
|
|
|
|
|
Title : program_dir |
|
399
|
|
|
|
|
|
|
Usage : ->program_dir() |
|
400
|
|
|
|
|
|
|
Function: returns the program directory, obtained from ENV variable. |
|
401
|
|
|
|
|
|
|
Returns: string |
|
402
|
|
|
|
|
|
|
Args : |
|
403
|
|
|
|
|
|
|
|
|
404
|
|
|
|
|
|
|
=cut |
|
405
|
|
|
|
|
|
|
|
|
406
|
|
|
|
|
|
|
sub program_dir { |
|
407
|
3
|
50
|
|
3
|
1
|
14
|
return Bio::Root::IO->catfile($ENV{SLRDIR}) if $ENV{SLRDIR}; |
|
408
|
|
|
|
|
|
|
} |
|
409
|
|
|
|
|
|
|
|
|
410
|
|
|
|
|
|
|
|
|
411
|
|
|
|
|
|
|
=head2 new |
|
412
|
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
Title : new |
|
414
|
|
|
|
|
|
|
Usage : my $obj = Bio::Tools::Run::Phylo::SLR->new(); |
|
415
|
|
|
|
|
|
|
Function: Builds a new Bio::Tools::Run::Phylo::SLR object |
|
416
|
|
|
|
|
|
|
Returns : Bio::Tools::Run::Phylo::SLR |
|
417
|
|
|
|
|
|
|
Args : -alignment => the Bio::Align::AlignI object |
|
418
|
|
|
|
|
|
|
-save_tempfiles => boolean to save the generated tempfiles and |
|
419
|
|
|
|
|
|
|
NOT cleanup after onesself (default FALSE) |
|
420
|
|
|
|
|
|
|
-tree => the Bio::Tree::TreeI object |
|
421
|
|
|
|
|
|
|
-params => a hashref of SLR parameters (all passed to set_parameter) |
|
422
|
|
|
|
|
|
|
-executable => where the SLR executable resides |
|
423
|
|
|
|
|
|
|
|
|
424
|
|
|
|
|
|
|
See also: L, L |
|
425
|
|
|
|
|
|
|
|
|
426
|
|
|
|
|
|
|
=cut |
|
427
|
|
|
|
|
|
|
|
|
428
|
|
|
|
|
|
|
sub new { |
|
429
|
1
|
|
|
1
|
1
|
98
|
my($class,@args) = @_; |
|
430
|
|
|
|
|
|
|
|
|
431
|
1
|
|
|
|
|
10
|
my $self = $class->SUPER::new(@args); |
|
432
|
1
|
|
|
|
|
17
|
my ($aln, $tree, $st, $params, $exe, |
|
433
|
|
|
|
|
|
|
$ubl) = $self->_rearrange([qw(ALIGNMENT TREE SAVE_TEMPFILES |
|
434
|
|
|
|
|
|
|
PARAMS EXECUTABLE)], |
|
435
|
|
|
|
|
|
|
@args); |
|
436
|
1
|
50
|
|
|
|
9
|
defined $aln && $self->alignment($aln); |
|
437
|
1
|
50
|
|
|
|
2
|
defined $tree && $self->tree($tree); |
|
438
|
1
|
50
|
|
|
|
2
|
defined $st && $self->save_tempfiles($st); |
|
439
|
1
|
50
|
|
|
|
2
|
defined $exe && $self->executable($exe); |
|
440
|
|
|
|
|
|
|
|
|
441
|
1
|
|
|
|
|
3
|
$self->set_default_parameters(); |
|
442
|
1
|
50
|
|
|
|
3
|
if( defined $params ) { |
|
443
|
0
|
0
|
|
|
|
0
|
if( ref($params) !~ /HASH/i ) { |
|
444
|
0
|
|
|
|
|
0
|
$self->warn("Must provide a valid hash ref for parameter -FLAGS"); |
|
445
|
|
|
|
|
|
|
} else { |
|
446
|
0
|
|
|
|
|
0
|
map { $self->set_parameter($_, $$params{$_}) } keys %$params; |
|
|
0
|
|
|
|
|
0
|
|
|
447
|
|
|
|
|
|
|
} |
|
448
|
|
|
|
|
|
|
} |
|
449
|
1
|
|
|
|
|
4
|
return $self; |
|
450
|
|
|
|
|
|
|
} |
|
451
|
|
|
|
|
|
|
|
|
452
|
|
|
|
|
|
|
|
|
453
|
|
|
|
|
|
|
=head2 prepare |
|
454
|
|
|
|
|
|
|
|
|
455
|
|
|
|
|
|
|
Title : prepare |
|
456
|
|
|
|
|
|
|
Usage : my $rundir = $slr->prepare($aln); |
|
457
|
|
|
|
|
|
|
Function: prepare the SLR analysis using the default or updated parameters |
|
458
|
|
|
|
|
|
|
the alignment parameter must have been set |
|
459
|
|
|
|
|
|
|
Returns : value of rundir |
|
460
|
|
|
|
|
|
|
Args : L object, |
|
461
|
|
|
|
|
|
|
L object |
|
462
|
|
|
|
|
|
|
|
|
463
|
|
|
|
|
|
|
=cut |
|
464
|
|
|
|
|
|
|
|
|
465
|
|
|
|
|
|
|
sub prepare{ |
|
466
|
0
|
|
|
0
|
1
|
0
|
my ($self,$aln,$tree) = @_; |
|
467
|
0
|
0
|
|
|
|
0
|
unless ( $self->save_tempfiles ) { |
|
468
|
|
|
|
|
|
|
# brush so we don't get plaque buildup ;) |
|
469
|
0
|
|
|
|
|
0
|
$self->cleanup(); |
|
470
|
|
|
|
|
|
|
} |
|
471
|
0
|
0
|
|
|
|
0
|
$tree = $self->tree unless $tree; |
|
472
|
0
|
0
|
|
|
|
0
|
$aln = $self->alignment unless $aln; |
|
473
|
0
|
0
|
|
|
|
0
|
if( ! $aln ) { |
|
474
|
0
|
|
|
|
|
0
|
$self->warn("must have supplied a valid alignment file in order to run SLR"); |
|
475
|
0
|
|
|
|
|
0
|
return 0; |
|
476
|
|
|
|
|
|
|
} |
|
477
|
0
|
0
|
|
|
|
0
|
if( ! $tree ) { |
|
478
|
0
|
|
|
|
|
0
|
$self->warn("must have supplied a valid tree file in order to run SLR"); |
|
479
|
0
|
|
|
|
|
0
|
return 0; |
|
480
|
|
|
|
|
|
|
} |
|
481
|
0
|
|
|
|
|
0
|
my ($tempdir) = $self->tempdir(); |
|
482
|
0
|
|
|
|
|
0
|
my ($tempseqFH,$tempseqfile); |
|
483
|
|
|
|
|
|
|
|
|
484
|
|
|
|
|
|
|
# Reorder the alignment according to the tree |
|
485
|
0
|
|
|
|
|
0
|
my $ct = 1; |
|
486
|
0
|
|
|
|
|
0
|
my %order; |
|
487
|
0
|
|
|
|
|
0
|
foreach my $node ($tree->get_leaf_nodes) { |
|
488
|
0
|
|
|
|
|
0
|
$order{$node->id_output} = $ct++; |
|
489
|
|
|
|
|
|
|
} |
|
490
|
0
|
|
|
|
|
0
|
my @seq; my @ids; |
|
491
|
0
|
|
|
|
|
0
|
foreach my $seq ( $aln->each_seq() ) { |
|
492
|
0
|
|
|
|
|
0
|
push @seq, $seq; |
|
493
|
0
|
|
|
|
|
0
|
push @ids, $seq->display_id; |
|
494
|
|
|
|
|
|
|
} |
|
495
|
|
|
|
|
|
|
# use the map-sort-map idiom: |
|
496
|
0
|
|
|
|
|
0
|
my @sorted = map { $_->[1] } sort { $a->[0] <=> $b->[0] } map { [$order{$_->id()}, $_] } @seq; |
|
|
0
|
|
|
|
|
0
|
|
|
|
0
|
|
|
|
|
0
|
|
|
|
0
|
|
|
|
|
0
|
|
|
497
|
0
|
|
|
|
|
0
|
my $sorted_aln = Bio::SimpleAlign->new(); |
|
498
|
0
|
|
|
|
|
0
|
foreach (@sorted) { |
|
499
|
0
|
|
|
|
|
0
|
$sorted_aln->add_seq($_); |
|
500
|
|
|
|
|
|
|
} |
|
501
|
|
|
|
|
|
|
|
|
502
|
|
|
|
|
|
|
# Rename the leaf nodes in the tree from 1 to n |
|
503
|
0
|
|
|
|
|
0
|
$ct = 1; |
|
504
|
0
|
|
|
|
|
0
|
foreach my $node ($tree->get_leaf_nodes) { |
|
505
|
0
|
|
|
|
|
0
|
$node->id($ct++); |
|
506
|
|
|
|
|
|
|
} |
|
507
|
|
|
|
|
|
|
|
|
508
|
0
|
0
|
|
|
|
0
|
($tempseqFH,$tempseqfile) = $self->io->tempfile |
|
509
|
|
|
|
|
|
|
('-dir' => $tempdir, |
|
510
|
|
|
|
|
|
|
UNLINK => ($self->save_tempfiles ? 0 : 1)); |
|
511
|
0
|
0
|
|
|
|
0
|
my $alnout = Bio::AlignIO->new('-format' => 'phylip', |
|
512
|
|
|
|
|
|
|
'-fh' => $tempseqFH, |
|
513
|
|
|
|
|
|
|
'-interleaved' => 0, |
|
514
|
|
|
|
|
|
|
'-idlinebreak' => 1, |
|
515
|
|
|
|
|
|
|
'-idlength' => $MINNAMELEN > $aln->maxdisplayname_length() ? $MINNAMELEN : $aln->maxdisplayname_length() +1); |
|
516
|
|
|
|
|
|
|
|
|
517
|
0
|
|
|
|
|
0
|
$alnout->write_aln($sorted_aln); |
|
518
|
0
|
|
|
|
|
0
|
$alnout->close(); |
|
519
|
0
|
|
|
|
|
0
|
undef $alnout; |
|
520
|
0
|
|
|
|
|
0
|
close($tempseqFH); |
|
521
|
|
|
|
|
|
|
|
|
522
|
0
|
|
|
|
|
0
|
my ($temptreeFH,$temptreefile); |
|
523
|
0
|
0
|
|
|
|
0
|
($temptreeFH,$temptreefile) = $self->io->tempfile |
|
524
|
|
|
|
|
|
|
('-dir' => $tempdir, |
|
525
|
|
|
|
|
|
|
UNLINK => ($self->save_tempfiles ? 0 : 1)); |
|
526
|
|
|
|
|
|
|
|
|
527
|
0
|
|
|
|
|
0
|
my $treeout = Bio::TreeIO->new('-format' => 'newick', |
|
528
|
|
|
|
|
|
|
'-fh' => $temptreeFH); |
|
529
|
|
|
|
|
|
|
|
|
530
|
|
|
|
|
|
|
# We need to add a line with the num of leaves ($ct-1) and the |
|
531
|
|
|
|
|
|
|
# num of trees (1) |
|
532
|
0
|
|
|
|
|
0
|
$treeout->_print(sprintf("%d 1\n",($ct-1))); |
|
533
|
0
|
|
|
|
|
0
|
$treeout->write_tree($tree); |
|
534
|
0
|
|
|
|
|
0
|
$treeout->close(); |
|
535
|
0
|
|
|
|
|
0
|
close($temptreeFH); |
|
536
|
|
|
|
|
|
|
|
|
537
|
|
|
|
|
|
|
# now let's print the ctl file. |
|
538
|
|
|
|
|
|
|
# many of the these programs are finicky about what the filename is |
|
539
|
|
|
|
|
|
|
# and won't even run without the properly named file. |
|
540
|
|
|
|
|
|
|
|
|
541
|
0
|
|
|
|
|
0
|
my ($treevolume,$treedirectories,$treefile) = File::Spec->splitpath( $temptreefile ); |
|
542
|
0
|
|
|
|
|
0
|
my ($alnvolume,$alndirectories,$alnfile) = File::Spec->splitpath( $tempseqfile ); |
|
543
|
0
|
|
|
|
|
0
|
my $slr_ctl = "$tempdir/slr.ctl"; |
|
544
|
0
|
0
|
|
|
|
0
|
open(SLR, ">$slr_ctl") or $self->throw("cannot open $slr_ctl for writing"); |
|
545
|
0
|
|
|
|
|
0
|
print SLR "seqfile\: $alnfile\n"; |
|
546
|
0
|
|
|
|
|
0
|
print SLR "treefile\: $treefile\n"; |
|
547
|
0
|
|
|
|
|
0
|
my $outfile = $self->outfile_name; |
|
548
|
0
|
|
|
|
|
0
|
print SLR "outfile\: $outfile\n"; |
|
549
|
|
|
|
|
|
|
|
|
550
|
0
|
|
|
|
|
0
|
my %params = $self->get_parameters; |
|
551
|
0
|
|
|
|
|
0
|
while( my ($param,$val) = each %params ) { |
|
552
|
0
|
0
|
|
|
|
0
|
next if $param eq 'outfile'; |
|
553
|
0
|
|
|
|
|
0
|
print SLR "$param\: $val\n"; |
|
554
|
|
|
|
|
|
|
} |
|
555
|
0
|
|
|
|
|
0
|
close(SLR); |
|
556
|
0
|
|
|
|
|
0
|
return $tempdir; |
|
557
|
|
|
|
|
|
|
} |
|
558
|
|
|
|
|
|
|
|
|
559
|
|
|
|
|
|
|
|
|
560
|
|
|
|
|
|
|
|
|
561
|
|
|
|
|
|
|
=head2 run |
|
562
|
|
|
|
|
|
|
|
|
563
|
|
|
|
|
|
|
Title : run |
|
564
|
|
|
|
|
|
|
Usage : my ($rc,$parser) = $slr->run($aln,$tree); |
|
565
|
|
|
|
|
|
|
Function: run the SLR analysis using the default or updated parameters |
|
566
|
|
|
|
|
|
|
the alignment parameter must have been set |
|
567
|
|
|
|
|
|
|
Returns : Return code, L |
|
568
|
|
|
|
|
|
|
Args : L object, |
|
569
|
|
|
|
|
|
|
L object |
|
570
|
|
|
|
|
|
|
|
|
571
|
|
|
|
|
|
|
|
|
572
|
|
|
|
|
|
|
=cut |
|
573
|
|
|
|
|
|
|
|
|
574
|
|
|
|
|
|
|
sub run { |
|
575
|
0
|
|
|
0
|
1
|
0
|
my ($self) = shift;; |
|
576
|
0
|
|
|
|
|
0
|
my $outfile = $self->outfile_name; |
|
577
|
0
|
|
|
|
|
0
|
my $tmpdir = $self->prepare(@_); |
|
578
|
|
|
|
|
|
|
|
|
579
|
|
|
|
|
|
|
#my ($rc,$parser) = (1); |
|
580
|
0
|
|
|
|
|
0
|
my ($rc,$results) = (1); |
|
581
|
|
|
|
|
|
|
{ |
|
582
|
0
|
|
|
|
|
0
|
my $cwd = cwd(); |
|
|
0
|
|
|
|
|
0
|
|
|
583
|
0
|
|
|
|
|
0
|
my $exit_status; |
|
584
|
0
|
|
|
|
|
0
|
chdir($tmpdir); |
|
585
|
0
|
|
|
|
|
0
|
my $slrexe = $self->executable(); |
|
586
|
0
|
0
|
0
|
|
|
0
|
$self->throw("unable to find or run executable for SLR") unless $slrexe && -e $slrexe && -x _; |
|
|
|
|
0
|
|
|
|
|
|
587
|
0
|
|
|
|
|
0
|
my $run; |
|
588
|
0
|
0
|
|
|
|
0
|
open($run, "$slrexe |") or $self->throw("Cannot open exe $slrexe"); |
|
589
|
0
|
|
|
|
|
0
|
my @output = <$run>; |
|
590
|
0
|
|
|
|
|
0
|
$exit_status = close($run); |
|
591
|
0
|
|
|
|
|
0
|
$self->error_string(join('',@output)); |
|
592
|
0
|
0
|
0
|
|
|
0
|
if( (grep { /\berr(or)?: /io } @output) || !$exit_status) { |
|
|
0
|
|
|
|
|
0
|
|
|
593
|
0
|
|
|
|
|
0
|
$self->warn("There was an error - see error_string for the program output"); |
|
594
|
0
|
|
|
|
|
0
|
$rc = 0; |
|
595
|
|
|
|
|
|
|
} |
|
596
|
0
|
|
|
|
|
0
|
eval { |
|
597
|
0
|
0
|
|
|
|
0
|
open RESULTS, "$tmpdir/$outfile" or die "couldnt open results file: $!\n"; |
|
598
|
0
|
|
|
|
|
0
|
my $okay = 0; |
|
599
|
0
|
|
|
|
|
0
|
my $sites; |
|
600
|
0
|
|
|
|
|
0
|
my $type = 'default'; |
|
601
|
0
|
|
|
|
|
0
|
while () { |
|
602
|
0
|
|
|
|
|
0
|
chomp $_; |
|
603
|
0
|
0
|
|
|
|
0
|
if ( /^\#/ ) {next;} |
|
|
0
|
|
|
|
|
0
|
|
|
604
|
0
|
0
|
|
|
|
0
|
if ( /\!/ ) {$type = 'random';} # random is last |
|
|
0
|
0
|
|
|
|
0
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
605
|
0
|
|
|
|
|
0
|
elsif ( /\+/ ) {$type = 'positive';} |
|
606
|
0
|
|
|
|
|
0
|
elsif ( /\-\s+/ ) {$type = 'negative';} |
|
607
|
0
|
|
|
|
|
0
|
elsif ( /Constant/ ) {$type = 'constant';} |
|
608
|
0
|
|
|
|
|
0
|
elsif ( /All gaps/ ) {$type = 'all_gaps';} |
|
609
|
0
|
|
|
|
|
0
|
elsif ( /Single character/ ) {$type = 'single_character';} |
|
610
|
0
|
|
|
|
|
0
|
elsif ( /Synonymous/ ) {$type = 'synonymous';} |
|
611
|
0
|
|
|
|
|
0
|
else {$type = 'default'} |
|
612
|
0
|
0
|
|
|
|
0
|
if ( /^\s+(\d+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/ ) { |
|
613
|
0
|
|
|
|
|
0
|
push @{$sites->{$type}}, [$1,$2,$3,$4,$5,$6,$7,$8,$9]; |
|
|
0
|
|
|
|
|
0
|
|
|
614
|
|
|
|
|
|
|
} else { |
|
615
|
0
|
|
|
|
|
0
|
$DB::single=1;1; |
|
|
0
|
|
|
|
|
0
|
|
|
616
|
|
|
|
|
|
|
} |
|
617
|
|
|
|
|
|
|
} |
|
618
|
0
|
|
|
|
|
0
|
$results = $sites; |
|
619
|
0
|
|
|
|
|
0
|
close RESULTS; |
|
620
|
|
|
|
|
|
|
# TODO: we could have a proper parser object |
|
621
|
|
|
|
|
|
|
# $parser = Bio::Tools::Phylo::SLR->new(-file => "$tmpdir/$outfile", |
|
622
|
|
|
|
|
|
|
# -dir => "$tmpdir"); |
|
623
|
|
|
|
|
|
|
}; |
|
624
|
0
|
0
|
|
|
|
0
|
if( $@ ) { |
|
625
|
0
|
|
|
|
|
0
|
$self->warn($self->error_string); |
|
626
|
|
|
|
|
|
|
} |
|
627
|
0
|
|
|
|
|
0
|
chdir($cwd); |
|
628
|
|
|
|
|
|
|
} |
|
629
|
|
|
|
|
|
|
# return ($rc,$parser); |
|
630
|
0
|
|
|
|
|
0
|
return ($rc,$results); |
|
631
|
|
|
|
|
|
|
} |
|
632
|
|
|
|
|
|
|
|
|
633
|
|
|
|
|
|
|
=head2 error_string |
|
634
|
|
|
|
|
|
|
|
|
635
|
|
|
|
|
|
|
Title : error_string |
|
636
|
|
|
|
|
|
|
Usage : $obj->error_string($newval) |
|
637
|
|
|
|
|
|
|
Function: Where the output from the last analysus run is stored. |
|
638
|
|
|
|
|
|
|
Returns : value of error_string |
|
639
|
|
|
|
|
|
|
Args : newvalue (optional) |
|
640
|
|
|
|
|
|
|
|
|
641
|
|
|
|
|
|
|
|
|
642
|
|
|
|
|
|
|
=cut |
|
643
|
|
|
|
|
|
|
|
|
644
|
|
|
|
|
|
|
sub error_string{ |
|
645
|
0
|
|
|
0
|
1
|
0
|
my ($self,$value) = @_; |
|
646
|
0
|
0
|
|
|
|
0
|
if( defined $value) { |
|
647
|
0
|
|
|
|
|
0
|
$self->{'error_string'} = $value; |
|
648
|
|
|
|
|
|
|
} |
|
649
|
0
|
|
|
|
|
0
|
return $self->{'error_string'}; |
|
650
|
|
|
|
|
|
|
|
|
651
|
|
|
|
|
|
|
} |
|
652
|
|
|
|
|
|
|
|
|
653
|
|
|
|
|
|
|
=head2 alignment |
|
654
|
|
|
|
|
|
|
|
|
655
|
|
|
|
|
|
|
Title : alignment |
|
656
|
|
|
|
|
|
|
Usage : $slr->align($aln); |
|
657
|
|
|
|
|
|
|
Function: Get/Set the L object |
|
658
|
|
|
|
|
|
|
Returns : L object |
|
659
|
|
|
|
|
|
|
Args : [optional] L |
|
660
|
|
|
|
|
|
|
Comment : We could potentially add support for running directly on a file |
|
661
|
|
|
|
|
|
|
but we shall keep it simple |
|
662
|
|
|
|
|
|
|
See also: L |
|
663
|
|
|
|
|
|
|
|
|
664
|
|
|
|
|
|
|
=cut |
|
665
|
|
|
|
|
|
|
|
|
666
|
|
|
|
|
|
|
sub alignment{ |
|
667
|
0
|
|
|
0
|
1
|
0
|
my ($self,$aln) = @_; |
|
668
|
|
|
|
|
|
|
|
|
669
|
0
|
0
|
|
|
|
0
|
if( defined $aln ) { |
|
670
|
0
|
0
|
0
|
|
|
0
|
if( -e $aln ) { |
|
|
|
0
|
|
|
|
|
|
|
671
|
0
|
|
|
|
|
0
|
$self->{'_alignment'} = $aln; |
|
672
|
|
|
|
|
|
|
} elsif( !ref($aln) || ! $aln->isa('Bio::Align::AlignI') ) { |
|
673
|
0
|
|
|
|
|
0
|
$self->warn("Must specify a valid Bio::Align::AlignI object to the alignment function not $aln"); |
|
674
|
0
|
|
|
|
|
0
|
return undef; |
|
675
|
|
|
|
|
|
|
} else { |
|
676
|
0
|
|
|
|
|
0
|
$self->{'_alignment'} = $aln; |
|
677
|
|
|
|
|
|
|
} |
|
678
|
|
|
|
|
|
|
} |
|
679
|
0
|
|
|
|
|
0
|
return $self->{'_alignment'}; |
|
680
|
|
|
|
|
|
|
} |
|
681
|
|
|
|
|
|
|
|
|
682
|
|
|
|
|
|
|
=head2 tree |
|
683
|
|
|
|
|
|
|
|
|
684
|
|
|
|
|
|
|
Title : tree |
|
685
|
|
|
|
|
|
|
Usage : $slr->tree($tree, %params); |
|
686
|
|
|
|
|
|
|
Function: Get/Set the L object |
|
687
|
|
|
|
|
|
|
Returns : L |
|
688
|
|
|
|
|
|
|
Args : [optional] $tree => L, |
|
689
|
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
Comment : We could potentially add support for running directly on a file |
|
691
|
|
|
|
|
|
|
but we shall keep it simple |
|
692
|
|
|
|
|
|
|
See also: L |
|
693
|
|
|
|
|
|
|
|
|
694
|
|
|
|
|
|
|
=cut |
|
695
|
|
|
|
|
|
|
|
|
696
|
|
|
|
|
|
|
sub tree { |
|
697
|
0
|
|
|
0
|
1
|
0
|
my ($self, $tree, %params) = @_; |
|
698
|
0
|
0
|
|
|
|
0
|
if( defined $tree ) { |
|
699
|
0
|
0
|
0
|
|
|
0
|
if( ! ref($tree) || ! $tree->isa('Bio::Tree::TreeI') ) { |
|
700
|
0
|
|
|
|
|
0
|
$self->warn("Must specify a valid Bio::Tree::TreeI object to the alignment function"); |
|
701
|
|
|
|
|
|
|
} |
|
702
|
0
|
|
|
|
|
0
|
$self->{'_tree'} = $tree; |
|
703
|
|
|
|
|
|
|
} |
|
704
|
0
|
|
|
|
|
0
|
return $self->{'_tree'}; |
|
705
|
|
|
|
|
|
|
} |
|
706
|
|
|
|
|
|
|
|
|
707
|
|
|
|
|
|
|
=head2 get_parameters |
|
708
|
|
|
|
|
|
|
|
|
709
|
|
|
|
|
|
|
Title : get_parameters |
|
710
|
|
|
|
|
|
|
Usage : my %params = $self->get_parameters(); |
|
711
|
|
|
|
|
|
|
Function: returns the list of parameters as a hash |
|
712
|
|
|
|
|
|
|
Returns : associative array keyed on parameter names |
|
713
|
|
|
|
|
|
|
Args : none |
|
714
|
|
|
|
|
|
|
|
|
715
|
|
|
|
|
|
|
|
|
716
|
|
|
|
|
|
|
=cut |
|
717
|
|
|
|
|
|
|
|
|
718
|
|
|
|
|
|
|
sub get_parameters{ |
|
719
|
0
|
|
|
0
|
1
|
0
|
my ($self) = @_; |
|
720
|
|
|
|
|
|
|
# we're returning a copy of this |
|
721
|
0
|
|
|
|
|
0
|
return %{ $self->{'_slrparams'} }; |
|
|
0
|
|
|
|
|
0
|
|
|
722
|
|
|
|
|
|
|
} |
|
723
|
|
|
|
|
|
|
|
|
724
|
|
|
|
|
|
|
|
|
725
|
|
|
|
|
|
|
=head2 set_parameter |
|
726
|
|
|
|
|
|
|
|
|
727
|
|
|
|
|
|
|
Title : set_parameter |
|
728
|
|
|
|
|
|
|
Usage : $slr->set_parameter($param,$val); |
|
729
|
|
|
|
|
|
|
Function: Sets a SLR parameter, will be validated against |
|
730
|
|
|
|
|
|
|
the valid values as set in the %VALIDVALUES class variable. |
|
731
|
|
|
|
|
|
|
The checks can be ignored if one turns off param checks like this: |
|
732
|
|
|
|
|
|
|
$slr->no_param_checks(1) |
|
733
|
|
|
|
|
|
|
Returns : boolean if set was success, if verbose is set to -1 |
|
734
|
|
|
|
|
|
|
then no warning will be reported |
|
735
|
|
|
|
|
|
|
Args : $param => name of the parameter |
|
736
|
|
|
|
|
|
|
$value => value to set the parameter to |
|
737
|
|
|
|
|
|
|
See also: L |
|
738
|
|
|
|
|
|
|
|
|
739
|
|
|
|
|
|
|
=cut |
|
740
|
|
|
|
|
|
|
|
|
741
|
|
|
|
|
|
|
sub set_parameter{ |
|
742
|
0
|
|
|
0
|
1
|
0
|
my ($self,$param,$value) = @_; |
|
743
|
0
|
0
|
0
|
|
|
0
|
unless (defined $self->{'no_param_checks'} && $self->{'no_param_checks'} == 1) { |
|
744
|
0
|
0
|
|
|
|
0
|
if ( ! defined $VALIDVALUES{$param} ) { |
|
745
|
0
|
|
|
|
|
0
|
$self->warn("unknown parameter $param will not be set unless you force by setting no_param_checks to true"); |
|
746
|
0
|
|
|
|
|
0
|
return 0; |
|
747
|
|
|
|
|
|
|
} |
|
748
|
0
|
0
|
0
|
|
|
0
|
if ( ref( $VALIDVALUES{$param}) =~ /ARRAY/i && |
|
749
|
0
|
|
|
|
|
0
|
scalar @{$VALIDVALUES{$param}} > 0 ) { |
|
750
|
|
|
|
|
|
|
|
|
751
|
0
|
0
|
|
|
|
0
|
unless ( grep { $value eq $_ } @{ $VALIDVALUES{$param} } ) { |
|
|
0
|
|
|
|
|
0
|
|
|
|
0
|
|
|
|
|
0
|
|
|
752
|
0
|
|
|
|
|
0
|
$self->warn("parameter $param specified value $value is not recognized, please see the documentation and the code for this module or set the no_param_checks to a true value"); |
|
753
|
0
|
|
|
|
|
0
|
return 0; |
|
754
|
|
|
|
|
|
|
} |
|
755
|
|
|
|
|
|
|
} |
|
756
|
|
|
|
|
|
|
} |
|
757
|
0
|
|
|
|
|
0
|
$self->{'_slrparams'}->{$param} = $value; |
|
758
|
0
|
|
|
|
|
0
|
return 1; |
|
759
|
|
|
|
|
|
|
} |
|
760
|
|
|
|
|
|
|
|
|
761
|
|
|
|
|
|
|
=head2 set_default_parameters |
|
762
|
|
|
|
|
|
|
|
|
763
|
|
|
|
|
|
|
Title : set_default_parameters |
|
764
|
|
|
|
|
|
|
Usage : $slr->set_default_parameters(0); |
|
765
|
|
|
|
|
|
|
Function: (Re)set the default parameters from the defaults |
|
766
|
|
|
|
|
|
|
(the first value in each array in the |
|
767
|
|
|
|
|
|
|
%VALIDVALUES class variable) |
|
768
|
|
|
|
|
|
|
Returns : none |
|
769
|
|
|
|
|
|
|
Args : boolean: keep existing parameter values |
|
770
|
|
|
|
|
|
|
|
|
771
|
|
|
|
|
|
|
|
|
772
|
|
|
|
|
|
|
=cut |
|
773
|
|
|
|
|
|
|
|
|
774
|
|
|
|
|
|
|
sub set_default_parameters{ |
|
775
|
1
|
|
|
1
|
1
|
2
|
my ($self,$keepold) = @_; |
|
776
|
1
|
50
|
|
|
|
3
|
$keepold = 0 unless defined $keepold; |
|
777
|
|
|
|
|
|
|
|
|
778
|
1
|
|
|
|
|
6
|
while( my ($param,$val) = each %VALIDVALUES ) { |
|
779
|
|
|
|
|
|
|
# skip if we want to keep old values and it is already set |
|
780
|
19
|
50
|
33
|
|
|
48
|
next if( defined $self->{'_slrparams'}->{$param} && $keepold); |
|
781
|
19
|
100
|
|
|
|
40
|
if(ref($val)=~/ARRAY/i ) { |
|
782
|
12
|
|
|
|
|
31
|
$self->{'_slrparams'}->{$param} = $val->[0]; |
|
783
|
|
|
|
|
|
|
} else { |
|
784
|
7
|
|
|
|
|
19
|
$self->{'_slrparams'}->{$param} = $val; |
|
785
|
|
|
|
|
|
|
} |
|
786
|
|
|
|
|
|
|
} |
|
787
|
|
|
|
|
|
|
} |
|
788
|
|
|
|
|
|
|
|
|
789
|
|
|
|
|
|
|
|
|
790
|
|
|
|
|
|
|
=head1 Bio::Tools::Run::WrapperBase methods |
|
791
|
|
|
|
|
|
|
|
|
792
|
|
|
|
|
|
|
=cut |
|
793
|
|
|
|
|
|
|
|
|
794
|
|
|
|
|
|
|
=head2 no_param_checks |
|
795
|
|
|
|
|
|
|
|
|
796
|
|
|
|
|
|
|
Title : no_param_checks |
|
797
|
|
|
|
|
|
|
Usage : $obj->no_param_checks($newval) |
|
798
|
|
|
|
|
|
|
Function: Boolean flag as to whether or not we should |
|
799
|
|
|
|
|
|
|
trust the sanity checks for parameter values |
|
800
|
|
|
|
|
|
|
Returns : value of no_param_checks |
|
801
|
|
|
|
|
|
|
Args : newvalue (optional) |
|
802
|
|
|
|
|
|
|
|
|
803
|
|
|
|
|
|
|
|
|
804
|
|
|
|
|
|
|
=cut |
|
805
|
|
|
|
|
|
|
|
|
806
|
|
|
|
|
|
|
sub no_param_checks{ |
|
807
|
0
|
|
|
0
|
1
|
0
|
my ($self,$value) = @_; |
|
808
|
0
|
0
|
|
|
|
0
|
if( defined $value) { |
|
809
|
0
|
|
|
|
|
0
|
$self->{'no_param_checks'} = $value; |
|
810
|
|
|
|
|
|
|
} |
|
811
|
0
|
|
|
|
|
0
|
return $self->{'no_param_checks'}; |
|
812
|
|
|
|
|
|
|
} |
|
813
|
|
|
|
|
|
|
|
|
814
|
|
|
|
|
|
|
|
|
815
|
|
|
|
|
|
|
=head2 save_tempfiles |
|
816
|
|
|
|
|
|
|
|
|
817
|
|
|
|
|
|
|
Title : save_tempfiles |
|
818
|
|
|
|
|
|
|
Usage : $obj->save_tempfiles($newval) |
|
819
|
|
|
|
|
|
|
Function: |
|
820
|
|
|
|
|
|
|
Returns : value of save_tempfiles |
|
821
|
|
|
|
|
|
|
Args : newvalue (optional) |
|
822
|
|
|
|
|
|
|
|
|
823
|
|
|
|
|
|
|
|
|
824
|
|
|
|
|
|
|
=cut |
|
825
|
|
|
|
|
|
|
|
|
826
|
|
|
|
|
|
|
=head2 outfile_name |
|
827
|
|
|
|
|
|
|
|
|
828
|
|
|
|
|
|
|
Title : outfile_name |
|
829
|
|
|
|
|
|
|
Usage : my $outfile = $slr->outfile_name(); |
|
830
|
|
|
|
|
|
|
Function: Get/Set the name of the output file for this run |
|
831
|
|
|
|
|
|
|
(if you wanted to do something special) |
|
832
|
|
|
|
|
|
|
Returns : string |
|
833
|
|
|
|
|
|
|
Args : [optional] string to set value to |
|
834
|
|
|
|
|
|
|
|
|
835
|
|
|
|
|
|
|
|
|
836
|
|
|
|
|
|
|
=cut |
|
837
|
|
|
|
|
|
|
|
|
838
|
|
|
|
|
|
|
sub outfile_name { |
|
839
|
0
|
|
|
0
|
1
|
0
|
my $self = shift; |
|
840
|
0
|
0
|
|
|
|
0
|
if( @_ ) { |
|
841
|
0
|
|
|
|
|
0
|
return $self->{'_slrparams'}->{'outfile'} = shift @_; |
|
842
|
|
|
|
|
|
|
} |
|
843
|
0
|
0
|
|
|
|
0
|
unless (defined $self->{'_slrparams'}->{'outfile'}) { |
|
844
|
0
|
|
|
|
|
0
|
$self->{'_slrparams'}->{'outfile'} = 'out.res'; |
|
845
|
|
|
|
|
|
|
} |
|
846
|
0
|
|
|
|
|
0
|
return $self->{'_slrparams'}->{'outfile'}; |
|
847
|
|
|
|
|
|
|
} |
|
848
|
|
|
|
|
|
|
|
|
849
|
|
|
|
|
|
|
=head2 tempdir |
|
850
|
|
|
|
|
|
|
|
|
851
|
|
|
|
|
|
|
Title : tempdir |
|
852
|
|
|
|
|
|
|
Usage : my $tmpdir = $self->tempdir(); |
|
853
|
|
|
|
|
|
|
Function: Retrieve a temporary directory name (which is created) |
|
854
|
|
|
|
|
|
|
Returns : string which is the name of the temporary directory |
|
855
|
|
|
|
|
|
|
Args : none |
|
856
|
|
|
|
|
|
|
|
|
857
|
|
|
|
|
|
|
|
|
858
|
|
|
|
|
|
|
=cut |
|
859
|
|
|
|
|
|
|
|
|
860
|
|
|
|
|
|
|
=head2 cleanup |
|
861
|
|
|
|
|
|
|
|
|
862
|
|
|
|
|
|
|
Title : cleanup |
|
863
|
|
|
|
|
|
|
Usage : $slr->cleanup(); |
|
864
|
|
|
|
|
|
|
Function: Will cleanup the tempdir directory after an SLR run |
|
865
|
|
|
|
|
|
|
Returns : none |
|
866
|
|
|
|
|
|
|
Args : none |
|
867
|
|
|
|
|
|
|
|
|
868
|
|
|
|
|
|
|
|
|
869
|
|
|
|
|
|
|
=cut |
|
870
|
|
|
|
|
|
|
|
|
871
|
|
|
|
|
|
|
=head2 io |
|
872
|
|
|
|
|
|
|
|
|
873
|
|
|
|
|
|
|
Title : io |
|
874
|
|
|
|
|
|
|
Usage : $obj->io($newval) |
|
875
|
|
|
|
|
|
|
Function: Gets a L object |
|
876
|
|
|
|
|
|
|
Returns : L |
|
877
|
|
|
|
|
|
|
Args : none |
|
878
|
|
|
|
|
|
|
|
|
879
|
|
|
|
|
|
|
|
|
880
|
|
|
|
|
|
|
=cut |
|
881
|
|
|
|
|
|
|
|
|
882
|
|
|
|
|
|
|
sub DESTROY { |
|
883
|
1
|
|
|
1
|
|
1768
|
my $self= shift; |
|
884
|
1
|
50
|
|
|
|
8
|
unless ( $self->save_tempfiles ) { |
|
885
|
1
|
|
|
|
|
10
|
$self->cleanup(); |
|
886
|
|
|
|
|
|
|
} |
|
887
|
1
|
|
|
|
|
5
|
$self->SUPER::DESTROY(); |
|
888
|
|
|
|
|
|
|
} |
|
889
|
|
|
|
|
|
|
|
|
890
|
|
|
|
|
|
|
1; |