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; |