line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package Bio::MUST::Core::SeqMask::Rates; |
2
|
|
|
|
|
|
|
# ABSTRACT: Evolutionary rates for sequence sites |
3
|
|
|
|
|
|
|
$Bio::MUST::Core::SeqMask::Rates::VERSION = '0.212530'; |
4
|
17
|
|
|
17
|
|
13398
|
use Moose; |
|
17
|
|
|
|
|
53
|
|
|
17
|
|
|
|
|
135
|
|
5
|
17
|
|
|
17
|
|
123934
|
use namespace::autoclean; |
|
17
|
|
|
|
|
66
|
|
|
17
|
|
|
|
|
153
|
|
6
|
|
|
|
|
|
|
|
7
|
17
|
|
|
17
|
|
1540
|
use autodie; |
|
17
|
|
|
|
|
46
|
|
|
17
|
|
|
|
|
144
|
|
8
|
17
|
|
|
17
|
|
93763
|
use feature qw(say); |
|
17
|
|
|
|
|
52
|
|
|
17
|
|
|
|
|
1458
|
|
9
|
|
|
|
|
|
|
|
10
|
|
|
|
|
|
|
# use Smart::Comments; |
11
|
|
|
|
|
|
|
|
12
|
17
|
|
|
17
|
|
138
|
use Carp; |
|
17
|
|
|
|
|
46
|
|
|
17
|
|
|
|
|
1400
|
|
13
|
17
|
|
|
17
|
|
188
|
use Const::Fast; |
|
17
|
|
|
|
|
41
|
|
|
17
|
|
|
|
|
189
|
|
14
|
17
|
|
|
17
|
|
1369
|
use List::AllUtils qw(each_arrayref); |
|
17
|
|
|
|
|
74
|
|
|
17
|
|
|
|
|
925
|
|
15
|
17
|
|
|
17
|
|
131
|
use POSIX; |
|
17
|
|
|
|
|
53
|
|
|
17
|
|
|
|
|
185
|
|
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
extends 'Bio::MUST::Core::SeqMask'; |
18
|
|
|
|
|
|
|
|
19
|
17
|
|
|
17
|
|
38550
|
use Bio::MUST::Core::Types; |
|
17
|
|
|
|
|
53
|
|
|
17
|
|
|
|
|
660
|
|
20
|
17
|
|
|
17
|
|
128
|
use Bio::MUST::Core::Constants qw(:files); |
|
17
|
|
|
|
|
48
|
|
|
17
|
|
|
|
|
3345
|
|
21
|
17
|
|
|
17
|
|
154
|
use aliased 'Bio::MUST::Core::SeqMask'; |
|
17
|
|
|
|
|
44
|
|
|
17
|
|
|
|
|
165
|
|
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
# override superclass' Bool type |
25
|
|
|
|
|
|
|
# Note: mask indices are as follow: [site] |
26
|
|
|
|
|
|
|
# mask values are rates |
27
|
|
|
|
|
|
|
has '+mask' => ( |
28
|
|
|
|
|
|
|
isa => 'ArrayRef[Num]', |
29
|
|
|
|
|
|
|
); |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
# TODO: mask non-applicable methods from superclass? (Liskov principle) |
32
|
|
|
|
|
|
|
|
33
|
|
|
|
|
|
|
|
34
|
|
|
|
|
|
|
|
35
|
|
|
|
|
|
|
sub min_rate { |
36
|
5
|
|
|
5
|
1
|
17
|
my $self = shift; |
37
|
5
|
|
|
|
|
11
|
return List::AllUtils::min @{ $self->mask }; |
|
5
|
|
|
|
|
179
|
|
38
|
|
|
|
|
|
|
} |
39
|
|
|
|
|
|
|
|
40
|
|
|
|
|
|
|
|
41
|
|
|
|
|
|
|
|
42
|
|
|
|
|
|
|
sub max_rate { |
43
|
5
|
|
|
5
|
1
|
17
|
my $self = shift; |
44
|
5
|
|
|
|
|
10
|
return List::AllUtils::max @{ $self->mask }; |
|
5
|
|
|
|
|
163
|
|
45
|
|
|
|
|
|
|
} |
46
|
|
|
|
|
|
|
|
47
|
|
|
|
|
|
|
|
48
|
|
|
|
|
|
|
|
49
|
|
|
|
|
|
|
sub delta_rates { |
50
|
1
|
|
|
1
|
1
|
12
|
my $self = shift; |
51
|
1
|
|
|
|
|
3
|
my $othr = shift; |
52
|
|
|
|
|
|
|
|
53
|
|
|
|
|
|
|
# check that both rates objects are the same length |
54
|
|
|
|
|
|
|
# potential bugs could come from constant sites etc |
55
|
1
|
|
|
|
|
50
|
my $s_width = $self->mask_len; |
56
|
1
|
|
|
|
|
38
|
my $o_width = $othr->mask_len; |
57
|
1
|
50
|
|
|
|
15
|
carp "[BMC] Warning: Rates widths do not match: $s_width vs. $o_width!" |
58
|
|
|
|
|
|
|
unless $s_width == $o_width; |
59
|
|
|
|
|
|
|
|
60
|
1
|
|
|
|
|
5
|
my @deltas; |
61
|
|
|
|
|
|
|
|
62
|
1
|
|
|
|
|
45
|
my $ea = each_arrayref [ $self->all_states ], [ $othr->all_states ]; |
63
|
1
|
|
|
|
|
3040
|
while (my ($s_rate, $o_rate) = $ea->() ) { |
64
|
22436
|
|
|
|
|
153401
|
push @deltas, 0 + ( sprintf "%.13f", abs( $s_rate - $o_rate ) ); |
65
|
|
|
|
|
|
|
} # Note: trick to get identical results across platforms |
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
# TODO: check that $self->new() is really correct |
68
|
1
|
|
|
|
|
106
|
return $self->new( mask => \@deltas ); |
69
|
|
|
|
|
|
|
} |
70
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
# Rates-based SeqMask factory methods |
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
|
75
|
|
|
|
|
|
|
# small delta for slightly increasing extreme bins |
76
|
|
|
|
|
|
|
const my $DELTA => 1e-13; |
77
|
|
|
|
|
|
|
|
78
|
|
|
|
|
|
|
sub bin_rates_masks { |
79
|
9
|
|
|
9
|
1
|
26
|
my $self = shift; |
80
|
9
|
|
|
|
|
22
|
my $bin_n = shift; |
81
|
9
|
|
50
|
|
|
32
|
my $args = shift // {}; # HashRef (should not be empty...) |
82
|
|
|
|
|
|
|
|
83
|
9
|
|
100
|
|
|
45
|
my $percentile = $args->{percentile} // 0; |
84
|
9
|
|
100
|
|
|
44
|
my $cumulative = $args->{cumulative} // 0; |
85
|
9
|
|
100
|
|
|
39
|
my $descending = $args->{descending} // 0; |
86
|
|
|
|
|
|
|
|
87
|
9
|
|
|
|
|
24
|
my @masks; |
88
|
|
|
|
|
|
|
|
89
|
|
|
|
|
|
|
# define bin bounds based on equal count (in terms of sites) |
90
|
9
|
100
|
|
|
|
28
|
if ($percentile) { |
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
# create rates-sorted index of sites (from slow to fast) |
93
|
|
|
|
|
|
|
my @index = sort { |
94
|
5
|
|
|
|
|
198
|
$self->state_at($a) <=> $self->state_at($b) |
|
239338
|
|
|
|
|
6838126
|
|
95
|
|
|
|
|
|
|
} 0..$self->mask_len-1; |
96
|
|
|
|
|
|
|
|
97
|
|
|
|
|
|
|
# optionally reverse index: higher values mean slower rates (TIGER) |
98
|
5
|
100
|
|
|
|
493
|
@index = reverse @index if $descending; |
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
# compute masks from index slices |
101
|
5
|
|
|
|
|
56
|
my $step = ceil( @index / $bin_n ); |
102
|
|
|
|
|
|
|
### $step |
103
|
5
|
|
|
|
|
36
|
for my $i (0..$bin_n-1) { |
104
|
32
|
100
|
|
|
|
148
|
my $min = $cumulative ? 0 : ($i * $step); |
105
|
32
|
|
|
|
|
98
|
my $max = (($i+1) * $step - 1); |
106
|
32
|
100
|
|
|
|
119
|
$max = $#index if $max > $#index; |
107
|
|
|
|
|
|
|
### min: $min |
108
|
|
|
|
|
|
|
### max: $max |
109
|
|
|
|
|
|
|
### rates: map { $self->state_at($_) } @index[ $min..$max ] |
110
|
32
|
|
|
|
|
1276
|
push @masks, SeqMask->custom_mask( |
111
|
|
|
|
|
|
|
$self->mask_len, [ @index[ $min..$max ] ] |
112
|
|
|
|
|
|
|
); |
113
|
|
|
|
|
|
|
|
114
|
|
|
|
|
|
|
} |
115
|
|
|
|
|
|
|
} |
116
|
|
|
|
|
|
|
|
117
|
|
|
|
|
|
|
# define bin bounds based on equal width (in terms of rates) |
118
|
|
|
|
|
|
|
else { |
119
|
4
|
|
|
|
|
10
|
my @bounds; |
120
|
|
|
|
|
|
|
|
121
|
|
|
|
|
|
|
# compute bin bounds from rate range |
122
|
4
|
|
|
|
|
20
|
my $min = $self->min_rate; |
123
|
4
|
|
|
|
|
22
|
my $max = $self->max_rate; |
124
|
4
|
|
|
|
|
19
|
my $step = ($max - $min) / $bin_n; |
125
|
4
|
|
|
|
|
21
|
for (my $i = $min + $step; $i <= $max; $i += $step) { |
126
|
40
|
|
|
|
|
79
|
push @bounds, $i; |
127
|
|
|
|
|
|
|
} |
128
|
|
|
|
|
|
|
# Note: did try to use Statistics::Descriptive with no luck |
129
|
|
|
|
|
|
|
|
130
|
|
|
|
|
|
|
# add lower bound for first bin... |
131
|
|
|
|
|
|
|
# ... and small delta to first/last bins for catching min/max values |
132
|
4
|
|
|
|
|
15
|
unshift @bounds, $min - $DELTA; |
133
|
4
|
|
|
|
|
10
|
$bounds[-1] += $DELTA; |
134
|
|
|
|
|
|
|
|
135
|
|
|
|
|
|
|
# optionally reverse bounds: higher values mean slower rates (TIGER) |
136
|
4
|
50
|
|
|
|
17
|
@bounds = reverse @bounds if $descending; |
137
|
|
|
|
|
|
|
### @bounds |
138
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
# derive masks from bin bounds |
140
|
4
|
|
|
|
|
20
|
for my $i (1..$#bounds) { |
141
|
40
|
100
|
|
|
|
447
|
push @masks, $self->rates_mask( |
142
|
|
|
|
|
|
|
$cumulative ? $bounds[0] : $bounds[$i - 1], # bin min |
143
|
|
|
|
|
|
|
$bounds[$i] # bin max |
144
|
|
|
|
|
|
|
); |
145
|
|
|
|
|
|
|
} |
146
|
|
|
|
|
|
|
} |
147
|
|
|
|
|
|
|
|
148
|
9
|
|
|
|
|
110
|
return @masks; |
149
|
|
|
|
|
|
|
} |
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
|
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
sub rates_mask { |
154
|
50
|
|
|
50
|
1
|
242
|
my $self = shift; |
155
|
50
|
|
|
|
|
120
|
my $min = shift; |
156
|
50
|
|
|
|
|
99
|
my $max = shift; |
157
|
|
|
|
|
|
|
|
158
|
|
|
|
|
|
|
# ensure min and max are correctly ordered |
159
|
50
|
|
|
|
|
1022
|
($min, $max) = sort ($min, $max); |
160
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
# filter out sites not within (bin) bounds |
162
|
|
|
|
|
|
|
my @mask = map { |
163
|
50
|
|
|
|
|
3547
|
( |
164
|
45270
|
|
100
|
|
|
1322052
|
$self->state_at($_) > $min |
165
|
|
|
|
|
|
|
&& $self->state_at($_) <= $max |
166
|
|
|
|
|
|
|
) + 0 # ensure proper numeric context (0 or 1) |
167
|
|
|
|
|
|
|
} 0..$self->mask_len-1; |
168
|
|
|
|
|
|
|
|
169
|
50
|
|
|
|
|
5125
|
return SeqMask->new( mask => \@mask ); |
170
|
|
|
|
|
|
|
} |
171
|
|
|
|
|
|
|
|
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
# I/O methods |
174
|
|
|
|
|
|
|
|
175
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
# TODO: avoid duplicating code with SeqMask class |
177
|
|
|
|
|
|
|
|
178
|
|
|
|
|
|
|
# TIGER format: |
179
|
|
|
|
|
|
|
# 0.651325757576 |
180
|
|
|
|
|
|
|
# 0.633143939394 |
181
|
|
|
|
|
|
|
# 0.488257575758 |
182
|
|
|
|
|
|
|
# 1.0 |
183
|
|
|
|
|
|
|
# ... |
184
|
|
|
|
|
|
|
|
185
|
|
|
|
|
|
|
# PhyloBayes (MPI) format: |
186
|
|
|
|
|
|
|
# 0 15.5174 |
187
|
|
|
|
|
|
|
# 1 16.4429 |
188
|
|
|
|
|
|
|
# 2 18.664 |
189
|
|
|
|
|
|
|
# 3 30.1682 |
190
|
|
|
|
|
|
|
# ... |
191
|
|
|
|
|
|
|
|
192
|
|
|
|
|
|
|
# PhyloBayes (serial) format: |
193
|
|
|
|
|
|
|
# site rate std error |
194
|
|
|
|
|
|
|
# |
195
|
|
|
|
|
|
|
# 1 5.86191 0.207244 |
196
|
|
|
|
|
|
|
# 2 7.35053 0.350185 |
197
|
|
|
|
|
|
|
# 3 6.09158 0.24426 |
198
|
|
|
|
|
|
|
# ... |
199
|
|
|
|
|
|
|
|
200
|
|
|
|
|
|
|
# IQ-TREE format: |
201
|
|
|
|
|
|
|
# # Site-specific subtitution rates determined by empirical Bayesian method |
202
|
|
|
|
|
|
|
# # This file can be read in MS Excel or in R with command: |
203
|
|
|
|
|
|
|
# # tab=read.table('replica-3-concat-modified.fasta.rate',header=TRUE) |
204
|
|
|
|
|
|
|
# # Columns are tab-separated with following meaning: |
205
|
|
|
|
|
|
|
# # Site: Alignment site ID |
206
|
|
|
|
|
|
|
# # Rate: Posterior mean site rate weighted by posterior probability |
207
|
|
|
|
|
|
|
# # Cat: Category with highest posterior (0=invariable, 1=slow, etc) |
208
|
|
|
|
|
|
|
# # C_Rate: Corresponding rate of highest category |
209
|
|
|
|
|
|
|
# Site Rate Cat C_Rate |
210
|
|
|
|
|
|
|
# 1 0.31154 1 0.30935 |
211
|
|
|
|
|
|
|
# 2 1.03326 4 1.91612 |
212
|
|
|
|
|
|
|
# 3 1.39822 4 1.91612 |
213
|
|
|
|
|
|
|
# ... |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
const my %COLUMN_FOR => ( |
216
|
|
|
|
|
|
|
1 => 0, # TIGER |
217
|
|
|
|
|
|
|
2 => 1, # PhyloBayes (MPI) |
218
|
|
|
|
|
|
|
3 => 1, # PhyloBayes (serial) |
219
|
|
|
|
|
|
|
4 => 1, # IQ-TREE |
220
|
|
|
|
|
|
|
); |
221
|
|
|
|
|
|
|
|
222
|
|
|
|
|
|
|
sub load { |
223
|
12
|
|
|
12
|
1
|
2275
|
my $class = shift; |
224
|
12
|
|
|
|
|
31
|
my $infile = shift; |
225
|
|
|
|
|
|
|
|
226
|
12
|
|
|
|
|
79
|
open my $in, '<', $infile; |
227
|
|
|
|
|
|
|
|
228
|
12
|
|
|
|
|
6177
|
my $mask = $class->new(); |
229
|
12
|
|
|
|
|
30
|
my $col; |
230
|
|
|
|
|
|
|
|
231
|
|
|
|
|
|
|
LINE: |
232
|
12
|
|
|
|
|
6295
|
while (my $line = <$in>) { |
233
|
71734
|
|
|
|
|
115130
|
chomp $line; |
234
|
|
|
|
|
|
|
|
235
|
|
|
|
|
|
|
# skip empty lines, header line and process comment lines |
236
|
71734
|
100
|
100
|
|
|
457097
|
next LINE if $line =~ $EMPTY_LINE |
|
|
|
100
|
|
|
|
|
237
|
|
|
|
|
|
|
|| $line =~ m/^site/xmsi # PhyloBayes/IQ-TREE |
238
|
|
|
|
|
|
|
|| $mask->is_comment($line); |
239
|
|
|
|
|
|
|
|
240
|
|
|
|
|
|
|
# try to split line on whitespace |
241
|
71723
|
|
|
|
|
239628
|
my @fields = split /\s+/xms, $line; |
242
|
71723
|
|
100
|
|
|
145228
|
$col //= $COLUMN_FOR{ scalar @fields }; |
243
|
|
|
|
|
|
|
|
244
|
|
|
|
|
|
|
# store either first or second field depending on split outcome |
245
|
|
|
|
|
|
|
# Note: we check this outcome only once for efficiency |
246
|
71723
|
|
|
|
|
2512269
|
$mask->add_state( $fields[$col] ); |
247
|
|
|
|
|
|
|
} |
248
|
|
|
|
|
|
|
|
249
|
12
|
|
|
|
|
437
|
return $mask; |
250
|
|
|
|
|
|
|
} |
251
|
|
|
|
|
|
|
|
252
|
|
|
|
|
|
|
__PACKAGE__->meta->make_immutable; |
253
|
|
|
|
|
|
|
1; |
254
|
|
|
|
|
|
|
|
255
|
|
|
|
|
|
|
__END__ |
256
|
|
|
|
|
|
|
|
257
|
|
|
|
|
|
|
=pod |
258
|
|
|
|
|
|
|
|
259
|
|
|
|
|
|
|
=head1 NAME |
260
|
|
|
|
|
|
|
|
261
|
|
|
|
|
|
|
Bio::MUST::Core::SeqMask::Rates - Evolutionary rates for sequence sites |
262
|
|
|
|
|
|
|
|
263
|
|
|
|
|
|
|
=head1 VERSION |
264
|
|
|
|
|
|
|
|
265
|
|
|
|
|
|
|
version 0.212530 |
266
|
|
|
|
|
|
|
|
267
|
|
|
|
|
|
|
=head1 SYNOPSIS |
268
|
|
|
|
|
|
|
|
269
|
|
|
|
|
|
|
# TODO |
270
|
|
|
|
|
|
|
|
271
|
|
|
|
|
|
|
=head1 DESCRIPTION |
272
|
|
|
|
|
|
|
|
273
|
|
|
|
|
|
|
# TODO |
274
|
|
|
|
|
|
|
|
275
|
|
|
|
|
|
|
=head1 METHODS |
276
|
|
|
|
|
|
|
|
277
|
|
|
|
|
|
|
=head2 min_rate |
278
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
=head2 max_rate |
280
|
|
|
|
|
|
|
|
281
|
|
|
|
|
|
|
=head2 delta_rates |
282
|
|
|
|
|
|
|
|
283
|
|
|
|
|
|
|
=head2 bin_rates_masks |
284
|
|
|
|
|
|
|
|
285
|
|
|
|
|
|
|
=head2 rates_mask |
286
|
|
|
|
|
|
|
|
287
|
|
|
|
|
|
|
=head2 load |
288
|
|
|
|
|
|
|
|
289
|
|
|
|
|
|
|
=head1 AUTHOR |
290
|
|
|
|
|
|
|
|
291
|
|
|
|
|
|
|
Denis BAURAIN <denis.baurain@uliege.be> |
292
|
|
|
|
|
|
|
|
293
|
|
|
|
|
|
|
=head1 COPYRIGHT AND LICENSE |
294
|
|
|
|
|
|
|
|
295
|
|
|
|
|
|
|
This software is copyright (c) 2013 by University of Liege / Unit of Eukaryotic Phylogenomics / Denis BAURAIN. |
296
|
|
|
|
|
|
|
|
297
|
|
|
|
|
|
|
This is free software; you can redistribute it and/or modify it under |
298
|
|
|
|
|
|
|
the same terms as the Perl 5 programming language system itself. |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
=cut |