line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package App::Anchr::Common; |
2
|
25
|
|
|
25
|
|
50847
|
use strict; |
|
25
|
|
|
|
|
60
|
|
|
25
|
|
|
|
|
745
|
|
3
|
25
|
|
|
25
|
|
121
|
use warnings; |
|
25
|
|
|
|
|
51
|
|
|
25
|
|
|
|
|
666
|
|
4
|
25
|
|
|
25
|
|
1044
|
use autodie; |
|
25
|
|
|
|
|
24622
|
|
|
25
|
|
|
|
|
121
|
|
5
|
|
|
|
|
|
|
|
6
|
25
|
|
|
25
|
|
114879
|
use 5.010001; |
|
25
|
|
|
|
|
383
|
|
7
|
|
|
|
|
|
|
|
8
|
25
|
|
|
25
|
|
132
|
use Carp qw(); |
|
25
|
|
|
|
|
49
|
|
|
25
|
|
|
|
|
371
|
|
9
|
25
|
|
|
25
|
|
8550
|
use File::ShareDir qw(); |
|
25
|
|
|
|
|
119471
|
|
|
25
|
|
|
|
|
596
|
|
10
|
25
|
|
|
25
|
|
12514
|
use Graph; |
|
25
|
|
|
|
|
1986726
|
|
|
25
|
|
|
|
|
928
|
|
11
|
25
|
|
|
25
|
|
1797
|
use GraphViz; |
|
0
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
use IO::Zlib; |
13
|
|
|
|
|
|
|
use IPC::Cmd qw(); |
14
|
|
|
|
|
|
|
use JSON qw(); |
15
|
|
|
|
|
|
|
use List::Util; |
16
|
|
|
|
|
|
|
use Path::Tiny qw(); |
17
|
|
|
|
|
|
|
use Statistics::Descriptive qw(); |
18
|
|
|
|
|
|
|
use Template; |
19
|
|
|
|
|
|
|
use YAML::Syck qw(); |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
use AlignDB::IntSpan; |
22
|
|
|
|
|
|
|
use AlignDB::Stopwatch; |
23
|
|
|
|
|
|
|
use App::RL::Common; |
24
|
|
|
|
|
|
|
use App::Fasops::Common; |
25
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
sub get_len_from_header { |
27
|
|
|
|
|
|
|
my $fa_fn = shift; |
28
|
|
|
|
|
|
|
|
29
|
|
|
|
|
|
|
my %len_of; |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
my $fa_fh; |
32
|
|
|
|
|
|
|
if ( lc $fa_fn eq 'stdin' ) { |
33
|
|
|
|
|
|
|
$fa_fh = *STDIN{IO}; |
34
|
|
|
|
|
|
|
} |
35
|
|
|
|
|
|
|
else { |
36
|
|
|
|
|
|
|
open $fa_fh, "<", $fa_fn; |
37
|
|
|
|
|
|
|
} |
38
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
while ( my $line = <$fa_fh> ) { |
40
|
|
|
|
|
|
|
if ( substr( $line, 0, 1 ) eq ">" ) { |
41
|
|
|
|
|
|
|
if ( $line =~ /\/(\d+)\/\d+_(\d+)/ ) { |
42
|
|
|
|
|
|
|
$len_of{$1} = $2; |
43
|
|
|
|
|
|
|
} |
44
|
|
|
|
|
|
|
} |
45
|
|
|
|
|
|
|
} |
46
|
|
|
|
|
|
|
|
47
|
|
|
|
|
|
|
close $fa_fh; |
48
|
|
|
|
|
|
|
|
49
|
|
|
|
|
|
|
return \%len_of; |
50
|
|
|
|
|
|
|
} |
51
|
|
|
|
|
|
|
|
52
|
|
|
|
|
|
|
sub get_replaces { |
53
|
|
|
|
|
|
|
my $fn = shift; |
54
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
my $replace_of = {}; |
56
|
|
|
|
|
|
|
my @lines = Path::Tiny::path($fn)->lines( { chomp => 1 } ); |
57
|
|
|
|
|
|
|
|
58
|
|
|
|
|
|
|
for my $line (@lines) { |
59
|
|
|
|
|
|
|
my @fields = split /\t/, $line; |
60
|
|
|
|
|
|
|
if ( @fields == 2 ) { |
61
|
|
|
|
|
|
|
if ( $fields[0] =~ /\/(\d+)\/\d+_\d+/ ) { |
62
|
|
|
|
|
|
|
$replace_of->{$1} = $fields[1]; |
63
|
|
|
|
|
|
|
} |
64
|
|
|
|
|
|
|
} |
65
|
|
|
|
|
|
|
} |
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
return $replace_of; |
68
|
|
|
|
|
|
|
} |
69
|
|
|
|
|
|
|
|
70
|
|
|
|
|
|
|
sub get_replaces2 { |
71
|
|
|
|
|
|
|
my $fn = shift; |
72
|
|
|
|
|
|
|
my $opt = shift; |
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
my $replace_of = {}; |
75
|
|
|
|
|
|
|
my @lines = Path::Tiny::path($fn)->lines( { chomp => 1 } ); |
76
|
|
|
|
|
|
|
|
77
|
|
|
|
|
|
|
for my $line (@lines) { |
78
|
|
|
|
|
|
|
my @fields = split /\t/, $line; |
79
|
|
|
|
|
|
|
if ( @fields == 2 ) { |
80
|
|
|
|
|
|
|
if ( defined $opt and ref $opt eq "HASH" and $opt->{reverse} ) { |
81
|
|
|
|
|
|
|
$replace_of->{ $fields[1] } = $fields[0]; |
82
|
|
|
|
|
|
|
} |
83
|
|
|
|
|
|
|
else { |
84
|
|
|
|
|
|
|
$replace_of->{ $fields[0] } = $fields[1]; |
85
|
|
|
|
|
|
|
} |
86
|
|
|
|
|
|
|
} |
87
|
|
|
|
|
|
|
} |
88
|
|
|
|
|
|
|
|
89
|
|
|
|
|
|
|
return $replace_of; |
90
|
|
|
|
|
|
|
} |
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
sub exec_cmd { |
93
|
|
|
|
|
|
|
my $cmd = shift; |
94
|
|
|
|
|
|
|
my $opt = shift; |
95
|
|
|
|
|
|
|
|
96
|
|
|
|
|
|
|
if ( defined $opt and ref $opt eq "HASH" and $opt->{verbose} ) { |
97
|
|
|
|
|
|
|
print STDERR "CMD: ", $cmd, "\n"; |
98
|
|
|
|
|
|
|
} |
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
system $cmd; |
101
|
|
|
|
|
|
|
} |
102
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
sub parse_ovlp_line { |
104
|
|
|
|
|
|
|
my $line = shift; |
105
|
|
|
|
|
|
|
|
106
|
|
|
|
|
|
|
chomp $line; |
107
|
|
|
|
|
|
|
my @fields = split "\t", $line; |
108
|
|
|
|
|
|
|
|
109
|
|
|
|
|
|
|
my $info = { |
110
|
|
|
|
|
|
|
f_id => $fields[0], |
111
|
|
|
|
|
|
|
g_id => $fields[1], |
112
|
|
|
|
|
|
|
ovlp_len => $fields[2], |
113
|
|
|
|
|
|
|
ovlp_idt => $fields[3], |
114
|
|
|
|
|
|
|
f_strand => $fields[4], |
115
|
|
|
|
|
|
|
f_B => $fields[5], |
116
|
|
|
|
|
|
|
f_E => $fields[6], |
117
|
|
|
|
|
|
|
f_len => $fields[7], |
118
|
|
|
|
|
|
|
g_strand => $fields[8], |
119
|
|
|
|
|
|
|
g_B => $fields[9], |
120
|
|
|
|
|
|
|
g_E => $fields[10], |
121
|
|
|
|
|
|
|
g_len => $fields[11], |
122
|
|
|
|
|
|
|
contained => $fields[12], |
123
|
|
|
|
|
|
|
}; |
124
|
|
|
|
|
|
|
|
125
|
|
|
|
|
|
|
return $info; |
126
|
|
|
|
|
|
|
} |
127
|
|
|
|
|
|
|
|
128
|
|
|
|
|
|
|
sub create_ovlp_line { |
129
|
|
|
|
|
|
|
my $info = shift; |
130
|
|
|
|
|
|
|
|
131
|
|
|
|
|
|
|
my $line = join "\t", |
132
|
|
|
|
|
|
|
$info->{f_id}, $info->{g_id}, $info->{ovlp_len}, $info->{ovlp_idt}, |
133
|
|
|
|
|
|
|
$info->{f_strand}, $info->{f_B}, $info->{f_E}, $info->{f_len}, |
134
|
|
|
|
|
|
|
$info->{g_strand}, $info->{g_B}, $info->{g_E}, $info->{g_len}, $info->{contained}; |
135
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
return $line; |
137
|
|
|
|
|
|
|
} |
138
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
sub parse_paf_line { |
140
|
|
|
|
|
|
|
my $line = shift; |
141
|
|
|
|
|
|
|
|
142
|
|
|
|
|
|
|
chomp $line; |
143
|
|
|
|
|
|
|
my @fields = split "\t", $line; |
144
|
|
|
|
|
|
|
|
145
|
|
|
|
|
|
|
my $info = { |
146
|
|
|
|
|
|
|
f_id => $fields[0], |
147
|
|
|
|
|
|
|
g_id => $fields[5], |
148
|
|
|
|
|
|
|
ovlp_len => $fields[10], |
149
|
|
|
|
|
|
|
ovlp_idt => sprintf( "%.3f", $fields[9] / $fields[10] ), |
150
|
|
|
|
|
|
|
f_strand => 0, |
151
|
|
|
|
|
|
|
f_B => $fields[2] + 1, |
152
|
|
|
|
|
|
|
f_E => $fields[3] + 1, |
153
|
|
|
|
|
|
|
f_len => $fields[1], |
154
|
|
|
|
|
|
|
$fields[4] eq "+" |
155
|
|
|
|
|
|
|
? ( g_strand => 0, |
156
|
|
|
|
|
|
|
g_B => $fields[7] + 1, |
157
|
|
|
|
|
|
|
g_E => $fields[8] + 1, |
158
|
|
|
|
|
|
|
) |
159
|
|
|
|
|
|
|
: ( g_strand => 1, |
160
|
|
|
|
|
|
|
g_B => $fields[8] + 1, |
161
|
|
|
|
|
|
|
g_E => $fields[7] + 1, |
162
|
|
|
|
|
|
|
), |
163
|
|
|
|
|
|
|
g_len => $fields[6], |
164
|
|
|
|
|
|
|
contained => 'overlap', |
165
|
|
|
|
|
|
|
}; |
166
|
|
|
|
|
|
|
|
167
|
|
|
|
|
|
|
for my $key (qw(f_B f_E g_B g_E)) { |
168
|
|
|
|
|
|
|
if ( $info->{$key} == 1 ) { |
169
|
|
|
|
|
|
|
$info->{$key} = 0; |
170
|
|
|
|
|
|
|
} |
171
|
|
|
|
|
|
|
} |
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
return $info; |
174
|
|
|
|
|
|
|
} |
175
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
sub beg_end { |
177
|
|
|
|
|
|
|
my $beg = shift; |
178
|
|
|
|
|
|
|
my $end = shift; |
179
|
|
|
|
|
|
|
|
180
|
|
|
|
|
|
|
if ( $beg > $end ) { |
181
|
|
|
|
|
|
|
( $beg, $end ) = ( $end, $beg ); |
182
|
|
|
|
|
|
|
} |
183
|
|
|
|
|
|
|
|
184
|
|
|
|
|
|
|
if ( $beg == 0 ) { |
185
|
|
|
|
|
|
|
$beg = 1; |
186
|
|
|
|
|
|
|
} |
187
|
|
|
|
|
|
|
|
188
|
|
|
|
|
|
|
return ( $beg, $end ); |
189
|
|
|
|
|
|
|
} |
190
|
|
|
|
|
|
|
|
191
|
|
|
|
|
|
|
sub bump_coverage { |
192
|
|
|
|
|
|
|
my $tier_of = shift; |
193
|
|
|
|
|
|
|
my $beg = shift; |
194
|
|
|
|
|
|
|
my $end = shift; |
195
|
|
|
|
|
|
|
my $coverage = shift || 2; |
196
|
|
|
|
|
|
|
|
197
|
|
|
|
|
|
|
return if $tier_of->{$coverage}->equals( $tier_of->{all} ); |
198
|
|
|
|
|
|
|
|
199
|
|
|
|
|
|
|
( $beg, $end ) = beg_end( $beg, $end ); |
200
|
|
|
|
|
|
|
|
201
|
|
|
|
|
|
|
my $new_set = AlignDB::IntSpan->new->add_pair( $beg, $end ); |
202
|
|
|
|
|
|
|
for my $i ( 1 .. $coverage ) { |
203
|
|
|
|
|
|
|
my $i_set = $tier_of->{$i}->intersect($new_set); |
204
|
|
|
|
|
|
|
$tier_of->{$i}->add($new_set); |
205
|
|
|
|
|
|
|
|
206
|
|
|
|
|
|
|
my $j = $i + 1; |
207
|
|
|
|
|
|
|
last if $j > $coverage; |
208
|
|
|
|
|
|
|
$new_set = $i_set->copy; |
209
|
|
|
|
|
|
|
} |
210
|
|
|
|
|
|
|
} |
211
|
|
|
|
|
|
|
|
212
|
|
|
|
|
|
|
sub serial2name { |
213
|
|
|
|
|
|
|
my $dazz_db = shift; |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
my $serials = shift; |
216
|
|
|
|
|
|
|
|
217
|
|
|
|
|
|
|
my $cmd = sprintf "DBshow -n %s %s ", $dazz_db, join( " ", @{$serials} ); |
218
|
|
|
|
|
|
|
my @lines = map { $_ =~ s/^>//; $_; } grep {defined} split /\n/, `$cmd`; |
219
|
|
|
|
|
|
|
|
220
|
|
|
|
|
|
|
my %name_of; |
221
|
|
|
|
|
|
|
for my $i ( 0 .. $#lines ) { |
222
|
|
|
|
|
|
|
$name_of{ $serials->[$i] } = $lines[$i]; |
223
|
|
|
|
|
|
|
} |
224
|
|
|
|
|
|
|
|
225
|
|
|
|
|
|
|
return \%name_of; |
226
|
|
|
|
|
|
|
} |
227
|
|
|
|
|
|
|
|
228
|
|
|
|
|
|
|
sub judge_distance { |
229
|
|
|
|
|
|
|
my $d_ref = shift; |
230
|
|
|
|
|
|
|
my $max_dis = shift || 5000; |
231
|
|
|
|
|
|
|
|
232
|
|
|
|
|
|
|
return 0 unless defined $d_ref; |
233
|
|
|
|
|
|
|
|
234
|
|
|
|
|
|
|
my $sum = 0; |
235
|
|
|
|
|
|
|
my $min = $d_ref->[0]; |
236
|
|
|
|
|
|
|
my $max = $min; |
237
|
|
|
|
|
|
|
for my $d ( @{$d_ref} ) { |
238
|
|
|
|
|
|
|
$sum += $d; |
239
|
|
|
|
|
|
|
if ( $d < $min ) { $min = $d; } |
240
|
|
|
|
|
|
|
if ( $d > $max ) { $max = $d; } |
241
|
|
|
|
|
|
|
} |
242
|
|
|
|
|
|
|
my $avg = $sum / scalar( @{$d_ref} ); |
243
|
|
|
|
|
|
|
return 0 if $avg > $max_dis; |
244
|
|
|
|
|
|
|
|
245
|
|
|
|
|
|
|
# max k-mer is 127. |
246
|
|
|
|
|
|
|
# For k-unitigs, overlaps are less than k-mer |
247
|
|
|
|
|
|
|
my $v = $max - $min; |
248
|
|
|
|
|
|
|
if ( $v < 20 or abs( $v / $avg ) < 0.2 ) { |
249
|
|
|
|
|
|
|
return 1; |
250
|
|
|
|
|
|
|
} |
251
|
|
|
|
|
|
|
else { |
252
|
|
|
|
|
|
|
return 0; |
253
|
|
|
|
|
|
|
} |
254
|
|
|
|
|
|
|
} |
255
|
|
|
|
|
|
|
|
256
|
|
|
|
|
|
|
sub g2gv { |
257
|
|
|
|
|
|
|
|
258
|
|
|
|
|
|
|
#@type Graph |
259
|
|
|
|
|
|
|
my $g = shift; |
260
|
|
|
|
|
|
|
my $fn = shift; |
261
|
|
|
|
|
|
|
|
262
|
|
|
|
|
|
|
my $gv = GraphViz->new( directed => 1 ); |
263
|
|
|
|
|
|
|
|
264
|
|
|
|
|
|
|
for my $v ( $g->vertices ) { |
265
|
|
|
|
|
|
|
$gv->add_node($v); |
266
|
|
|
|
|
|
|
} |
267
|
|
|
|
|
|
|
|
268
|
|
|
|
|
|
|
for my $e ( $g->edges ) { |
269
|
|
|
|
|
|
|
if ( $g->has_edge_weight( @{$e} ) ) { |
270
|
|
|
|
|
|
|
$gv->add_edge( @{$e}, label => $g->get_edge_weight( @{$e} ) ); |
271
|
|
|
|
|
|
|
} |
272
|
|
|
|
|
|
|
else { |
273
|
|
|
|
|
|
|
$gv->add_edge( @{$e} ); |
274
|
|
|
|
|
|
|
} |
275
|
|
|
|
|
|
|
} |
276
|
|
|
|
|
|
|
|
277
|
|
|
|
|
|
|
Path::Tiny::path($fn)->spew_raw( $gv->as_png ); |
278
|
|
|
|
|
|
|
} |
279
|
|
|
|
|
|
|
|
280
|
|
|
|
|
|
|
sub g2gv0 { |
281
|
|
|
|
|
|
|
|
282
|
|
|
|
|
|
|
#@type Graph |
283
|
|
|
|
|
|
|
my $g = shift; |
284
|
|
|
|
|
|
|
my $fn = shift; |
285
|
|
|
|
|
|
|
|
286
|
|
|
|
|
|
|
my $gv = GraphViz->new( directed => 0 ); |
287
|
|
|
|
|
|
|
|
288
|
|
|
|
|
|
|
for my $v ( $g->vertices ) { |
289
|
|
|
|
|
|
|
$gv->add_node($v); |
290
|
|
|
|
|
|
|
} |
291
|
|
|
|
|
|
|
|
292
|
|
|
|
|
|
|
for my $e ( $g->edges ) { |
293
|
|
|
|
|
|
|
$gv->add_edge( @{$e} ); |
294
|
|
|
|
|
|
|
} |
295
|
|
|
|
|
|
|
|
296
|
|
|
|
|
|
|
Path::Tiny::path($fn)->spew_raw( $gv->as_png ); |
297
|
|
|
|
|
|
|
} |
298
|
|
|
|
|
|
|
|
299
|
|
|
|
|
|
|
sub transitive_reduction { |
300
|
|
|
|
|
|
|
|
301
|
|
|
|
|
|
|
#@type Graph |
302
|
|
|
|
|
|
|
my $g = shift; |
303
|
|
|
|
|
|
|
|
304
|
|
|
|
|
|
|
my $count = 0; |
305
|
|
|
|
|
|
|
my $prev_count; |
306
|
|
|
|
|
|
|
while (1) { |
307
|
|
|
|
|
|
|
last if defined $prev_count and $prev_count == $count; |
308
|
|
|
|
|
|
|
$prev_count = $count; |
309
|
|
|
|
|
|
|
|
310
|
|
|
|
|
|
|
for my $v ( $g->vertices ) { |
311
|
|
|
|
|
|
|
next if $g->out_degree($v) < 2; |
312
|
|
|
|
|
|
|
|
313
|
|
|
|
|
|
|
my @s = sort { $a cmp $b } $g->successors($v); |
314
|
|
|
|
|
|
|
for my $i ( 0 .. $#s ) { |
315
|
|
|
|
|
|
|
for my $j ( 0 .. $#s ) { |
316
|
|
|
|
|
|
|
next if $i == $j; |
317
|
|
|
|
|
|
|
if ( $g->is_reachable( $s[$i], $s[$j] ) ) { |
318
|
|
|
|
|
|
|
$g->delete_edge( $v, $s[$j] ); |
319
|
|
|
|
|
|
|
|
320
|
|
|
|
|
|
|
$count++; |
321
|
|
|
|
|
|
|
} |
322
|
|
|
|
|
|
|
} |
323
|
|
|
|
|
|
|
} |
324
|
|
|
|
|
|
|
} |
325
|
|
|
|
|
|
|
} |
326
|
|
|
|
|
|
|
|
327
|
|
|
|
|
|
|
return $count; |
328
|
|
|
|
|
|
|
} |
329
|
|
|
|
|
|
|
|
330
|
|
|
|
|
|
|
sub poa_consensus { |
331
|
|
|
|
|
|
|
my $seq_refs = shift; |
332
|
|
|
|
|
|
|
|
333
|
|
|
|
|
|
|
my $aln_prog = "poa"; |
334
|
|
|
|
|
|
|
|
335
|
|
|
|
|
|
|
# temp in and out |
336
|
|
|
|
|
|
|
my $temp_in = Path::Tiny->tempfile("seq_in_XXXXXXXX"); |
337
|
|
|
|
|
|
|
my $temp_out = Path::Tiny->tempfile("seq_out_XXXXXXXX"); |
338
|
|
|
|
|
|
|
|
339
|
|
|
|
|
|
|
# msa may change the order of sequences |
340
|
|
|
|
|
|
|
my @indexes = 0 .. scalar( @{$seq_refs} - 1 ); |
341
|
|
|
|
|
|
|
{ |
342
|
|
|
|
|
|
|
my $fh = $temp_in->openw(); |
343
|
|
|
|
|
|
|
for my $i (@indexes) { |
344
|
|
|
|
|
|
|
printf {$fh} ">seq_%d\n", $i; |
345
|
|
|
|
|
|
|
printf {$fh} "%s\n", $seq_refs->[$i]; |
346
|
|
|
|
|
|
|
} |
347
|
|
|
|
|
|
|
close $fh; |
348
|
|
|
|
|
|
|
} |
349
|
|
|
|
|
|
|
|
350
|
|
|
|
|
|
|
my @args; |
351
|
|
|
|
|
|
|
{ |
352
|
|
|
|
|
|
|
push @args, "-hb"; |
353
|
|
|
|
|
|
|
push @args, "-read_fasta " . $temp_in->absolute->stringify; |
354
|
|
|
|
|
|
|
push @args, "-clustal " . $temp_out->absolute->stringify; |
355
|
|
|
|
|
|
|
push @args, File::ShareDir::dist_file( 'App-Anchr', 'poa-blosum80.mat' ); |
356
|
|
|
|
|
|
|
} |
357
|
|
|
|
|
|
|
|
358
|
|
|
|
|
|
|
my $cmd_line = join " ", ( $aln_prog, @args ); |
359
|
|
|
|
|
|
|
my $ok = IPC::Cmd::run( command => $cmd_line ); |
360
|
|
|
|
|
|
|
|
361
|
|
|
|
|
|
|
if ( !$ok ) { |
362
|
|
|
|
|
|
|
Carp::confess("Calling [$aln_prog] failed\n"); |
363
|
|
|
|
|
|
|
} |
364
|
|
|
|
|
|
|
|
365
|
|
|
|
|
|
|
my $consensus = join "", grep {/^CONSENS0/} $temp_out->lines( { chomp => 1, } ); |
366
|
|
|
|
|
|
|
$consensus =~ s/CONSENS0//g; |
367
|
|
|
|
|
|
|
$consensus =~ s/\s//g; |
368
|
|
|
|
|
|
|
$consensus =~ s/-//g; |
369
|
|
|
|
|
|
|
|
370
|
|
|
|
|
|
|
return $consensus; |
371
|
|
|
|
|
|
|
} |
372
|
|
|
|
|
|
|
|
373
|
|
|
|
|
|
|
# https://metacpan.org/source/GSULLIVAN/String-LCSS-1.00/lib/String/LCSS.pm |
374
|
|
|
|
|
|
|
# `undef` is returned if the susbstring length is one char or less. |
375
|
|
|
|
|
|
|
# In scalar context, returns the substring. |
376
|
|
|
|
|
|
|
# In array context, returns the index of the match root in the two args. |
377
|
|
|
|
|
|
|
sub lcss { |
378
|
|
|
|
|
|
|
my $solns0 = ( _lcss( $_[0], $_[1] ) )[0]; |
379
|
|
|
|
|
|
|
return unless $solns0; |
380
|
|
|
|
|
|
|
my @match = @{$solns0}; |
381
|
|
|
|
|
|
|
return if length $match[0] == 1; |
382
|
|
|
|
|
|
|
return wantarray ? @match : $match[0]; |
383
|
|
|
|
|
|
|
} |
384
|
|
|
|
|
|
|
|
385
|
|
|
|
|
|
|
sub _lcss { |
386
|
|
|
|
|
|
|
|
387
|
|
|
|
|
|
|
# Return array-of-arrays of longest substrings and indices |
388
|
|
|
|
|
|
|
my ( $r1, $r2 ) = @_; |
389
|
|
|
|
|
|
|
my ( $l1, $l2, ) = ( length $r1, length $r2, ); |
390
|
|
|
|
|
|
|
( $r1, $r2, $l1, $l2, ) = ( $r2, $r1, $l2, $l1, ) if $l1 > $l2; |
391
|
|
|
|
|
|
|
|
392
|
|
|
|
|
|
|
my ( $best, @solns ) = 0; |
393
|
|
|
|
|
|
|
for my $start ( 0 .. $l2 - 1 ) { |
394
|
|
|
|
|
|
|
for my $l ( reverse 1 .. $l1 - $start ) { |
395
|
|
|
|
|
|
|
my $substr = substr( $r1, $start, $l ); |
396
|
|
|
|
|
|
|
my $o = index( $r2, $substr ); |
397
|
|
|
|
|
|
|
next if $o < 0; |
398
|
|
|
|
|
|
|
if ( $l > $best ) { |
399
|
|
|
|
|
|
|
$best = length $substr; |
400
|
|
|
|
|
|
|
@solns = [ $substr, $start, $o ]; |
401
|
|
|
|
|
|
|
} |
402
|
|
|
|
|
|
|
elsif ( $l == $best ) { |
403
|
|
|
|
|
|
|
push @solns, [ $substr, $start, $o ]; |
404
|
|
|
|
|
|
|
} |
405
|
|
|
|
|
|
|
} |
406
|
|
|
|
|
|
|
} |
407
|
|
|
|
|
|
|
return @solns; |
408
|
|
|
|
|
|
|
} |
409
|
|
|
|
|
|
|
|
410
|
|
|
|
|
|
|
1; |