[Bio] / FigKernelPackages / NewGenome.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/NewGenome.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.75 - (view) (download) (as text)

1 : overbeek 1.25 # -*- perl -*-
2 : gdpusch 1.49 ########################################################################
3 : overbeek 1.1 # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 : gdpusch 1.49 ########################################################################
18 : overbeek 1.1
19 :     package NewGenome;
20 :    
21 :     use strict;
22 :     use FIG;
23 :     my $fig = new FIG;
24 :    
25 : overbeek 1.2 use Carp qw(:DEFAULT cluck);
26 : overbeek 1.1 use Data::Dumper;
27 :    
28 : overbeek 1.29 # use Memoize;
29 :     # memoize('search_for_upstream_start');
30 :     # memoize('search_for_upstream_stop');
31 :     # memoize('search_for_downstream_start');
32 : gdpusch 1.75 # memoize('extend_to_downstream_stop');
33 : overbeek 1.26
34 : overbeek 1.1 use vars '$AUTOLOAD';
35 :     sub AUTOLOAD
36 :     {
37 :     my ($self, @args) = @_;
38 :    
39 : overbeek 1.6 confess "\nObject accessed via unknown method $AUTOLOAD,\n"
40 : overbeek 1.1 , "with args: ", join(", ", @args), "\n\n"
41 :     , Dumper($self);
42 :     }
43 :    
44 :     sub DESTROY
45 :     {
46 :     #...Currently, nothing needs to be done...
47 :     }
48 :    
49 :     sub what_methods
50 :     {
51 :     foreach my $symname (sort keys %NewGenome::)
52 :     {
53 :     local *name = $NewGenome::{$symname};
54 :     print STDERR "$symname\n" if defined(&name);
55 :     }
56 :     }
57 :    
58 :     use constant _CONTIG => 0;
59 :     use constant _LOC_END => 1;
60 :     use constant _STRAND => 2;
61 :     use constant _LENGTH => 3;
62 :     use constant _ORF_LEN => 4;
63 :     use constant _ORF_TRANS => 5;
64 :     use constant _ORF_SIMS => 6;
65 :    
66 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
67 :     # This is the constructor. Presumably, $class is 'NewGenome'.
68 :     # usage: my $newG = new NewGenome($dir)
69 :     #-----------------------------------------------------------------------
70 :     sub new
71 :     {
72 : overbeek 1.4 my ($class, $dir, @features) = @_;
73 : overbeek 1.1
74 :     my $newG = {};
75 : overbeek 1.6 bless $newG, $class;
76 :    
77 : overbeek 1.1 $newG->{_dir} = $dir;
78 : overbeek 1.25 ($dir =~ /^(.*\/)?(\d+\.\d+)$/) || die "Skeleton path $dir does not end in a well-formed Org-ID";
79 : overbeek 1.1 my $taxid = $newG->{_taxid} = $2;
80 :    
81 :     (-s "$dir/contigs") || die "$dir/contigs does not exist or has zero size";
82 :     (-s "$dir/GENOME") || die "$dir/GENOME does not exist or has zero size";
83 :     (-s "$dir/PROJECT") || die "$dir/PROJECT does not exist or has zero size";
84 :     (-s "$dir/TAXONOMY") || die "$dir/TAXONOMY does not exist or has zero size";
85 :    
86 :     &FIG::verify_dir("$dir/Features");
87 :    
88 :     $newG->{_genome} = `head -1 $dir/GENOME` or die "Could not read $dir/GENOME";
89 :     chomp $newG->{_genome};
90 :    
91 :     $newG->{_project} = `cat $dir/PROJECT` or die "Could not read $dir/PROJECT";
92 :     chomp $newG->{_project};
93 :    
94 :     $newG->{_taxonomy} = `head -1 $dir/TAXONOMY` or die "Could not read $dir/TAXONOMY";
95 :     chomp $newG->{_taxonomy};
96 :    
97 : overbeek 1.47 if (-s "$dir/GENETIC_CODE") {
98 :     my $code = `head -1 $dir/GENETIC_CODE`;
99 :     chomp $code;
100 :     if ($code =~ m/^(\d+)/o) {
101 : overbeek 1.48 print STDERR "Using genetic code $1\n" if $ENV{VERBOSE};
102 : gdpusch 1.49 $newG->{_genetic_code_number} = $1;
103 :     $newG->{_translation_table} = &FIG::genetic_code($1);
104 : overbeek 1.47 }
105 :     else {
106 :     die "Sorry, cannot handle non-numeric genetic code $code";
107 :     }
108 :     }
109 :     else {
110 :     #...Default to standard code:
111 : overbeek 1.48 print STDERR "Using standard genetic code\n" if $ENV{VERBOSE};
112 : gdpusch 1.49 $newG->{_genetic_code_number} = 11;
113 :     $newG->{_translation_table} = &FIG::standard_genetic_code();
114 : overbeek 1.47 }
115 :    
116 : overbeek 1.26 my $contig_seqs = {};
117 :     my $contig_lens = {};
118 : overbeek 1.1 open(CONTIGS,"<$dir/contigs") || die "could not read-open $dir/contigs";
119 :    
120 : overbeek 1.36 my ($id, $seqP);
121 :     my $GC_count = 0;
122 :     my $AT_count = 0;
123 : overbeek 1.45 my $trouble = 0;
124 : overbeek 1.36 while (($id, $seqP) = &FIG::read_fasta_record(\*CONTIGS))
125 : overbeek 1.1 {
126 : overbeek 1.45 if (defined($contig_lens->{$id})) {
127 :     $trouble = 1;
128 :     warn "Duplicate contig ID: $id\n";
129 :     next;
130 :     }
131 :    
132 : overbeek 1.26 $contig_seqs->{$id} = lc($$seqP);
133 :     $contig_lens->{$id} = length($$seqP);
134 : overbeek 1.36
135 :     $GC_count += ($$seqP =~ tr/gcGC//);
136 :     $AT_count += ($$seqP =~ tr/atAT//);
137 : overbeek 1.1 }
138 :     close(CONTIGS);
139 : overbeek 1.45
140 :     if ($trouble) {
141 :     die "\nAborted due to duplicate contig IDs\n\n";
142 :     }
143 :    
144 : overbeek 1.26 $newG->{_contig_seqs} = $contig_seqs;
145 :     $newG->{_contig_lens} = $contig_lens;
146 : overbeek 1.36 $newG->{_GC_content} = ($GC_count + 1) / ($GC_count + $AT_count + 2);
147 : gdpusch 1.62
148 : gdpusch 1.61
149 : gdpusch 1.62 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
150 :     #...Auxilliary hashes to support insertion and deletion of features:
151 :     #-----------------------------------------------------------------------
152 :     $newG->{_used_list} = {};
153 :     # $newG->{_sort_left} = {};
154 :     # $newG->{_sort_right} = {};
155 :     # $newG->{_overlaps} = {};
156 :    
157 :    
158 :    
159 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
160 :     #...Quit here if no features will be imported:
161 :     #-----------------------------------------------------------------------
162 :     if (grep { m/^none/o } @features) {
163 :     return $newG;
164 :     }
165 :    
166 : overbeek 1.1
167 : gdpusch 1.62 ########################################################################
168 :     ########################################################################
169 :     ########################################################################
170 :    
171 :    
172 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
173 :     #...Import existing features:
174 :     #-----------------------------------------------------------------------
175 :    
176 :     if (-d "$dir/Features/rna") {
177 :     print STDERR "Directory $dir/Features/rna exists, and will be loaded\n" if $ENV{VERBOSE};
178 : overbeek 1.6 }
179 : gdpusch 1.64 # else {
180 :     # if (grep { m/^(rna|all)$/io } @features) {
181 :     # foreach my $entry ( $newG->find_rna() )
182 :     # {
183 :     # chomp $entry;
184 :     # my (undef, $seed_loc, $func) = split /\t/, $entry;
185 :     # $newG->add_feature( -type => 'rna'
186 :     # , -loc => $newG->from_seed_loc($seed_loc)
187 :     # , -func => $func )
188 :     # || die "Could not add RNA entry $entry";
189 :     # }
190 :     # }
191 :     # }
192 : overbeek 1.6
193 : overbeek 1.33 if (not @features) { $features[0] = 'all'; }
194 : overbeek 1.17 if ($features[0] eq 'all') {
195 : overbeek 1.19 print STDERR "Loading all features: " if $ENV{VERBOSE};
196 : overbeek 1.17 if (-d "$dir/Features") {
197 : overbeek 1.6 opendir(FEATURES, "$dir/Features") || confess "Could not opendir $dir/Features";
198 : overbeek 1.31 @features = grep { (-d "$dir/Features/$_") && !m/^\./o } readdir(FEATURES);
199 : overbeek 1.6 closedir(FEATURES) || confess "Could not closedir $dir/Features";
200 : gdpusch 1.49
201 : gdpusch 1.70 if ($ENV{VERBOSE}) {
202 :     if (@features) {
203 :     print STDERR join(", ", @features);
204 :     }
205 :     else {
206 :     print STDERR "No feature directories found.";
207 :     }
208 :    
209 :     print STDERR "\n";
210 : gdpusch 1.49 }
211 : overbeek 1.6 }
212 : overbeek 1.17 else {
213 : overbeek 1.33 confess "No directory '$dir/Features' found";
214 : overbeek 1.17 }
215 : overbeek 1.6 }
216 :    
217 :     if (@features)
218 :     {
219 :     if ($newG->import_features(@features)) {
220 :     print STDERR ("Imported ", join(', ', @features), "\n") if $ENV{VERBOSE};
221 :     }
222 :     }
223 : overbeek 1.4
224 :     return $newG;
225 :     }
226 :    
227 :     sub import_features
228 :     {
229 :     my ($self, @types) = @_;
230 : overbeek 1.26 my (%orf_lens, $orf_len);
231 :    
232 :     my $dir = $self->get_genome_dir();
233 : overbeek 1.7 my $taxid = $self->get_taxid();
234 : overbeek 1.4
235 :     if (not @types)
236 : overbeek 1.1 {
237 : overbeek 1.4 print STDERR "No feature-types specified --- defaulting to all subdirectories in $dir/Features\n"
238 :     if $ENV{VERBOSE};
239 :    
240 :     opendir(FEATURES, "$dir/Features")
241 :     || confess "Could not opendir $dir/Features";
242 :    
243 : overbeek 1.31 (@types = grep { (not m/^\./o) && (-d "$dir/Features/$_") } readdir(FEATURES))
244 : overbeek 1.4 || confess "Could not find any subdirectories in $dir/Features";
245 :    
246 :     closedir(FEATURES) || confess "Could not closedir $dir/Features";
247 :     }
248 :    
249 : gdpusch 1.70
250 : gdpusch 1.68 foreach my $type (@types) {
251 : overbeek 1.4 my $feature_dir = "$dir/Features/$type";
252 :     print STDERR "Importing features from $feature_dir\n" if $ENV{VERBOSE};
253 :    
254 : gdpusch 1.70 if (!-d "$feature_dir") {
255 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1)) {
256 :     cluck "WARNING: Feature directory $feature_dir does not exist";
257 :     }
258 :     else {
259 :     warn "WARNING: Feature directory $feature_dir does not exist";
260 :     }
261 :     next;
262 :     }
263 :    
264 :     if (!-s "$feature_dir/tbl") {
265 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1)) {
266 :     cluck "WARNING: Tbl file $feature_dir/tbl does not exist";
267 :     }
268 :     else {
269 :     warn "WARNING: Tbl file $feature_dir/tbl does not exist";
270 :     }
271 :     next;
272 :     }
273 :    
274 :    
275 : gdpusch 1.68 if (($type eq 'orf') && (-s "$feature_dir/orf_len")) {
276 : overbeek 1.26 open(ORF_LEN, "<$feature_dir/orf_len")
277 :     || confess "Could not read-open $feature_dir/orf_len";
278 :    
279 : overbeek 1.31 %orf_lens = map { m/^(\S+)\t(\d+)/o ? ($1 => $2)
280 :     : do { warn("Malformed ORF-length entry $_");
281 :     ();
282 :     }
283 :     } <ORF_LEN>;
284 : gdpusch 1.73 close(ORF_LEN) || die "Could not close $feature_dir/orf_len";
285 : overbeek 1.26 }
286 :    
287 : gdpusch 1.70
288 :     my %tmp_seq_of;
289 :     undef %tmp_seq_of;
290 :    
291 :     if (!-s "$feature_dir/fasta") {
292 :     confess "Fasta file $feature_dir/fasta does not exist";
293 : overbeek 1.4 }
294 : gdpusch 1.68 else {
295 : gdpusch 1.70 open(FASTA, "<$feature_dir/fasta") || confess "Could not read-open $feature_dir/fasta";
296 :     my $entry;
297 :     while (my ($fid, $seqP) = &FIG::read_fasta_record(\*FASTA)) {
298 :     if (length($$seqP)) {
299 :     print STDERR "Caching sequence for feature $fid\n" if $ENV{VERBOSE};
300 :    
301 :     if (($type eq 'peg') || ($type eq 'orf')) {
302 :     $$seqP =~ s/\*$//o;
303 :     }
304 :     $tmp_seq_of{$fid} = $$seqP;
305 :     }
306 :     else {
307 :     warn "Sequence for $fid has zero-length FASTA entry --- not caching\n";
308 :     }
309 :     }
310 : gdpusch 1.73 close(FASTA) || die "Could not close $feature_dir/fasta";
311 : gdpusch 1.70 }
312 : gdpusch 1.73
313 : gdpusch 1.70 open(TBL, "<$feature_dir/tbl") || confess "Could not read-open $feature_dir/tbl";
314 :    
315 :     my ($entry, $fid, $locus, $rest);
316 :     my ($contig_id, $beg, $end, $loc);
317 :     while (defined($entry = <TBL>)) {
318 :     chomp $entry;
319 :     if ($entry =~ m/^(\S+)\t(\S+)/o) {
320 :     ($fid, $locus) = ($1, $2);
321 :     } else {
322 :     confess "Malformed TBL entry: $entry";
323 :     }
324 : gdpusch 1.68
325 : gdpusch 1.70 if ($fid =~ m/^fig\|(\d+\.\d+)\.([^\/]+)\.(\d+)/o) {
326 :     confess "Taxon-ID $taxid does not match FID $fid" unless ($1 eq $taxid);
327 :     confess "Type $type does not match FID $fid" unless ($2 eq $type);
328 : gdpusch 1.68 }
329 :     else {
330 : gdpusch 1.70 confess "Malformed FID in $type entry $entry";
331 :     }
332 :    
333 :     $loc = $self->from_seed_loc($locus);
334 :    
335 :     my $func = "";
336 :     my $aliases = [];
337 :     if (($type eq 'rna') && ($entry =~ m/^\S+\t\S+\t(\S.*)$/o)) {
338 :     $func = $1;
339 :     }
340 :     elsif ($entry =~ m/^\S+\t\S+\t(\S.*\S)$/o) {
341 :     @$aliases = [ split /\t/, $1 ];
342 :     }
343 :    
344 :     if (($type eq 'peg') || ($type eq 'orf')) {
345 :     my $dna_seq = $self->get_dna_subseq($loc);
346 :     my $start_codon = substr($dna_seq, 0, 3);
347 :     if (($start_codon !~ m/[agt]tg/io) && (not $self->possibly_truncated($locus)))
348 :     {
349 :     if ($fid) {
350 : gdpusch 1.72 warn "WARNING: Feature $fid has invalid START codon '$start_codon'.\n";
351 : gdpusch 1.68 }
352 :     else {
353 : gdpusch 1.72 warn "WARNING: Feature at loc="
354 : gdpusch 1.70 , $self->flatten_dumper($loc),
355 : gdpusch 1.72 , " has invalid START codon '$start_codon'.\n";
356 : gdpusch 1.68 }
357 : gdpusch 1.72 # next;
358 : gdpusch 1.68 }
359 :     }
360 :    
361 : gdpusch 1.70
362 :     if (not defined($tmp_seq_of{$fid})) {
363 :     warn "Attempt to add feature $fid without corresponding FASTA entry --- skipping\n";
364 :     next;
365 :     }
366 :    
367 :     if ($type eq 'orf') {
368 :     if (!-s "$feature_dir/orf_len") {
369 :     $loc = $self->search_for_upstream_stop($loc);
370 : overbeek 1.17 }
371 :     else {
372 : gdpusch 1.70 if ($orf_len = $orf_lens{$fid}) {
373 :     $self->set_orf_length($loc, $orf_len);
374 : overbeek 1.31 }
375 :     else {
376 : gdpusch 1.70 confess "FID $fid has an undefined ORF length";
377 : overbeek 1.13 }
378 : overbeek 1.1 }
379 :     }
380 : gdpusch 1.70
381 :     $self->add_feature( -type => $type,
382 :     -fid => $fid,
383 :     -loc => $loc,
384 :     -seq => $tmp_seq_of{$fid},
385 :     -func => $func,
386 :     -aliases => $aliases
387 :     )
388 :     || warn "Could not add feature for $feature_dir/tbl entry $entry";
389 :    
390 :     delete $tmp_seq_of{$fid};
391 :     }
392 : gdpusch 1.73 close(TBL) || die "Could not close $feature_dir/tbl";
393 :    
394 : gdpusch 1.70
395 :     if (my @leftovers = keys %tmp_seq_of) {
396 :     die "Unprocessed FASTA entries left in sequence cache --- aborting: "
397 :     , join(qq(, ), @leftovers);
398 : overbeek 1.1 }
399 :     }
400 :    
401 : gdpusch 1.70
402 : overbeek 1.15 if (-s "$dir/assigned_functions")
403 :     {
404 : overbeek 1.26 print STDERR "Opening file $dir/assigned_functions\n" if $ENV{VERBOSE};
405 : overbeek 1.15 open(FUNC, "<$dir/assigned_functions")
406 :     || confess "Could not read-open $dir/assigned_functions";
407 :    
408 :     my $entry;
409 :     while (defined($entry = <FUNC>))
410 :     {
411 : overbeek 1.26 print STDERR "Reading entry: $entry" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
412 :    
413 : overbeek 1.15 chomp $entry;
414 : overbeek 1.31 if ($entry =~ m/^(\S+)\t(.*)$/o)
415 : overbeek 1.15 {
416 : gdpusch 1.55 my ($fid, $func) = ($1, $2);
417 :     if ($self->is_feature($fid)) {
418 :     $self->set_function($fid, $func)
419 :     || confess "Could not set function for FID $fid to $func";
420 :     }
421 :     else {
422 :     warn "Skipping attempt to set function of undefined feature $fid" if $ENV{VERBOSE};
423 :     }
424 : overbeek 1.15 }
425 :     else
426 :     {
427 :     confess "Malformed entry in $dir/assigned_functions:\n$entry";
428 :     }
429 :     }
430 : gdpusch 1.73 close(FUNC) || die "Could not close $dir/assigned_functions";
431 : overbeek 1.15 }
432 : overbeek 1.26 else {
433 :     print STDERR "File $dir/assigned_functions is empty\n" if $ENV{VERBOSE};
434 :     }
435 : overbeek 1.15
436 : overbeek 1.4 return 1;
437 :     }
438 :    
439 :     sub export_features
440 :     {
441 :     my ($self, @types) = @_;
442 :     my $dir = $self->get_genome_dir();
443 :    
444 : overbeek 1.13 if ((not @types) || ($types[0] =~ m/all/o)) { @types = $self->get_feature_types(); }
445 : overbeek 1.15
446 :     open(FUNC, ">$dir/assigned_functions")
447 :     || confess "Could not write-open $dir/assigned_functions";
448 :    
449 : overbeek 1.4 foreach my $type (@types)
450 :     {
451 :     my $feature_dir = "$dir/Features/$type";
452 : overbeek 1.15
453 : overbeek 1.17 my @fids = $self->get_fids_for_type($type);
454 : overbeek 1.20 if (@fids > 0) {
455 :     print STDERR "Exporting ", (scalar @fids), " $type features\n" if $ENV{VERBOSE};
456 : gdpusch 1.50
457 :     use File::Path;
458 :     (-d $feature_dir)
459 :     || mkpath($feature_dir, 0, 0777)
460 :     || confess "Could not create path $feature_dir";
461 :    
462 :     open(TBL, ">$feature_dir/tbl")
463 :     || confess "Could not write-open $feature_dir/tbl";
464 :    
465 :     open(FASTA, ">$feature_dir/fasta")
466 :     || confess "Could not write-open $feature_dir/fasta";
467 :    
468 : gdpusch 1.69 if ($type eq 'orf') {
469 : gdpusch 1.50 open(ORF_LEN, ">$feature_dir/orf_len")
470 :     || confess "Could not write-open $feature_dir/orf_len";
471 : gdpusch 1.69
472 :     foreach my $fid (@fids) {
473 :     my $orf_loc = $self->get_feature_loc($fid);
474 :     my $orf_len = $self->get_orf_length($orf_loc);
475 :     if ((not defined($orf_len)) || (not $orf_len)) {
476 :     confess "Could not fetch ORF boundaries for $fid:\n"
477 :     , Dumper($self->{_features}->{$fid});
478 :     }
479 :     print ORF_LEN "$fid\t$orf_len\n";
480 :     }
481 :     close(ORF_LEN) || die "Could not close $feature_dir/orf_len";
482 : gdpusch 1.50 }
483 :    
484 :     foreach my $fid (@fids)
485 :     {
486 :     my $seed_loc = $self->get_seed_loc($fid);
487 :     $seed_loc || confess "Could not get SEED-format locus for $fid";
488 :     print STDERR "Exporting FID $fid,"
489 :     , " loc=", &flatten_dumper($self->get_feature_loc($fid))
490 :     , " seed_loc=$seed_loc\n"
491 :     if $ENV{VERBOSE};
492 :    
493 :     my $seq = $self->get_feature_sequence($fid);
494 :     $seq || confess "Could not get sequence for $fid";
495 :    
496 :     print TBL "$fid\t$seed_loc\t\n";
497 :     &FIG::display_id_and_seq($fid, \$seq, \*FASTA);
498 :    
499 :     my $func = $self->get_function($fid);
500 :     print FUNC "$fid\t$func\n" if $func;
501 :     }
502 :     close(FASTA) || die "Could not close $feature_dir/fasta";
503 :     close(TBL) || die "Could not close feature_dir/tbl";
504 : overbeek 1.20 }
505 :     else {
506 : overbeek 1.25 if (-d $feature_dir) {
507 : gdpusch 1.50 print STDERR "No exportable features of type $type --- removing directory $feature_dir\n"
508 :     if $ENV{VERBOSE};
509 :     system("/bin/rm -fRv $feature_dir")
510 : gdpusch 1.49 && warn "Could not remove directory $feature_dir";
511 : overbeek 1.25 }
512 : overbeek 1.20 }
513 : overbeek 1.4 }
514 : gdpusch 1.73 close(FUNC) || die "Could not close $dir/assigned_functions";
515 : overbeek 1.4
516 :     return 1;
517 : overbeek 1.1 }
518 :    
519 :    
520 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
521 :     # ... Accessor functions ...
522 :     #-----------------------------------------------------------------------
523 :    
524 :     sub get_taxid
525 :     {
526 :     my ($self) = @_;
527 :     return $self->{_taxid};
528 :     }
529 :    
530 :     sub get_genome_dir
531 :     {
532 :     my ($self) = @_;
533 :     return $self->{_dir};
534 :     }
535 :    
536 :     sub get_genome_name
537 :     {
538 :     my ($self) = @_;
539 :     return $self->{_genome};
540 :     }
541 :    
542 :     sub get_project
543 :     {
544 :     my ($self) = @_;
545 :     return $self->{_project};
546 :     }
547 :    
548 :     sub get_taxonomy
549 :     {
550 :     my ($self) = @_;
551 :     return $self->{_taxonomy};
552 :     }
553 :    
554 : gdpusch 1.49 sub get_genetic_code_number {
555 :     my ($self) = @_;
556 :     return $self->{_genetic_code_number};
557 :     }
558 :    
559 :     sub get_translation_table {
560 :     my ($self) = @_;
561 :     return $self->{_translation_table};
562 :     }
563 :    
564 : overbeek 1.36 sub get_GC_content {
565 :     my ($self) = @_;
566 :     return $self->{_GC_content};
567 :     }
568 :    
569 : overbeek 1.1 sub get_contig_ids
570 :     {
571 :     my ($self) = @_;
572 : overbeek 1.26 return ( sort keys % { $self->{_contig_lens} } );
573 : overbeek 1.1 }
574 :    
575 : overbeek 1.7 sub get_contig_length
576 : overbeek 1.1 {
577 :     my ($self, $contig_id) = @_;
578 :     confess "No contig_id" unless $contig_id;
579 :    
580 :     my $len;
581 : overbeek 1.26 if ($len = $self->{_contig_lens}->{$contig_id})
582 : overbeek 1.1 {
583 :     return $len;
584 :     }
585 :     else
586 :     {
587 :     confess "No sequence for contig_id $contig_id";
588 : overbeek 1.26 return undef;
589 : overbeek 1.1 }
590 :     }
591 :    
592 : overbeek 1.26 sub get_contig_seqP
593 : overbeek 1.1 {
594 :     my ($self, $contig_id) = @_;
595 :     confess "No contig_id" unless $contig_id;
596 :    
597 : overbeek 1.26 if (defined($self->{_contig_seqs}->{$contig_id}))
598 : overbeek 1.1 {
599 : overbeek 1.26 return \do{ $self->{_contig_seqs}->{$contig_id} };
600 : overbeek 1.1 }
601 :     else
602 :     {
603 :     confess "No sequence for contig_id $contig_id";
604 : overbeek 1.26 return undef;
605 : overbeek 1.1 }
606 :     }
607 :    
608 : gdpusch 1.74
609 :    
610 : overbeek 1.1 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
611 :     # ... Feature accessors ...
612 :     #-----------------------------------------------------------------------
613 :    
614 : overbeek 1.2 sub is_feature
615 :     {
616 :     my ($self,$fid) = @_;
617 :    
618 :     return defined($self->{_features}->{$fid});
619 :     }
620 :    
621 : overbeek 1.18 sub is_used {
622 :     my ($self, $loc) = @_;
623 :    
624 :     my $contig = $loc->[_CONTIG];
625 :     my $strand = $loc->[_LOC_END];
626 :     my $end = $loc->[_STRAND];
627 :     my $key = $contig.$strand.$end;
628 :    
629 :     return $self->{_used_list}->{$key};
630 :     }
631 :    
632 : overbeek 1.4 sub get_feature_types
633 :     {
634 :     my ($self) = @_;
635 :    
636 :     return keys % { $self->{_features}->{_maxnum} };
637 :     }
638 :    
639 : overbeek 1.2 sub get_fids_for_type
640 :     {
641 : overbeek 1.25 my ($self, @types) = @_;
642 :     my $type = join('|', @types);
643 :    
644 : overbeek 1.13 my $patt;
645 :     if ((not $type) || ($type =~ m/all/o)) {
646 : overbeek 1.31 $patt = qr/^[^_]\w+$/o;
647 : overbeek 1.13 } else {
648 : overbeek 1.15 $patt = qr/^$type$/;
649 : overbeek 1.13 }
650 : overbeek 1.2
651 : overbeek 1.25 my $featureP = $self->{_features};
652 : overbeek 1.17 if (wantarray) {
653 :     return sort { &FIG::by_fig_id($a, $b) }
654 :     grep { $featureP->{$_}->{_type} =~ m/$patt/ }
655 :     keys %$featureP;
656 :     }
657 :     else {
658 :     return (scalar grep { $featureP->{$_}->{_type} =~ m/$patt/ } keys %$featureP);
659 :     }
660 : overbeek 1.2 }
661 :    
662 : overbeek 1.1 sub get_feature_object
663 :     {
664 :     my ($self, $fid) = @_;
665 :    
666 : overbeek 1.7 return undef if not $self->is_feature($fid);
667 : overbeek 1.1 return $self->{_features}->{$fid};
668 :     }
669 :    
670 :     sub get_feature_type
671 :     {
672 :     my ($self, $fid) = @_;
673 :    
674 : overbeek 1.7 return undef if not $self->is_feature($fid);
675 : overbeek 1.1 return $self->{_features}->{$fid}->{_type};
676 :     }
677 :    
678 :     sub get_feature_loc
679 :     {
680 : overbeek 1.4 my ($self, $fid) = @_;
681 : overbeek 1.1
682 : overbeek 1.7 return undef if not $self->is_feature($fid);
683 : overbeek 1.1 return $self->{_features}->{$fid}->{_loc};
684 :     }
685 :    
686 : overbeek 1.7 sub get_feature_length
687 :     {
688 :     my ($self, $fid) = @_;
689 :    
690 :     return undef if not $self->is_feature($fid);
691 :    
692 :     my $beg = $self->get_feature_beginpoint($fid);
693 :     my $end = $self->get_feature_endpoint($fid);
694 :    
695 :     return (1 + abs($end-$beg));
696 :     }
697 :    
698 : overbeek 1.19
699 :     sub set_feature_aliases
700 :     {
701 :     my ($self, $fid, @aliases) = @_;
702 :    
703 :     return undef if not $self->is_feature($fid);
704 :    
705 :     return $self->{_features}->{$fid}->{_aliases} = [ @aliases ];
706 :     }
707 :    
708 :     sub get_feature_aliases
709 :     {
710 :     my ($self, $fid) = @_;
711 :    
712 :     return undef if not $self->is_feature($fid);
713 :    
714 :     return $self->{_features}->{$fid}->{_aliases};
715 :     }
716 :    
717 :    
718 : gdpusch 1.63 sub from_seed_loc {
719 : overbeek 1.6 my ($self, $loc) = @_;
720 :    
721 :     my @loc = ();
722 :     my @seed_loc = split /,/, $loc;
723 : gdpusch 1.63 foreach my $exon (@seed_loc) {
724 :     if ($exon =~ m/(\S+)_(\d+)_(\d+)/o) {
725 : overbeek 1.6 my ($contig_id, $beg, $end) = ($1, $2, $3);
726 :     my $len = 1 + abs($end - $beg);
727 : overbeek 1.26 my $strand = ($end > $beg) ? qq(+) : qq(-) ;
728 : overbeek 1.6
729 :     push @loc, [$contig_id, $end, $strand, $len];
730 :     }
731 : gdpusch 1.63 else {
732 : overbeek 1.6 confess "Could not parse exon $exon";
733 :     }
734 :     }
735 :    
736 :     return [@loc];
737 :     }
738 :    
739 : gdpusch 1.63 sub to_seed_loc {
740 : gdpusch 1.59 my ($self, $loc) = @_;
741 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
742 :    
743 :     my @exon_locs = ();
744 : gdpusch 1.63 foreach my $exon (@$loc) {
745 : gdpusch 1.59 my ($contig_id, $end, $strand, $len) = @$exon;
746 :     my $sign = ($strand eq qq(+)) ? +1 : -1 ;
747 :     my $beg = $end - $sign*($len-1);
748 :     push @exon_locs, join(qq(_), ($contig_id, $beg, $end));
749 :     }
750 :    
751 :     return join(qq(,), @exon_locs);
752 :     }
753 :    
754 :    
755 : gdpusch 1.63 sub get_seed_loc {
756 : overbeek 1.4 my ($self, $fid) = @_;
757 : overbeek 1.7 return undef if not $self->is_feature($fid);
758 : overbeek 1.1
759 :     my $seed_loc = "";
760 : gdpusch 1.69 my @loc = @ { $self->get_feature_loc($fid) };
761 :     foreach my $exon (@loc) {
762 : overbeek 1.1 my ($contig_id, $end, $strand, $len) = @$exon;
763 : overbeek 1.26 my $sign = ($strand eq qq(+)) ? +1 : -1 ;
764 : overbeek 1.1 my $beg = $end - $sign*($len-1);
765 : gdpusch 1.67 $exon = join('_', ($contig_id, $beg, $end));
766 : overbeek 1.1 }
767 : gdpusch 1.69 $seed_loc = join(",", @loc);
768 : overbeek 1.1
769 :     return $seed_loc;
770 :     }
771 :    
772 : gdpusch 1.63 sub set_feature_loc {
773 : overbeek 1.1 my ($self, $fid, $loc) = @_;
774 : gdpusch 1.63 if (not $self->is_feature($fid)) {
775 : overbeek 1.1 confess "Attempt to set location of undefined feature $fid";
776 :     }
777 :    
778 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
779 :    
780 : overbeek 1.12 return ($self->{_features}->{$fid}->{_loc} = $loc);
781 : overbeek 1.1 }
782 :    
783 : gdpusch 1.63 sub get_feature_contig {
784 : overbeek 1.1 my ($self, $fid) = @_;
785 : overbeek 1.7 return undef if not $self->is_feature($fid);
786 : overbeek 1.1
787 : overbeek 1.7 my $loc = $self->get_feature_loc($fid);
788 :     return $loc->[0]->[_CONTIG];
789 : overbeek 1.1 }
790 :    
791 : overbeek 1.7
792 : gdpusch 1.63 sub get_feature_strand {
793 : overbeek 1.7 my ($self, $fid) = @_;
794 :     return undef if not $self->is_feature($fid);
795 : overbeek 1.1
796 : overbeek 1.7 my $loc = $self->get_feature_loc($fid);
797 :     return $loc->[0]->[_STRAND];
798 : overbeek 1.1 }
799 :    
800 : gdpusch 1.63 sub get_feature_sign {
801 : overbeek 1.9 my ($self, $fid) = @_;
802 :     return undef if not $self->is_feature($fid);
803 :    
804 :     my $loc = $self->get_feature_loc($fid);
805 : overbeek 1.26 return (($loc->[0]->[_STRAND] eq qq(+)) ? +1 : -1);
806 : overbeek 1.9 }
807 :    
808 : overbeek 1.1
809 : gdpusch 1.63 sub get_feature_endpoint {
810 : overbeek 1.1 my ($self, $fid) = @_;
811 : overbeek 1.12 return undef if not $self->is_feature($fid);
812 : overbeek 1.26 confess "Problem with feature $fid:\n", &flatten_dumper($self->{_features}->{$fid})
813 : overbeek 1.12 unless ( $self->{_features}->{$fid}->{_loc}
814 :     && $self->{_features}->{$fid}->{_loc}->[-1]
815 :     && $self->{_features}->{$fid}->{_loc}->[-1]->[_LOC_END]);
816 : overbeek 1.1
817 : overbeek 1.7 return $self->{_features}->{$fid}->{_loc}->[-1]->[_LOC_END];
818 : overbeek 1.1 }
819 :    
820 : gdpusch 1.63 sub set_feature_endpoint {
821 : overbeek 1.1 my ($self, $fid, $end) = @_;
822 : overbeek 1.7 if (not $self->is_feature($fid)) {
823 : overbeek 1.1 $self->{_features}->{$fid} = [];
824 :     }
825 :    
826 : overbeek 1.7 $self->{_features}->{$fid}->{_loc}->[-1]->[_LOC_END] = $end;
827 : overbeek 1.1 }
828 :    
829 :    
830 : gdpusch 1.63 sub get_feature_beginpoint {
831 : overbeek 1.1 my ($self, $fid) = @_;
832 :    
833 : overbeek 1.7 return undef if not $self->is_feature($fid);
834 :    
835 : overbeek 1.12 my $loc = $self->get_feature_loc($fid);
836 : overbeek 1.26 my $sign = ($self->get_feature_strand($fid) eq qq(+)) ? +1 : -1 ;
837 : overbeek 1.7 my $contig_length = $self->get_contig_length($self->get_feature_contig($fid));
838 :    
839 : overbeek 1.26 my $end = $self->get_exon_end($loc->[-1]);
840 :     my $len = $self->get_exon_length($loc->[-1]);
841 : overbeek 1.7 my $beg = &FIG::max( 1, &FIG::min( ($end - $sign*($len-1) ), $contig_length ) );
842 :    
843 :     # print STDERR "BEG: fid=$fid\tsign=$sign\tend=$end\tlen=$len\tbeg=$beg\tcontig_length=$contig_length\n";
844 :     return $beg;
845 : overbeek 1.1 }
846 :    
847 : overbeek 1.7
848 : gdpusch 1.63 sub get_feature_leftbound {
849 : overbeek 1.7 my ($self, $fid) = @_;
850 :    
851 :     return undef if not $self->is_feature($fid);
852 : overbeek 1.1
853 : overbeek 1.7 my $beg = $self->get_feature_beginpoint($fid);
854 :     my $end = $self->get_feature_endpoint($fid);
855 :    
856 :     return &FIG::min($beg, $end);
857 :     }
858 :    
859 : gdpusch 1.63 sub get_feature_rightbound {
860 : overbeek 1.7 my ($self, $fid) = @_;
861 :    
862 :     return undef if not $self->is_feature($fid);
863 :    
864 :     my $beg = $self->get_feature_beginpoint($fid);
865 :     my $end = $self->get_feature_endpoint($fid);
866 : overbeek 1.1
867 : overbeek 1.7 return &FIG::max($beg, $end);
868 : overbeek 1.1 }
869 :    
870 :    
871 : gdpusch 1.63 sub get_feature_dna {
872 : overbeek 1.1 my ($self, $fid) = @_;
873 :    
874 : overbeek 1.7 return undef if not $self->is_feature($fid);
875 : overbeek 1.1 return $self->get_dna_subseq($self->{_features}->{$fid}->{_loc});
876 :     }
877 :    
878 :    
879 : gdpusch 1.63 sub get_feature_sequence {
880 : overbeek 1.1 my ($self, $fid) = @_;
881 :    
882 : overbeek 1.7 return undef if not $self->is_feature($fid);
883 : overbeek 1.1 return undef if not $self->{_features}->{$fid}->{_seq};
884 :    
885 : gdpusch 1.75 my $seq = $self->{_features}->{$fid}->{_seq};
886 :    
887 :     if ($self->get_feature_type($fid) eq qq(orf)) {
888 :     my $len = $self->get_feature_length($fid);
889 :     # warn (qq(ORF $fid:\t), $len, qq(\t), ($len - 3)/3, qq(\n));
890 :     $seq = substr($seq, -($len - 3)/3);
891 :     $seq =~ s/^[MLV]/M/;
892 :     }
893 :    
894 :     return $seq;
895 : overbeek 1.1 }
896 :    
897 : gdpusch 1.63 sub set_feature_sequence {
898 : overbeek 1.1 my ($self, $fid, $seq) = @_;
899 : overbeek 1.4 my $type = $self->get_feature_type($fid);
900 : overbeek 1.1
901 : gdpusch 1.63 if (not defined($type)) {
902 : overbeek 1.1 confess "Attempt to set sequence of undefined feature $fid";
903 :     }
904 :    
905 : gdpusch 1.63 if (($type eq 'peg') || ($type eq 'orf')) {
906 : gdpusch 1.65 if ($seq =~ m/[^ABCDEFGHIKLMNPQRSTVWXxYZ]/o) {
907 : overbeek 1.1 $seq =~ tr/ABCDEFGHIKLMNPQRSTVWXYZ/abcdefghiklmnpqrstvwxyz/;
908 :     confess "Translation for FID $fid contains invalid chars:\n"
909 :     , $seq;
910 :     }
911 :     }
912 :    
913 :     $self->{_features}->{$fid}->{_seq} = $seq;
914 :     }
915 :    
916 : gdpusch 1.63 sub get_function {
917 : overbeek 1.15 my ($self, $fid) = @_;
918 :    
919 : overbeek 1.20 if (not $self->is_feature($fid)) {
920 :     print STDERR "WARNING: Attempt to get function of undefined feature $fid\n";
921 :     return undef;
922 :     }
923 : overbeek 1.15
924 :     return $self->{_features}->{$fid}->{_func};
925 :     }
926 :    
927 : gdpusch 1.63 sub set_function {
928 : overbeek 1.15 my ($self, $fid, $func) = @_;
929 : overbeek 1.26 print STDERR "Setting function of $fid to $func\n"
930 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
931 : overbeek 1.15
932 : overbeek 1.20 if (not $self->is_feature($fid)) {
933 :     print STDERR "WARNING: Attempt to set function of undefined feature $fid to $func\n";
934 :     return undef;
935 :     }
936 : overbeek 1.15
937 :     return ($self->{_features}->{$fid}->{_func} = $func);
938 :     }
939 :    
940 : gdpusch 1.63 sub get_exon_contig {
941 : overbeek 1.18 my ($self, $loc) = @_;
942 :    
943 :     return $loc->[_CONTIG];
944 :     }
945 :    
946 : gdpusch 1.63 sub get_exon_end {
947 : overbeek 1.12 my ($self, $loc) = @_;
948 : overbeek 1.7
949 : overbeek 1.12 return $loc->[_LOC_END];
950 : overbeek 1.7 }
951 :    
952 : gdpusch 1.63 sub get_exon_length {
953 : overbeek 1.12 my ($self, $loc) = @_;
954 : overbeek 1.1
955 : overbeek 1.12 return $loc->[_LENGTH];
956 : overbeek 1.1 }
957 :    
958 : gdpusch 1.63 sub set_exon_length {
959 : overbeek 1.12 my ($self, $loc, $len) = @_;
960 : overbeek 1.26 my $orf_length;
961 :    
962 :     if (defined($orf_length = $self->get_orf_length($loc))) {
963 :     if ($len <= $orf_length) {
964 : overbeek 1.29 return ($loc->[_LENGTH] = $len);
965 : overbeek 1.26 }
966 :     else {
967 :     return undef;
968 :     }
969 :     }
970 : overbeek 1.1
971 : overbeek 1.29 return ($loc->[_LENGTH] = $len);
972 : overbeek 1.1 }
973 :    
974 :    
975 : gdpusch 1.63 sub get_orf_length {
976 : overbeek 1.12 my ($self, $loc) = @_;
977 :     my $orf_len;
978 :    
979 : gdpusch 1.63 if (ref($loc) eq 'ARRAY') {
980 : overbeek 1.12 if (ref($loc->[-1]) ne 'ARRAY') { $loc = [$loc]; }
981 :    
982 : gdpusch 1.63 if (defined($orf_len = $loc->[-1]->[_ORF_LEN])) {
983 : overbeek 1.12 return $orf_len;
984 :     }
985 :     }
986 : gdpusch 1.63 else {
987 : overbeek 1.26 confess "Not a location object:\n", &flatten_dumper($loc);
988 : overbeek 1.1 }
989 :    
990 : overbeek 1.26 return undef;
991 : overbeek 1.1 }
992 :    
993 : gdpusch 1.63 sub set_orf_length {
994 : overbeek 1.12 my ($self, $loc, $orf_len) = @_;
995 :    
996 : gdpusch 1.63 if (ref($loc) eq 'ARRAY') {
997 : overbeek 1.12 if (ref($loc->[-1]) eq 'ARRAY') {
998 : overbeek 1.19 print STDERR "Setting ORF-len RefRef\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
999 : overbeek 1.12 $loc->[-1]->[_ORF_LEN] = $orf_len;
1000 :     } else {
1001 : overbeek 1.19 print STDERR "Setting ORF-len Ref\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
1002 : overbeek 1.12 $loc->[_ORF_LEN] = $orf_len;
1003 :     }
1004 :     }
1005 : gdpusch 1.63 else {
1006 : overbeek 1.12 confess STDERR "Attempt to set ORF length to $orf_len for a non-location object:\n"
1007 : overbeek 1.26 , &flatten_dumper($loc);
1008 : overbeek 1.1 }
1009 :    
1010 : overbeek 1.15 return $loc;
1011 : overbeek 1.1 }
1012 :    
1013 :    
1014 : gdpusch 1.63 sub get_orf_sequence {
1015 : overbeek 1.1 my ($self, $fid) = @_;
1016 :    
1017 :     return undef if not $self->{_features}->{$fid};
1018 :     return $self->{_features}->{$fid}->[_ORF_TRANS];
1019 :     }
1020 :    
1021 : gdpusch 1.63 sub set_orf_sequence {
1022 : overbeek 1.1 my ($self, $fid, $trans) = @_;
1023 :     if (not $self->{_features}->{$fid}) {
1024 :     confess STDERR "Undefined feature $fid";
1025 : overbeek 1.26 return undef;
1026 : overbeek 1.1 }
1027 :    
1028 : gdpusch 1.63 if ($trans =~ m/[^ABCDEFGHIKLMNPQRSTVWXYZ]/o) {
1029 : overbeek 1.1 $trans =~ tr/ABCDEFGHIKLMNPQRSTVWXYZ/abcdefghiklmnpqrstvwxyz/;
1030 :     croak "Translation for FID $fid contains invalid chars:\n"
1031 :     , $trans;
1032 :     }
1033 :    
1034 :     $self->{_features}->{$fid}->[_ORF_TRANS] = $trans;
1035 :     }
1036 :    
1037 :    
1038 : gdpusch 1.63 sub get_orf_sims {
1039 : overbeek 1.1 my ($self, $fid) = @_;
1040 :    
1041 :     return undef if not $self->{_features}->{$fid};
1042 :     return $self->{_features}->{$fid}->[_ORF_SIMS];
1043 :     }
1044 :    
1045 : gdpusch 1.63 sub set_orf_sims {
1046 : overbeek 1.1 my ($self, $fid, $simobj) = @_;
1047 :     if (not $self->{_features}->{$fid}) {
1048 :     $self->{_features}->{$fid} = [];
1049 :     }
1050 :    
1051 :     $self->{_features}->{$fid}->[_ORF_SIMS] = $simobj;
1052 :     }
1053 :    
1054 : gdpusch 1.63 sub get_all_orfs {
1055 : overbeek 1.8 my ($self, $min_len) = @_;
1056 :    
1057 :     my $orfs = {};
1058 :     %$orfs = map { $_ => NewGenome::ORF->new($self, $_) } $self->get_fids_for_type('orf');
1059 :    
1060 :     return $orfs;
1061 :     }
1062 :    
1063 :    
1064 : overbeek 1.12 # sub get_overlapping_fids
1065 :     # {
1066 :     # my ($self, $fid) = @_;
1067 :     #
1068 :     # if (not $self->is_feature($fid))
1069 :     # {
1070 :     # confess "Attempt to find overlaps for non-existent feature $fid";
1071 :     # }
1072 :     #
1073 :     # my $overlaps = $self->{_overlaps};
1074 :     # return keys % { $overlaps->{$fid} };
1075 :     # }
1076 :    
1077 :     # sub get_overlap
1078 :     # {
1079 :     # my ($self, $fid1, $fid2) = @_;
1080 :     #
1081 :     # if (not $self->is_feature($fid1))
1082 :     # {
1083 :     # confess "Attempt to find overlaps for non-existent feature $fid1";
1084 :     # }
1085 :     #
1086 :     # if (not $self->is_feature($fid2))
1087 :     # {
1088 :     # confess "Attempt to find overlaps for non-existent feature $fid2";
1089 :     # }
1090 :     #
1091 :     # my $overlaps = $self->{_overlaps};
1092 :     # if (defined($overlaps->{$fid1}))
1093 :     # {
1094 :     # if (defined($overlaps->{$fid1}->{$fid2}))
1095 :     # {
1096 :     # return $overlaps->{$fid1}->{$fid2};
1097 :     # }
1098 :     # }
1099 :     #
1100 :     # return 0;
1101 :     # }
1102 : overbeek 1.10
1103 :    
1104 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1105 :     # ... Region Comparisons ...
1106 :     #-----------------------------------------------------------------------
1107 :    
1108 : gdpusch 1.63 sub contains_embedded {
1109 : overbeek 1.10 my ($self, $fid1, $fid2) = @_;
1110 :    
1111 :     my $beg1 = $self->get_feature_beginpoint( $fid1 );
1112 :     my $end1 = $self->get_feature_endpoint( $fid1 );
1113 :    
1114 :     my $beg2 = $self->get_feature_beginpoint( $fid2 );
1115 :     my $end2 = $self->get_feature_endpoint( $fid2 );
1116 :    
1117 :     if (&FIG::between($beg1, $beg2, $end1) && &FIG::between($beg1, $end2, $end1))
1118 :     {
1119 :     return $self->get_feature_length($fid2);
1120 :     }
1121 :    
1122 :     return 0;
1123 :     }
1124 :    
1125 :    
1126 :    
1127 : overbeek 1.1 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1128 :     # ... Utility functions ...
1129 :     #-----------------------------------------------------------------------
1130 :    
1131 : gdpusch 1.63 sub is_rna {
1132 : overbeek 1.12 my ($self, $fid) = @_;
1133 :    
1134 :     return ($self->is_feature($fid) && ($self->{_features}->{$fid}->{_type} eq 'rna'));
1135 :     }
1136 :    
1137 : gdpusch 1.55 sub translate {
1138 :     my ($self, $seq) = @_;
1139 :    
1140 :     return &FIG::translate($seq, $self->get_translation_table);
1141 :     }
1142 :    
1143 : gdpusch 1.63 sub _parse_exon {
1144 : overbeek 1.1 my ($exon) = @_;
1145 :    
1146 : gdpusch 1.63 if ($exon =~ m/^(\S+):(\d+)([+-])(\d+)$/o) {
1147 : overbeek 1.1 return ($1, $2, $3, $4);
1148 :     }
1149 : gdpusch 1.63 else {
1150 : overbeek 1.1 confess "Invalid exon string $exon";
1151 :     }
1152 :     }
1153 :    
1154 : gdpusch 1.63 sub make_exon {
1155 : overbeek 1.12 my ($self, @args) = @_;
1156 :    
1157 :     return [ @args ];
1158 : overbeek 1.1 }
1159 :    
1160 : gdpusch 1.63 sub make_loc {
1161 : overbeek 1.12 my ($self, @args) = @_;
1162 :    
1163 :     return [ [@args] ];
1164 :     }
1165 :    
1166 : gdpusch 1.63 sub copy_loc {
1167 : overbeek 1.12 my ($self, $loc) = @_;
1168 :     my $new_loc = [];
1169 :    
1170 : gdpusch 1.63 if (ref($loc->[0]) eq 'ARRAY') {
1171 : overbeek 1.12 foreach my $exon (@$loc) { push @$new_loc, [ @$exon ]; }
1172 :     }
1173 : gdpusch 1.63 else {
1174 : overbeek 1.12 @$new_loc = @$loc;
1175 :     }
1176 :    
1177 :     return $new_loc;
1178 : overbeek 1.1 }
1179 :    
1180 : gdpusch 1.63 sub parse_loc {
1181 : overbeek 1.1 my ($self, $loc_string) = @_;
1182 :    
1183 :     my $loc = [];
1184 :     my @exons = split /,/, $loc_string;
1185 : gdpusch 1.63 foreach my $exon (@exons) {
1186 : overbeek 1.1 push @$loc, [ &_parse_exon($exon) ];
1187 :     }
1188 :    
1189 :     return $loc;
1190 :     }
1191 :    
1192 : gdpusch 1.63 sub compare_loc {
1193 : overbeek 1.7 my ($self, $a, $b) = @_;
1194 :    
1195 : overbeek 1.12 return $self->compare_left($a, $b);
1196 :     }
1197 :    
1198 : gdpusch 1.63 sub compare_left {
1199 : overbeek 1.12 my ($self, $a, $b) = @_;
1200 :    
1201 : overbeek 1.7 my $A_contig = $self->get_feature_contig($a);
1202 :     my $B_contig = $self->get_feature_contig($b);
1203 :    
1204 :     my $A_strand = $self->get_feature_strand($a);
1205 :     my $B_strand = $self->get_feature_strand($b);
1206 :    
1207 :     my $A_left = $self->get_feature_leftbound($a);
1208 :     my $B_left = $self->get_feature_leftbound($b);
1209 :    
1210 :     my $A_right = $self->get_feature_rightbound($a);
1211 :     my $B_right = $self->get_feature_rightbound($b);
1212 :    
1213 :     my $A_len = $self->get_feature_length($a);
1214 :     my $B_len = $self->get_feature_length($b);
1215 :    
1216 : gdpusch 1.63 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2)) {
1217 : overbeek 1.7 print STDERR "A = ($a, $A_contig, $A_strand, $A_left, $A_right, $A_len)\n";
1218 :     print STDERR "B = ($b, $B_contig, $B_strand, $B_left, $B_right, $B_len)\n";
1219 :     }
1220 :    
1221 :     return ( ($A_contig cmp $B_contig)
1222 :     || ($A_left <=> $B_left)
1223 :     || ($B_len <=> $A_len)
1224 :     || ($A_strand cmp $B_strand)
1225 :     );
1226 :     }
1227 :    
1228 : overbeek 1.12
1229 : gdpusch 1.63 sub compare_right {
1230 : overbeek 1.12 my ($self, $a, $b) = @_;
1231 :    
1232 :     my $A_contig = $self->get_feature_contig($a);
1233 :     my $B_contig = $self->get_feature_contig($b);
1234 :    
1235 :     my $A_strand = $self->get_feature_strand($a);
1236 :     my $B_strand = $self->get_feature_strand($b);
1237 :    
1238 :     my $A_left = $self->get_feature_leftbound($a);
1239 :     my $B_left = $self->get_feature_leftbound($b);
1240 :    
1241 :     my $A_right = $self->get_feature_rightbound($a);
1242 :     my $B_right = $self->get_feature_rightbound($b);
1243 :    
1244 :     my $A_len = $self->get_feature_length($a);
1245 :     my $B_len = $self->get_feature_length($b);
1246 :    
1247 : gdpusch 1.63 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2)) {
1248 : overbeek 1.12 print STDERR "A = ($a, $A_contig, $A_strand, $A_left, $A_right, $A_len)\n";
1249 :     print STDERR "B = ($b, $B_contig, $B_strand, $B_left, $B_right, $B_len)\n";
1250 :     }
1251 :    
1252 :     return ( ($A_contig cmp $B_contig)
1253 :     || ($A_right <=> $B_right)
1254 :     || ($B_len <=> $A_len)
1255 :     || ($A_strand cmp $B_strand)
1256 :     );
1257 :     }
1258 :    
1259 :    
1260 : gdpusch 1.70 sub possibly_truncated {
1261 :     my($self, $loc) = @_;
1262 :    
1263 :     my $seed_loc;
1264 :     if (ref($loc) eq qq(ARRAY)) {
1265 :     $seed_loc = $self->to_seed_loc($loc);
1266 :     }
1267 :     else {
1268 :     $seed_loc = $loc;
1269 :     }
1270 :    
1271 : gdpusch 1.71 my ($contig_id, $beg, $end) = $fig->boundaries_of( $seed_loc );
1272 : gdpusch 1.70 my $contig_len = $self->get_contig_length( $contig_id );
1273 :    
1274 :     if ((! $self->_near_end($contig_len, $beg)) && (! $self->_near_end($contig_len, $end)))
1275 :     {
1276 :     return 0;
1277 :     }
1278 :     return 1;
1279 :     }
1280 :    
1281 :     sub _near_end {
1282 :     my($self, $contig_len, $x) = @_;
1283 :    
1284 :     return (($x < 300) || ($x > ($contig_len - 300)));
1285 :     }
1286 :    
1287 :    
1288 : gdpusch 1.63 sub overlap_between {
1289 : overbeek 1.7 my ($self, $a, $b) = @_;
1290 :    
1291 :     my $A_contig = $self->get_feature_contig($a);
1292 :     my $B_contig = $self->get_feature_contig($b);
1293 :     return 0 if ($A_contig ne $B_contig);
1294 :    
1295 :     my $A_left = $self->get_feature_leftbound($a);
1296 :     my $B_left = $self->get_feature_leftbound($b);
1297 :    
1298 :     my $A_right = $self->get_feature_rightbound($a);
1299 :     my $B_right = $self->get_feature_rightbound($b);
1300 :    
1301 : gdpusch 1.63 if ($A_left > $B_left) {
1302 : overbeek 1.7 ($A_left, $B_left) = ($B_left, $A_left);
1303 :     ($A_right, $B_right) = ($B_right, $A_right);
1304 :     }
1305 :    
1306 :     my $overlap = 0;
1307 :     if ($A_right >= $B_left) { $overlap = &FIG::min($A_right, $B_right) - $B_left + 1; }
1308 :    
1309 :     return $overlap;
1310 :     }
1311 :    
1312 : gdpusch 1.63 sub on_forward_strand {
1313 : overbeek 1.1 my ($self, $loc) = @_;
1314 :    
1315 : overbeek 1.12 warn "on_forward_strand: ", join(', ', @$loc) if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 2));
1316 : overbeek 1.26 return ($loc->[_STRAND] eq qq(+));
1317 : overbeek 1.1 }
1318 :    
1319 : gdpusch 1.63 sub on_reverse_strand {
1320 : overbeek 1.1 my ($self, $loc) = @_;
1321 :    
1322 : overbeek 1.12 warn "on_reverse_strand: ", join(', ', @$loc) if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 2));
1323 : overbeek 1.26 return ($loc->[_STRAND] eq qq(-));
1324 : overbeek 1.1 }
1325 :    
1326 : gdpusch 1.63 sub check_bounds {
1327 : overbeek 1.1 my ($self, $loc) = @_;
1328 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
1329 :    
1330 :     my $ok = 1;
1331 : gdpusch 1.63 foreach my $exon (@$loc) {
1332 : overbeek 1.1 my ($contig_id, $end, $strand, $len) = @$exon;
1333 : overbeek 1.12 print STDERR "Checking exon $contig_id, $end, $strand, $len\n"
1334 :     if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 2));
1335 : overbeek 1.1
1336 : overbeek 1.26 my $contig_length;
1337 : gdpusch 1.63 if (defined($contig_length = $self->get_contig_length($contig_id))) {
1338 :     if ($self->on_forward_strand($exon)) {
1339 : overbeek 1.1 my $beg = ($end - $len) + 1;
1340 : overbeek 1.7 if (($beg <= 0) || ($end > $contig_length)) {
1341 : overbeek 1.1 $ok = 0;
1342 : overbeek 1.14 print STDERR "\tOut-of-bounds plus-strand exon coordinates $contig_id\_$beg\_$end"
1343 : overbeek 1.12 , " (contig_length = $contig_length)\n" if $ENV{VERBOSE};
1344 : overbeek 1.1 }
1345 :     }
1346 : gdpusch 1.63 else {
1347 : overbeek 1.1 my $beg = ($end + $len) - 1;
1348 : overbeek 1.7 if (($end <= 0) || ($beg > $contig_length)) {
1349 : overbeek 1.1 $ok = 0;
1350 : overbeek 1.12 print STDERR "\tOut-of-bounds minus-strand exon coordinates $contig_id\_$beg\_$end"
1351 :     , " (contig_length = $contig_length)\n" if $ENV{VERBOSE};
1352 : overbeek 1.1 }
1353 :     }
1354 :     }
1355 :     }
1356 :    
1357 :     return $ok;
1358 :     }
1359 :    
1360 : gdpusch 1.63 sub get_dna_subseq {
1361 : overbeek 1.1 my ($self, $loc) = @_;
1362 : overbeek 1.26 my ($contig_seqP, $beg);
1363 : overbeek 1.1
1364 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
1365 :    
1366 :     my $dna_seq = "";
1367 : gdpusch 1.63 foreach my $exon (@$loc) {
1368 : overbeek 1.12 my ($contig_id, $end, $strand, $len) = @$exon;
1369 :    
1370 : gdpusch 1.63 if (not $self->check_bounds($exon)) {
1371 : overbeek 1.26 my $x = &flatten_dumper($loc);
1372 : overbeek 1.25 cluck "Exon $contig_id:$end$strand$len is out of bounds, loc: $x\n"
1373 : overbeek 1.14 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
1374 : overbeek 1.26 return undef;
1375 : overbeek 1.12 }
1376 : overbeek 1.1
1377 :     print STDERR "Getting exon $contig_id, $end, $strand, $len\n"
1378 : overbeek 1.12 if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 2));
1379 : overbeek 1.1
1380 : gdpusch 1.63 if (defined($self->get_contig_length($contig_id))) {
1381 : overbeek 1.26 $contig_seqP = $self->get_contig_seqP($contig_id);
1382 : overbeek 1.1
1383 : gdpusch 1.63 if ($self->on_forward_strand($exon)) {
1384 : overbeek 1.1 $beg = ($end - $len) + 1;
1385 : overbeek 1.26 $dna_seq .= substr($$contig_seqP, $beg-1, $len);
1386 : overbeek 1.1 }
1387 : gdpusch 1.63 else {
1388 : overbeek 1.1 $beg = ($end + $len) - 1;
1389 : overbeek 1.26 $dna_seq .= &FIG::reverse_comp(substr($$contig_seqP, $end-1, $len));
1390 : overbeek 1.1 }
1391 :     }
1392 : gdpusch 1.63 else {
1393 : overbeek 1.1 print STDERR "Invalid contig $contig_id\n";
1394 :     }
1395 :     }
1396 : overbeek 1.12 print STDERR "dna_seq = $dna_seq\n" if ($ENV{VERBOSE} && ($ENV{VERBOSE} > 2));
1397 : overbeek 1.1
1398 :     return $dna_seq;
1399 :     }
1400 :    
1401 : overbeek 1.26 sub flatten_dumper {
1402 : overbeek 1.25
1403 :     my $x = Dumper($_[0]);
1404 :    
1405 : overbeek 1.26 $x =~ s/\$VAR\d+\s+\=\s+//o;
1406 : overbeek 1.25 $x =~ s/\n//gso;
1407 :     $x =~ s/\s+/ /go;
1408 :     $x =~ s/\'//go;
1409 : overbeek 1.26 # $x =~ s/^[^\(\[\{]+//o;
1410 :     # $x =~ s/[^\)\]\}]+$//o;
1411 : overbeek 1.25
1412 :     return $x;
1413 :     }
1414 :    
1415 :    
1416 :    
1417 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1418 :     # ... Feature manipulation ...
1419 :     #-----------------------------------------------------------------------
1420 :    
1421 : gdpusch 1.63 sub delete_features_of_type {
1422 : overbeek 1.25 my ($self, @types) = @_;
1423 :    
1424 :     if ((not @types) || ($types[0] =~ m/all/o)) { @types = $self->get_feature_types(); }
1425 :     my $type = join(", ", @types);
1426 :    
1427 :     my @fids = $self->get_fids_for_type(@types);
1428 : overbeek 1.15 my $s = (@fids != 1) ? 's' : '' ;
1429 : overbeek 1.8 print STDERR "Preparing to delete ", (scalar @fids), " feature$s of type $type\n"
1430 :     if $ENV{VERBOSE};
1431 :    
1432 : gdpusch 1.63 foreach my $fid (@fids) {
1433 : overbeek 1.8 $self->delete_feature($fid) || confess "Could not delete feature $fid";
1434 :     }
1435 :    
1436 :     return 1;
1437 :     }
1438 :    
1439 : gdpusch 1.63 sub delete_feature {
1440 : overbeek 1.1 my ($self, $fid) = @_;
1441 : overbeek 1.8 print STDERR "Deleting feature $fid\n" if $ENV{VERBOSE};
1442 : overbeek 1.1
1443 : gdpusch 1.63 if (not $self->is_feature($fid)) {
1444 : overbeek 1.7 confess "Attempt to delete non-existent feature $fid";
1445 :     }
1446 :    
1447 :     my $contig_id = $self->get_feature_contig($fid);
1448 :     my $strand = $self->get_feature_strand($fid);
1449 :     my $end = $self->get_feature_endpoint($fid);
1450 :    
1451 : overbeek 1.12 # my $overlaps = $self->{_overlaps};
1452 :     # foreach my $x (keys % { $overlaps->{$fid} } )
1453 :     # {
1454 :     # if (delete $overlaps->{$x}->{$fid})
1455 :     # {
1456 :     # print STDERR " Deleting $x --> $fid from _overlap structure\n" if $ENV{VERBOSE};
1457 :     # }
1458 :     # else
1459 :     # {
1460 :     # confess "Could not delete $x --> $fid from _overlap structure";
1461 :     # }
1462 :     # }
1463 :     #
1464 :     # if (delete $overlaps->{$fid})
1465 :     # {
1466 :     # print STDERR "Deleting $fid from _overlap structure\n" if $ENV{VERBOSE};
1467 :     # }
1468 :     # else
1469 :     # {
1470 :     # confess "Could not delete $fid from _overlap structure";
1471 :     # }
1472 : overbeek 1.7
1473 : overbeek 1.15 my $key = "$contig_id$strand$end";
1474 :     my $used = $self->{_used_list}->{$key};
1475 :     my $ok = "";
1476 :     for (my $i=0; $i < @$used; ++$i) {
1477 :     if ($used->[$i] eq $fid) { $ok = splice(@$used, $i, 1); last; }
1478 :     }
1479 :    
1480 : gdpusch 1.63 if (not $ok) {
1481 : overbeek 1.15 confess "Could not delete $fid from _used structure, key $key, _used list = "
1482 :     , join(", ", @$used);
1483 : overbeek 1.7 }
1484 :    
1485 :    
1486 : overbeek 1.12 # my ($i, $sort);
1487 :     # $sort = $self->{_sort_left}->{$contig_id};
1488 :     # # print STDERR &Dumper($sort);
1489 :     # if (not @$sort) { confess "_sort_left structure for $contig_id is empty\n", &Dumper($self); }
1490 :     # for ($i=0; ($i < @$sort) && ($sort->[$i] ne $fid); $i++) {}
1491 :     # if ($i < @$sort) {
1492 :     # splice(@$sort, $i, 1) || confess "Could not delete $fid from _sort_right at i=$i";
1493 :     # } else {
1494 :     # confess "Could not delete $fid\n";
1495 :     # }
1496 :     #
1497 :     # $sort = $self->{_sort_right}->{$contig_id};
1498 :     # # print STDERR &Dumper($sort);
1499 :     # if (not @$sort) { confess "_sort_right structure for $contig_id is empty\n", &Dumper($self); }
1500 :     # for ($i=0; ($i < @$sort) && ($sort->[$i] ne $fid); $i++) {}
1501 :     # if ($i < @$sort) {
1502 :     # splice(@$sort, $i, 1) || confess "Could not delete $fid from _sort_right at i=$i";
1503 :     # } else {
1504 :     # confess "Could not delete $fid\n";
1505 :     # }
1506 :    
1507 :    
1508 :     if (not delete $self->{_features}->{$fid})
1509 :     {
1510 : overbeek 1.7 confess "Could not delete $fid from _features structure";
1511 :     }
1512 :    
1513 :     return 1;
1514 : overbeek 1.1 }
1515 :    
1516 : gdpusch 1.63 sub add_feature {
1517 : overbeek 1.7 my ($self, %args) = @_;
1518 :    
1519 : overbeek 1.31 my $type = $args{-type} || confess "No feature type given";
1520 : overbeek 1.1
1521 : overbeek 1.19 my $fid = $args{-fid};
1522 :     my $func = $args{-func};
1523 :     my $seq = $args{-seq};
1524 :     my $aliases = $args{-aliases};
1525 : gdpusch 1.74 my $annot = $args{-annot};
1526 : overbeek 1.7
1527 : overbeek 1.19 my $tax_id = $self->get_taxid;
1528 : overbeek 1.31 my $base = "fig|$tax_id.$type";
1529 : overbeek 1.7
1530 : overbeek 1.31 my $loc = $args{-loc} || confess "No feature loc given";
1531 :     if ($ENV{VERBOSE}) {
1532 :     print STDERR qq(Attempting to add '$type' feature);
1533 : overbeek 1.32 if ($ENV{VERBOSE} > 2) {
1534 :     print STDERR qq(\n args = ) . &flatten_dumper(\%args) . qq(\n);
1535 :     }
1536 :     else {
1537 : overbeek 1.31 print STDERR qq(, loc=) . &flatten_dumper($loc) . qq(\n);
1538 :     }
1539 :     }
1540 : gdpusch 1.74
1541 :     if (ref($loc) ne 'ARRAY') {
1542 :     confess "Non-ARRAYref loc object passed to add_feature --- \%args =\n"
1543 :     , Dumper(\%args);
1544 :     }
1545 :    
1546 : overbeek 1.31 if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
1547 :     my ($contig_id, $end, $strand, $len) = @ { $loc->[-1] };
1548 : gdpusch 1.55
1549 :     $len = 0;
1550 : gdpusch 1.63 foreach my $exon (@$loc) {
1551 : gdpusch 1.55 $len += $self->get_exon_length($exon);
1552 :     }
1553 :    
1554 :     if ( (($type eq 'peg') || ($type eq 'orf')) && (($len % 3) != 0) ) {
1555 :     warn "Triality test fails for loc: " . &flatten_dumper($loc) . "\n";
1556 :     return undef;
1557 :     }
1558 : overbeek 1.1
1559 : gdpusch 1.63 if (not defined($self->{_features})) {
1560 : overbeek 1.1 $self->{_features} = {};
1561 : overbeek 1.31 $self->{_features}->{_maxnum} = {};
1562 : overbeek 1.4 print STDERR "No previous features --- creating feature-pointer\n" if $ENV{VERBOSE};
1563 : overbeek 1.1 }
1564 :    
1565 : gdpusch 1.63 if (not defined($self->{_features}->{_maxnum}->{$type})) {
1566 : overbeek 1.1 $self->{_features}->{_maxnum}->{$type} = 0;
1567 :     print STDERR "No previous features of type $type\n" if $ENV{VERBOSE};
1568 :     }
1569 : overbeek 1.7 my $numref = \$self->{_features}->{_maxnum}->{$type};
1570 : overbeek 1.1
1571 : overbeek 1.31 if (not $fid) {
1572 :     $fid = $base . qq(.) . ++$$numref;
1573 :     }
1574 :     else {
1575 :     if ($fid =~ m/\.(\d+)$/o) {
1576 :     $$numref = &FIG::max($$numref, $1);
1577 :     }
1578 :     else {
1579 :     confess "Malformed FID $fid";
1580 :     }
1581 :     }
1582 :    
1583 : gdpusch 1.63 if (defined($self->{_features}->{$fid})) {
1584 : overbeek 1.31 confess "FATAL ERROR: Redefining existing feature $fid";
1585 :     }
1586 : overbeek 1.12
1587 : gdpusch 1.63 if ($contig_id && $end && $strand && $len) {
1588 : overbeek 1.1 my $new_feature = {};
1589 :    
1590 :     $new_feature->{_type} = $type;
1591 :     $new_feature->{_loc} = $loc;
1592 :    
1593 : overbeek 1.5 if (defined($func)) { $new_feature->{_func} = $func; }
1594 : overbeek 1.1
1595 : gdpusch 1.63 if ($seq) {
1596 : overbeek 1.6 $new_feature->{_seq} = $seq;
1597 :     }
1598 : overbeek 1.38 else {
1599 :     if (($type eq 'peg') || ($type eq 'orf')) {
1600 : overbeek 1.39 my $tmp_loc = $self->copy_loc($loc);
1601 : overbeek 1.38 if ($type eq 'orf') {
1602 :     my $orf_len = $self->get_orf_length($tmp_loc->[-1]);
1603 :     $self->set_exon_length($tmp_loc->[-1], $orf_len)
1604 :     || confess "For FID $fid, Could not set len=orf_len=$orf_len for tmp_loc="
1605 :     , &flatten_dumper($tmp_loc);
1606 :     }
1607 :    
1608 : overbeek 1.40 my $dna = $self->get_dna_subseq($tmp_loc);
1609 : gdpusch 1.49 if (my $pep = &FIG::translate($dna, $self->get_translation_table, 1)) {
1610 : overbeek 1.38 $pep =~ s/\*$//o;
1611 : gdpusch 1.55 if ($pep =~ m/\*/o) {
1612 :     warn "Translation of $fid contains internal STOPs --- skipping\n";
1613 :     return undef;
1614 :     }
1615 : overbeek 1.38 $new_feature->{_seq} = $pep;
1616 : gdpusch 1.63 }
1617 :     else {
1618 : overbeek 1.38 confess "Translation failed for $dna";
1619 :     }
1620 :     }
1621 :     else {
1622 :     $new_feature->{_seq} = $self->get_dna_subseq($loc);
1623 : overbeek 1.1 }
1624 :     }
1625 :    
1626 : overbeek 1.19 if (defined($aliases)) {
1627 :     if (ref($aliases) eq 'ARRAY') {
1628 :     $new_feature->{_aliases} = $aliases;
1629 :     }
1630 :     else {
1631 : overbeek 1.26 confess "-aliases field is not an ARRAY ref for feature: "
1632 :     . &flatten_dumper($aliases, $new_feature) . "\n";
1633 : overbeek 1.19 }
1634 :     }
1635 :    
1636 : gdpusch 1.74
1637 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1638 :     #... This is a crude hack, since it only allows immediately writing
1639 :     # an annotation.
1640 :     #... Properly implemented, annotations should be handled as a queue,
1641 :     # with proper access methods.
1642 :     #-----------------------------------------------------------------------
1643 :     if (defined($annot)) {
1644 :     if (ref($annot) eq 'ARRAY') {
1645 :     my $timestamp = time();
1646 :     my $annot_file = $self->get_genome_dir() . qq(/annotations);
1647 :     my $annot_user = $annot->[0];
1648 :     my $annot_text = $annot->[1];
1649 :    
1650 :     open( ANNOT, ">>$annot_file")
1651 :     || confess "Could not append-open $annot_file";
1652 :    
1653 :     print ANNOT "$fid\n$timestamp\n$annot_user\n$annot_text\n//\n";
1654 :    
1655 :     close(ANNOT)
1656 :     || confess "Could not close $annot_file";
1657 :     }
1658 :     else {
1659 :     confess "-annot field is not an ARRAY ref for feature: "
1660 :     . &flatten_dumper($annot, $new_feature) . "\n";
1661 :     }
1662 :     }
1663 :    
1664 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1665 :     #...And finally, store the new feature in the NewGenome object struct:
1666 :     #-----------------------------------------------------------------------
1667 : overbeek 1.1 $self->{_features}->{$fid} = $new_feature;
1668 : gdpusch 1.74 #=======================================================================
1669 :    
1670 : overbeek 1.1 }
1671 : gdpusch 1.63 else {
1672 : gdpusch 1.74 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1673 :     #...Object creation is being aborted because one or more of
1674 :     # $contig_id, $end, $strand, or $len are missing or nil;
1675 :     # signal failure by returning 'undef'.
1676 :     #------------------------------------------------------------------------
1677 : overbeek 1.1 $fid = undef;
1678 :     }
1679 : overbeek 1.10
1680 : overbeek 1.26 if ($ENV{VERBOSE}) {
1681 :     if ($ENV{VERBOSE} == 1) {
1682 :     print STDERR "Added feature $fid at $contig_id:$end$strand$len"
1683 :     , (($#{$loc->[0]} == 4) ? (defined($loc->[0]->[4]) ? ",$loc->[0]->[4]" : ',undef') : "")
1684 :     , " func=$func"
1685 :     , "\n";
1686 :     }
1687 :     else {
1688 :     print STDERR "Added feature $fid at $contig_id:$end$strand$len:\n "
1689 :     , &flatten_dumper($self->{_features}->{$fid}), "\n\n";
1690 :     }
1691 :     }
1692 : overbeek 1.7
1693 :     my $key = "$contig_id$strand$end";
1694 : overbeek 1.15 if (not defined($self->{_used_list}->{$key})) { $self->{_used_list}->{$key} = []; }
1695 :     push @ { $self->{_used_list}->{$key} }, $fid;
1696 : overbeek 1.7
1697 : overbeek 1.12 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1698 :     #... Update dynamic overlap structures ...
1699 :     #-----------------------------------------------------------------------
1700 :     # my $sort;
1701 :     #
1702 :     # my $mid = 0;
1703 :     # if (not defined($self->{_sort_left}->{$contig_id}))
1704 :     # {
1705 :     # $self->{_sort_left}->{$contig_id} = [$fid];
1706 :     # }
1707 :     # else
1708 :     # {
1709 :     # $sort = $self->{_sort_left}->{$contig_id};
1710 :     #
1711 :     # my $low = 0;
1712 :     # my $high = $#$sort;
1713 :     # my $mid;
1714 :     # my $cmp;
1715 :     #
1716 :     # while ($low < $high)
1717 :     # {
1718 :     # $mid = int( ($low+$high) / 2 );
1719 :     #
1720 :     # $cmp = $self->compare_left($fid, $sort->[$mid]);
1721 :     # # print STDERR "low=$low, high=$high, mid=$mid, cmp=$cmp\n";
1722 :     #
1723 :     # if ( $cmp > 0 ) { $low = $mid+1; } else { $high = $mid; }
1724 :     # }
1725 :     #
1726 :     # $cmp = $self->compare_left($fid,$sort->[$high]);
1727 :     # if ($cmp == 0)
1728 :     # {
1729 :     # # print STDERR "IT IS AT $high\n";
1730 :     # splice @$sort, $high, 0, $fid;
1731 :     # }
1732 :     # elsif ($cmp > 0)
1733 :     # {
1734 :     # # print STDERR "INSERT after $high\n";
1735 :     # splice @$sort, $high+1, 0, $fid;
1736 :     # }
1737 :     # else
1738 :     # {
1739 :     # # print STDERR "INSERT before $high\n";
1740 :     # splice @$sort, $high-1, 0, $fid;
1741 :     # }
1742 :     # }
1743 :     #
1744 :     # $mid = 0;
1745 :     # if (not defined($self->{_sort_right}->{$contig_id}))
1746 :     # {
1747 :     # $self->{_sort_right}->{$contig_id} = [$fid];
1748 :     # }
1749 :     # else
1750 :     # {
1751 :     # $sort = $self->{_sort_right}->{$contig_id};
1752 :     #
1753 :     # my $low = 0;
1754 :     # my $high = $#$sort;
1755 :     # my $mid;
1756 :     # my $cmp;
1757 :     #
1758 :     # while ($low < $high)
1759 :     # {
1760 :     # $mid = int( ($low+$high) / 2 );
1761 :     #
1762 :     # $cmp = $self->compare_right($fid, $sort->[$mid]);
1763 :     # # print STDERR "low=$low, high=$high, mid=$mid, cmp=$cmp\n";
1764 :     #
1765 :     # if ( $cmp > 0 ) { $low = $mid+1; } else { $high = $mid; }
1766 :     # }
1767 :     #
1768 :     # $cmp = $self->compare_right($fid,$sort->[$high]);
1769 :     # if ($cmp == 0)
1770 :     # {
1771 :     # # print STDERR "IT IS AT $high\n";
1772 :     # splice @$sort, $high, 0, $fid;
1773 :     # }
1774 :     # elsif ($cmp > 0)
1775 : overbeek 1.7 # {
1776 : overbeek 1.12 # # print STDERR "INSERT after $high\n";
1777 :     # splice @$sort, $high+1, 0, $fid;
1778 : overbeek 1.7 # }
1779 : overbeek 1.12 # else
1780 : overbeek 1.8 # {
1781 : overbeek 1.12 # # print STDERR "INSERT before $high\n";
1782 :     # splice @$sort, $high-1, 0, $fid;
1783 : overbeek 1.8 # }
1784 : overbeek 1.12 # }
1785 :     #
1786 :     # my $overlaps = $self->{_overlaps};
1787 :     # if (not defined($overlaps->{$sort->[$mid]})) { $overlaps->{$sort->[$mid]} = {}; }
1788 :     #
1789 :     # $sort = $self->{_sort_left}->{$contig_id};
1790 :     # for (my $i=$mid+1; $i <= $#$sort; ++$i)
1791 :     # {
1792 : overbeek 1.28 # my $lap = $self->overlap_between($sort->[$mid], $sort->[$i]);
1793 : overbeek 1.12 # last unless $lap;
1794 :     #
1795 :     # if (not defined($overlaps->{$sort->[$i]})) { $overlaps->{$sort->[$i]} = {}; }
1796 :     #
1797 :     # $overlaps->{$sort->[$mid]}->{$sort->[$i]} = $lap;
1798 :     # $overlaps->{$sort->[$i]}->{$sort->[$mid]} = $lap;
1799 :     #
1800 :     # print STDERR "Adding overlaps between $sort->[$mid] <--> $sort->[$i]\n"
1801 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
1802 :     # }
1803 :     #
1804 :     # $sort = $self->{_sort_right}->{$contig_id};
1805 :     # for (my $i=$mid-1; $i >= 0; --$i)
1806 :     # {
1807 : overbeek 1.28 # my $lap = $self->overlap_between($sort->[$mid], $sort->[$i]);
1808 : overbeek 1.12 # last unless $lap;
1809 :     #
1810 :     # if (not defined($overlaps->{$sort->[$i]})) { $overlaps->{$sort->[$i]} = {}; }
1811 :     #
1812 :     # $overlaps->{$sort->[$mid]}->{$sort->[$i]} = $lap;
1813 :     # $overlaps->{$sort->[$i]}->{$sort->[$mid]} = $lap;
1814 :     #
1815 :     # print STDERR "Adding overlaps between $sort->[$mid] <--> $sort->[$i]\n"
1816 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
1817 :     # }
1818 : overbeek 1.1
1819 :     return $fid;
1820 :     }
1821 :    
1822 : gdpusch 1.63 sub remove_overlaps {
1823 : overbeek 1.12 my ($self, $closest_org, $min_len) = @_;
1824 :     my $dir = $self->get_genome_dir();
1825 : overbeek 1.10
1826 : overbeek 1.12 my $rna = (-s "$dir/Features/rna/tbl") ? "$dir/Features/rna/tbl" : '/dev/null';
1827 :     my $peg = (-s "$dir/Features/peg/tbl") ? "$dir/Features/peg/tbl" : '/dev/null';
1828 :     my $orf = (-s "$dir/Features/orf/tbl") ? "$dir/Features/orf/tbl" : '/dev/null';
1829 : overbeek 1.42 print STDERR "In remove_overlaps:\n rna=$rna\n peg=$peg\n orf=$orf\n"
1830 :     if $ENV{VERBOSE};
1831 : overbeek 1.10
1832 : overbeek 1.12 $self->export_features('all');
1833 : overbeek 1.10
1834 : overbeek 1.12 my @tbl;
1835 : overbeek 1.42 if ($ENV{VERBOSE}) {
1836 : olson 1.44 @tbl = `$FIG_Config::bin/filter_overlaps 90 20 50 150 150 $rna $peg < $orf 2>> $dir/overlap_removal.log`;
1837 : overbeek 1.42 }
1838 :     else {
1839 : olson 1.44 @tbl = `$FIG_Config::bin/filter_overlaps 90 20 50 150 150 $rna $peg < $orf`;
1840 : overbeek 1.10 }
1841 : overbeek 1.17 print STDERR "Keeping ", (scalar @tbl), " ORFs\n" if $ENV{VERBOSE};
1842 :     # die "aborting";
1843 : overbeek 1.10
1844 : overbeek 1.31 my %keep = map { m/^(\S+)/o; $1 => 1 } @tbl;
1845 : overbeek 1.12
1846 :     my $num_deleted = 0;
1847 : gdpusch 1.63 foreach my $fid ($self->get_fids_for_type('orf')) {
1848 :     if (not $keep{$fid}) {
1849 : overbeek 1.17 print STDERR "Attempting to delete $fid\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
1850 : overbeek 1.12 $self->delete_feature($fid) || confess "Could not delete $fid";
1851 :     ++$num_deleted;
1852 : overbeek 1.10 }
1853 :     }
1854 :    
1855 : overbeek 1.12 $self->export_features('all');
1856 :    
1857 : overbeek 1.10 return $num_deleted;
1858 :     }
1859 :    
1860 : overbeek 1.12 # sub remove_overlapping_fids
1861 :     # {
1862 :     # my ($self, %args) = @_;
1863 :     # my ($type, $olap, $type2);
1864 :     #
1865 :     # my $num_deleted = 0;
1866 :     #
1867 :     # my $fid = $args{-fid} || confess "No FID given to ".__PACKAGE__."::remove_overlaps";
1868 :     # my $rna_maxlap = $args{-rna} || 20;
1869 :     #
1870 :     # my @overlapping_fids = $self->get_overlapping_fids($fid);
1871 :     # return (0) unless @overlapping_fids;
1872 :     #
1873 :     # $type = $self->get_feature_type($fid);
1874 :     #
1875 :     # if ($type eq 'rna')
1876 :     # {
1877 :     # foreach my $fid2 (@overlapping_fids)
1878 :     # {
1879 :     # $type2 = $self->get_feature_type($fid2);
1880 :     # $olap = $self->get_overlap($fid, $fid2);
1881 :     #
1882 :     # if (($olap > $rna_maxlap) && (($type2 eq 'peg') || ($type2 eq 'orf')))
1883 :     # {
1884 :     # $self->delete_feature($fid2) || confess "Could not delete $fid2 overlapping $fid by $olap";
1885 :     # ++$num_deleted;
1886 :     # }
1887 :     # }
1888 :     #
1889 :     # return $num_deleted;
1890 :     # }
1891 :     #
1892 :     # if ($type eq 'peg')
1893 :     # {
1894 :     # foreach my $fid2 (@overlapping_fids)
1895 :     # {
1896 :     # $type2 = $self->get_feature_type($fid2);
1897 : overbeek 1.28 # if ($self->contains_embedded($fid, $fid2) && ($type2 eq 'orf'))
1898 : overbeek 1.12 # {
1899 :     # $self->delete_feature($fid2) || confess "Could not delete $fid2 overlapping $fid by $olap";
1900 :     # ++$num_deleted;
1901 :     # }
1902 :     # }
1903 :     # }
1904 :     #
1905 :     # return $num_deleted;
1906 :     # }
1907 :    
1908 :    
1909 : gdpusch 1.63 sub find_rna {
1910 : overbeek 1.6 my ($self) = @_;
1911 :    
1912 : olson 1.58 my $default_tool = "/vol/search_for_rnas-2007-0625/search_for_rnas";
1913 :     # my $default_tool = "/vol/search_for_rnas/bin/search_for_rnas"; # old version
1914 :     # my $default_tool = "$ENV{HOME}/rna_search/search_for_rnas";
1915 : olson 1.56
1916 :     my $tool = $ENV{RP_SEARCH_FOR_RNAS};
1917 : gdpusch 1.63 if (!$tool) {
1918 : olson 1.58 $tool = $default_tool;
1919 : olson 1.56 }
1920 :    
1921 : gdpusch 1.63 if (-x $tool) {
1922 : overbeek 1.6 my $taxid = $self->get_taxid();
1923 :     my $dir = $self->get_genome_dir();
1924 : olson 1.57
1925 : overbeek 1.6 my $bioname = $self->get_genome_name();
1926 :     my ($genus, $species);
1927 : overbeek 1.31 if ($bioname =~ m/^(\S+)\s+(\S+)/o) {
1928 : overbeek 1.6 ($genus, $species) = ($1, $2);
1929 :     } else {
1930 :     confess "Could not parse bioname $bioname";
1931 :     }
1932 :    
1933 :     my $taxonomy = $self->get_taxonomy();
1934 :     my $domain = uc(substr($taxonomy, 0, 1));
1935 :    
1936 : olson 1.57 my $cmd;
1937 : overbeek 1.19 if ($ENV{VERBOSE}) {
1938 : olson 1.57 $cmd = "$tool --contigs=$dir/contigs --orgid=$taxid --domain=$domain --genus=$genus --species=$species --log=$dir/rna.log";
1939 :     warn "Finding RNAs with $cmd\n";
1940 : overbeek 1.19 }
1941 :     else {
1942 : olson 1.57 $cmd = "$tool --contigs=$dir/contigs --orgid=$taxid --domain=$domain --genus=$genus --species=$species --log=$dir/rna.log 2> /dev/null";
1943 :     }
1944 :     my @tbl = `$cmd`;
1945 : gdpusch 1.63 if ($?) {
1946 : olson 1.57 my($rc, $sig, $msg) = &FIG::interpret_error_code($?);
1947 :     confess "$msg: $cmd";
1948 : overbeek 1.19 }
1949 : overbeek 1.6
1950 :     return @tbl;
1951 :     }
1952 : gdpusch 1.63 else {
1953 : overbeek 1.19 cluck "No RNA-finder tool found at $tool\n";
1954 : overbeek 1.6 return ();
1955 :     }
1956 :     }
1957 :    
1958 : gdpusch 1.63 sub possible_orfs {
1959 : overbeek 1.1 my ($self, $min_len) = @_;
1960 :     my $class = ref($self);
1961 :     my $taxon = $self->get_taxid;
1962 :    
1963 : overbeek 1.5 if (not $min_len)
1964 :     {
1965 :     $min_len = 90;
1966 :     cluck "Setting min_len = $min_len bp by default \n" if $ENV{VERBOSE};
1967 :     }
1968 :    
1969 : gdpusch 1.49 print STDERR "Beginning call to GLIMMER2\n" if $ENV{VERBOSE};
1970 : overbeek 1.5
1971 : gdpusch 1.49 my $org = $self->get_taxid();
1972 :     my $dir = $self->get_genome_dir();
1973 : overbeek 1.5 my @tmp_tbl;
1974 : gdpusch 1.49 my $code = $self->get_genetic_code_number();
1975 :     my $cmd = "$FIG_Config::bin/run_glimmer2 $org $dir/contigs -code=$code";
1976 : olson 1.44 my $err_tmp = "$FIG_Config::temp/glimmer.$$.stderr";
1977 :    
1978 :     open(GLIM, "$cmd 2> $err_tmp |") or confess "Error opening glimmer pipe $cmd: $!\n";
1979 :    
1980 : gdpusch 1.63 while (<GLIM>) {
1981 : olson 1.44 push(@tmp_tbl, $_);
1982 :     }
1983 :    
1984 : gdpusch 1.63 if (!close(GLIM)) {
1985 : olson 1.44 open(ERR, "<$err_tmp");
1986 : gdpusch 1.63 while (<ERR>) {
1987 : olson 1.44 print STDERR $_;
1988 :     }
1989 :     close(ERR);
1990 :     confess "Glimmer pipe close failed \$?=$? \$!=$!\n";
1991 :     }
1992 : overbeek 1.5
1993 : overbeek 1.38 my ($loc, $tmp_loc);
1994 : gdpusch 1.63 foreach my $entry (@tmp_tbl) {
1995 : overbeek 1.5 chomp $entry;
1996 :     my (undef, $locus) = split /\t/, $entry;
1997 :    
1998 : gdpusch 1.63 if ($locus =~ m/(\S+)_(\d+)_(\d+)/o) {
1999 : overbeek 1.5 my ($contig_id, $beg, $end) = ($1, $2, $3);
2000 : overbeek 1.31 my $len = 1 + abs($end - $beg);
2001 : overbeek 1.26 my $strand = ($end > $beg) ? qq(+) : qq(-) ;
2002 : overbeek 1.5
2003 : overbeek 1.31 $loc = $self->make_loc($contig_id, $end, $strand, $len)
2004 :     || confess "Could not create loc=$contig_id:$end$strand$len";
2005 :     print STDERR "Created loc=", &flatten_dumper($loc), "\n"
2006 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2007 :    
2008 :     $tmp_loc = $self->search_for_upstream_stop($loc)
2009 :     || confess "Could not find ORF boundary for loc=", &flatten_dumper($loc);
2010 :     print STDERR "ORF upstream STOP loc=", &flatten_dumper($tmp_loc), "\n"
2011 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2012 :    
2013 : overbeek 1.38 # my $orf_len = $self->get_orf_length($tmp_loc)
2014 :     # || confess "Could not extract ORF length from tmp_loc=", &flatten_dumper($tmp_loc);
2015 :     #
2016 :     # $self->set_exon_length($tmp_loc, $orf_len)
2017 :     # || confess "Could not set length field to ORF length for loc=", &flatten_dumper($tmp_loc);
2018 :     # print STDERR "Reset len=$orf_len, tmp_loc=", &flatten_dumper($tmp_loc), "\n"
2019 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2020 :     #
2021 :     # $tmp_loc = $self->search_for_downstream_start($tmp_loc)
2022 :     # || confess "Could not find first START in ORF loc=", &flatten_dumper($tmp_loc);
2023 :     # print STDERR "First START is at tmp_loc=", &flatten_dumper($tmp_loc), "\n"
2024 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2025 : overbeek 1.14
2026 : overbeek 1.31 $self->add_feature( -type => 'orf', -loc => $tmp_loc )
2027 :     || confess "Could not add 'orf' feature at loc=", &flatten_dumper($tmp_loc);
2028 : overbeek 1.5 }
2029 : overbeek 1.31 else {
2030 : overbeek 1.5 confess "Could not parse ORF header: $locus";
2031 :     }
2032 :     }
2033 :     print STDERR "Completed call to GLIMMER2 --- returning ORFs\n" if $ENV{VERBOSE};
2034 :    
2035 :     return { map { $_ => $self->{_features}->{$_} } $self->get_fids_for_type('orf') };
2036 :     }
2037 :    
2038 : gdpusch 1.63 sub recall_orfs {
2039 : overbeek 1.12 my ($self, $closest_org, $min_len) = @_;
2040 : overbeek 1.7 my $class = ref($self);
2041 :     my $taxon = $self->get_taxid;
2042 :    
2043 : gdpusch 1.63 if (not $min_len) {
2044 : overbeek 1.7 $min_len = 90;
2045 : overbeek 1.8 # cluck "Setting min_len = $min_len bp by default \n" if $ENV{VERBOSE};
2046 : overbeek 1.7 }
2047 :    
2048 : overbeek 1.26 $self->delete_features_of_type('orf') || confess "Could not delete existing ORFs";
2049 : overbeek 1.8
2050 : overbeek 1.7 print STDERR "Begining call to GLIMMER2\n" if $ENV{VERBOSE};
2051 : overbeek 1.13 print STDERR ("pwd=", `pwd`) if $ENV{VERBOSE};
2052 : overbeek 1.7
2053 : gdpusch 1.49 my $org = $self->get_taxid();
2054 :     my $dir = $self->get_genome_dir();
2055 :     my $code = $self->get_genetic_code_number();
2056 : overbeek 1.8 my $peg_tbl = "$dir/Features/peg/tbl";
2057 :     confess "Empty or non-existent training data $peg_tbl, size=", (-s $peg_tbl)
2058 :     unless (-s "$dir/Features/peg/tbl");
2059 :    
2060 : overbeek 1.7 my @tmp_tbl;
2061 : gdpusch 1.63 if ($ENV{VERBOSE}) {
2062 : gdpusch 1.49 (@tmp_tbl = `$FIG_Config::bin/run_glimmer2 $org $dir/contigs -train=$peg_tbl -code=$code`)
2063 : overbeek 1.8 || confess "GLIMMER found no candidate ORFs in $dir/contigs";
2064 :     }
2065 : gdpusch 1.63 else {
2066 : gdpusch 1.49 (@tmp_tbl = `$FIG_Config::bin/run_glimmer2 $org $dir/contigs -train=$peg_tbl -code=$code 2> /dev/null`)
2067 : overbeek 1.8 || confess "GLIMMER found no candidate ORFs in $dir/contigs";
2068 :     }
2069 : overbeek 1.26
2070 : overbeek 1.38 my ($loc, $tmp_loc);
2071 : gdpusch 1.63 foreach my $entry (@tmp_tbl) {
2072 : overbeek 1.7 chomp $entry;
2073 :     my (undef, $locus) = split /\t/, $entry;
2074 :    
2075 : gdpusch 1.63 if ($locus =~ m/(\S+)_(\d+)_(\d+)/o) {
2076 : overbeek 1.7 my ($contig_id, $beg, $end) = ($1, $2, $3);
2077 : overbeek 1.26 my $len = 1 + abs($end - $beg);
2078 :     my $strand = ($end > $beg) ? qq(+) : qq(-) ;
2079 : overbeek 1.7
2080 : overbeek 1.31 $loc = $self->make_loc($contig_id, $end, $strand, $len)
2081 :     || confess "Could not create loc=$contig_id:$end$strand$len";
2082 :     print STDERR "Created loc=", &flatten_dumper($loc), "\n"
2083 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2084 :    
2085 :     $tmp_loc = $self->search_for_upstream_stop($loc)
2086 :     || confess "Could not find ORF boundary for loc=", &flatten_dumper($loc);
2087 :     print STDERR "ORF upstream STOP loc=", &flatten_dumper($tmp_loc), "\n"
2088 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2089 :    
2090 : overbeek 1.38 # my $orf_len = $self->get_orf_length($tmp_loc)
2091 :     # || confess "Could not extract ORF length from tmp_loc=", &flatten_dumper($tmp_loc);
2092 :     #
2093 :     # $self->set_exon_length($tmp_loc, $orf_len)
2094 :     # || confess "Could not set length field to ORF length for loc=", &flatten_dumper($tmp_loc);
2095 :     # print STDERR "Reset len=$orf_len, tmp_loc=", &flatten_dumper($tmp_loc), "\n"
2096 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2097 :     #
2098 :     # $tmp_loc = $self->search_for_downstream_start($tmp_loc)
2099 :     # || confess "Could not find first START in ORF loc=", &flatten_dumper($tmp_loc);
2100 :     # print STDERR "First START is at tmp_loc=", &flatten_dumper($tmp_loc), "\n"
2101 :     # if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2102 : overbeek 1.31
2103 :     $self->add_feature( -type => 'orf', -loc => $tmp_loc )
2104 :     || confess "Could not add 'orf' feature at loc=", &flatten_dumper($tmp_loc);
2105 : overbeek 1.7 }
2106 : overbeek 1.31 else {
2107 : overbeek 1.7 confess "Could not parse ORF header: $locus";
2108 :     }
2109 :     }
2110 : overbeek 1.17 print STDERR "Found ", (scalar $self->get_fids_for_type('orf')), " new ORFs\n" if $ENV{VERBOSE};
2111 :    
2112 : overbeek 1.12 print STDERR "Completed call to GLIMMER2 --- removing overlaps\n" if $ENV{VERBOSE};
2113 :     $self->remove_overlaps($closest_org, $min_len);
2114 : overbeek 1.14
2115 : overbeek 1.7 return { map { $_ => $self->{_features}->{$_} } $self->get_fids_for_type('orf') };
2116 :     }
2117 :    
2118 : gdpusch 1.63 sub pull_orfs {
2119 : overbeek 1.5 my ($self, $min_len) = @_;
2120 :     my $class = ref($self);
2121 :     my $taxon = $self->get_taxid;
2122 :    
2123 : overbeek 1.1 my $tmp_orfs = "$FIG_Config::temp/tmp$$\_$class.fasta";
2124 :    
2125 : gdpusch 1.63 if (not $min_len) {
2126 : overbeek 1.2 $min_len = 90;
2127 : gdpusch 1.63 # cluck "Setting min_len = $min_len bp by default \n" if $ENV{VERBOSE};
2128 : overbeek 1.2 }
2129 :    
2130 : overbeek 1.1 open(PULL, "| pull_orfs $min_len > $tmp_orfs")
2131 :     || die "Could not open pipe-out through pull_orfs $min_len > $tmp_orfs";
2132 :    
2133 : overbeek 1.26 my $contig_seqP;
2134 : gdpusch 1.63 foreach my $contig_id ( $self->get_contig_ids ) {
2135 : overbeek 1.2 print STDERR "Pulling ORFs for contig $contig_id ...\n" if $ENV{VERBOSE};
2136 : gdpusch 1.63 if (defined($contig_seqP = $self->get_contig_seqP($contig_id))) {
2137 : overbeek 1.26 print PULL ">$contig_id\n$$contig_seqP\n";
2138 : overbeek 1.1 }
2139 :     }
2140 :     close(PULL) || die "Could not close pipe-out through pull_orfs $min_len > $tmp_orfs";
2141 : overbeek 1.2 confess "No ORFs found for ".$self->get_genome_dir() if (not -s $tmp_orfs);
2142 :    
2143 : overbeek 1.1 open(ORFS, "<$tmp_orfs") || die "Could not read-open $tmp_orfs";
2144 : gdpusch 1.63 while (my ($locus, $transP) = &FIG::read_fasta_record(\*ORFS)) {
2145 : overbeek 1.12 print STDERR "Adding ORF $locus\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2146 :    
2147 : gdpusch 1.63 if ($locus =~ m/(\S+)_(\d+)_(\d+)/o) {
2148 : overbeek 1.1 my ($contig_id, $beg, $end) = ($1, $2, $3);
2149 : overbeek 1.12 my $contig_len = $self->get_contig_length($contig_id);
2150 : overbeek 1.26 my $strand = ($end > $beg) ? qq(+) : qq(-) ;
2151 :     my $sign = ($end > $beg) ? +1 : -1 ;
2152 : overbeek 1.12 if ( (($end + $sign*3) > 0) && (($end + $sign*3) <= $contig_len) ) { $end += $sign*3; }
2153 : overbeek 1.1 my $orf_len = 1 + abs($end - $beg);
2154 : overbeek 1.12
2155 :    
2156 :     my $codon = $self->get_dna_subseq([$contig_id, $end, $strand, 3]);
2157 : gdpusch 1.63 if (&FIG::translate($codon, $self->get_translation_table) ne '*') {
2158 : overbeek 1.12 print STDERR " WARNING: $contig_id:$end$strand$orf_len does not end with STOP ($codon)\n"
2159 :     if $ENV{VERBOSE};
2160 :     }
2161 :    
2162 :     my $b;
2163 :     if ( (($beg-$sign*3) > 0) && (($beg-$sign*3) <= $contig_len) ) {
2164 :     $b = $beg - $sign;
2165 :     } else {
2166 :     $b = $beg + $sign*2;
2167 :     }
2168 :     $codon = $self->get_dna_subseq([$contig_id, $b, $strand, 3]);
2169 : gdpusch 1.63 if (&FIG::translate($codon, $self->get_translation_table) ne '*') {
2170 : overbeek 1.12 print STDERR " WARNING: $contig_id:$end$strand$orf_len does not begin with STOP ($codon)\n"
2171 :     if $ENV{VERBOSE};
2172 :     }
2173 : overbeek 1.1
2174 : overbeek 1.7 $self->add_feature( -type => 'orf'
2175 :     , -loc => [$contig_id, $end, $strand, $orf_len, $orf_len]
2176 :     , -seq => $$transP)
2177 : overbeek 1.1 || confess "Could not add 'orf' feature at $contig_id:$end$strand$orf_len";
2178 :     }
2179 : gdpusch 1.63 else {
2180 : overbeek 1.1 confess "Could not parse ORF header: $locus";
2181 :     }
2182 :     }
2183 :     unlink $tmp_orfs or confess "Could not remove $tmp_orfs";
2184 :    
2185 :     return { map { $_ => $self->{_features}->{$_} } $self->get_fids_for_type('orf') };
2186 :     }
2187 :    
2188 : olson 1.51 #
2189 :     # Blast the given query sequence against a database consisting of
2190 :     # the ORFs or pegs of this genome.
2191 :     #
2192 : gdpusch 1.63 sub candidate_orfs {
2193 : overbeek 1.1 my ($self, %args) = @_;
2194 :    
2195 : overbeek 1.10 my $query_id = $args{-qid} || 'query_seq';
2196 :     my $query_seq = $args{-seq} || $fig->get_translation($args{-qid}) || confess "No query sequence provided";
2197 :     my $cutoff = $args{-cutoff} || 1.00e-10;
2198 :     my $maxhits = $args{-maxhits} || 99999999;
2199 :    
2200 :     my $use_pegs = $args{-use_pegs} || '';
2201 : overbeek 1.1
2202 : overbeek 1.26 print STDERR "Processing rep_seq $query_seq\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2203 : overbeek 1.5
2204 : overbeek 1.1 my $query_file = "$FIG_Config::temp/tmp_query.$$.fasta";
2205 :     open(TMP, ">$query_file") || confess "Could not write-open $query_file";
2206 :     &FIG::display_id_and_seq($query_id, \$query_seq, \*TMP);
2207 :     close(TMP)
2208 : overbeek 1.26 || confess "Could not close query-file $query_file --- args:\n" . &flatten_dumper(\%args) . "\n";
2209 : overbeek 1.1 (-s $query_file)
2210 : overbeek 1.26 || confess "Could not write query sequence to $query_file --- args:\n" . &flatten_dumper(\%args) . "\n";
2211 : overbeek 1.1
2212 :     my $db = "$FIG_Config::temp/tmp_orfs.$$.fasta";
2213 : gdpusch 1.63 if (!-s $db) {
2214 : overbeek 1.1 open(DB, ">$db") || confess "Could not write-open sequence file $db";
2215 : overbeek 1.2
2216 : overbeek 1.10 my @fids;
2217 :     if ($use_pegs) {
2218 :     @fids = $self->get_fids_for_type('peg');
2219 :     } else {
2220 :     @fids = $self->get_fids_for_type('orf');
2221 :     }
2222 : overbeek 1.15 print STDERR "Writing ", (scalar @fids), " FIDs to $db\n"
2223 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2224 : overbeek 1.10
2225 : gdpusch 1.63 foreach my $fid (@fids) {
2226 : overbeek 1.5 print STDERR "Writing $fid to $db\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2227 :     &FIG::display_id_and_seq($fid, \$self->get_feature_sequence($fid), \*DB);
2228 : overbeek 1.1 }
2229 :     close(DB) || die "Could not close $db";
2230 : overbeek 1.13
2231 : overbeek 1.1 }
2232 : overbeek 1.13 (-s $db) || confess "ORF file $db is empty";
2233 :    
2234 : gdpusch 1.63 if ((!-s "$db.psq") || ((-M "$db.psq") > (-M $db))) {
2235 : overbeek 1.13 &FIG::run("formatdb -i $db -p T");
2236 :     }
2237 :     (-s "$db.psq") || confess "formatdb of $db failed";
2238 : overbeek 1.2
2239 : olson 1.44 my @sims = `$FIG_Config::ext_bin/blastall -i $query_file -d $db -p blastp -m8 -e $cutoff`;
2240 : overbeek 1.1
2241 : overbeek 1.2 if (@sims == 0) {
2242 : overbeek 1.13 print STDERR "No sims found for query $query_id against $db\n" if $ENV{VERBOSE};
2243 : overbeek 1.2 } else {
2244 :     my $num_sims = @sims;
2245 : gdpusch 1.75 print STDERR "Found $num_sims sim".(($num_sims == 1) ? qq() : qq(s))
2246 :     , " against candidate ",
2247 :     , ($use_pegs ? qq(PEGs) : qq(ORFs))
2248 :     , ":\n"
2249 : overbeek 1.26 , join("", @sims)
2250 :     if $ENV{VERBOSE};
2251 : overbeek 1.2 }
2252 : overbeek 1.1 chomp @sims;
2253 :    
2254 :     my %seen;
2255 :     my $sim;
2256 :     my @orfs;
2257 : gdpusch 1.63 foreach $sim (@sims) {
2258 :     if ((@orfs < $maxhits) && ($sim =~ /^\S+\t(\S+)/o) && (! $seen{$1}) && $self->is_feature($1)) {
2259 : overbeek 1.2 my $orf_id = $1;
2260 :     $seen{$orf_id} = 1;
2261 :     push @orfs, NewGenome::ORF->new($self, $orf_id);
2262 : overbeek 1.1 }
2263 :     }
2264 : overbeek 1.10
2265 : gdpusch 1.75 my $num_orfs = (scalar @orfs);
2266 :     print STDERR "Found $num_orfs candidate "
2267 :     , ($use_pegs ? qq(PEG) : qq(ORF)).(($num_orfs == 1) ? qq() : qq(s))
2268 :     , "\n"
2269 :     if $ENV{VERBOSE};
2270 :    
2271 : overbeek 1.1 return @orfs;
2272 :     }
2273 :    
2274 : gdpusch 1.63 sub mark_overlapping_orfs {
2275 : overbeek 1.7 my ($self, $fid) = @_;
2276 :    
2277 :     return 1;
2278 :     }
2279 :    
2280 : overbeek 1.19 sub promote_remaining_orfs {
2281 : gdpusch 1.74 my ($self, $called_by_fh) = @_;
2282 : overbeek 1.19
2283 : overbeek 1.25 my @fids = $self->get_fids_for_type(qw(rna peg));
2284 : overbeek 1.19
2285 :     my @existing_features = ();
2286 :     foreach my $fid (@fids) {
2287 :     my $loc = $self->get_seed_loc($fid);
2288 :     confess "No loc for FID $fid:\n", &Dumper($self) if (not defined($loc));
2289 :    
2290 :     my ($contig, $beg, $end) = $fig->boundaries_of($loc);
2291 :     my $aliases = $self->get_feature_aliases($fid);
2292 :     my $entry = [$fid, $loc, $aliases,
2293 :     $contig, $beg, $end,
2294 :     &FIG::min($beg,$end),
2295 :     &FIG::max($beg,$end),
2296 :     ($beg < $end) ? "+" : "-",
2297 :     'peg'
2298 :     ];
2299 :     push @existing_features, $entry;
2300 :     }
2301 : gdpusch 1.74
2302 : overbeek 1.19 my @orf_ids = sort { $self->get_feature_length($b) <=> $self->get_feature_length($a) }
2303 :     $self->get_fids_for_type('orf');
2304 :    
2305 :     my @new_orfs = ();
2306 :     foreach my $orf_id (@orf_ids) {
2307 :     my $loc = $self->get_seed_loc($orf_id);
2308 :     confess "No loc for FID $orf_id:\n", &Dumper($self) if (not defined($loc));
2309 :    
2310 :     my ($contig, $beg, $end) = $fig->boundaries_of($loc);
2311 :     my $aliases = $self->get_feature_aliases($orf_id);
2312 :     my $entry = [$orf_id, $loc, $aliases,
2313 :     $contig, $beg, $end,
2314 :     &FIG::min($beg,$end),
2315 :     &FIG::max($beg,$end),
2316 :     ($beg < $end) ? "+" : "-",
2317 :     'peg'
2318 :     ];
2319 :    
2320 :     push @new_orfs, $entry;
2321 :     }
2322 : overbeek 1.25
2323 :    
2324 : overbeek 1.19 my $parms = { min_peg_ln => 90,
2325 :     max_RNA_overlap => 20,
2326 :     max_convergent_overlap => 50,
2327 :     max_divergent_overlap => 150,
2328 :     max_same_strand_overlap => 120
2329 :     };
2330 :    
2331 : gdpusch 1.74 return &update_features($self, \@existing_features, \@new_orfs, $parms,
2332 :     qq(promote_remaining_orfs (no sims)), $called_by_fh
2333 :     );
2334 : overbeek 1.19 }
2335 :    
2336 :     sub update_features {
2337 : gdpusch 1.74 my ($self, $tbl, $new_entries, $parms, $calling_method, $called_by_fh) = @_;
2338 : overbeek 1.25
2339 :     print STDERR "\nEntered update_features\n" if $ENV{VERBOSE};
2340 : overbeek 1.19
2341 :     # constants for positions in "keep" lists
2342 :     use constant UPD_FID => 0;
2343 :     use constant UPD_LOC => 1;
2344 :     use constant UPD_ALIASES => 2;
2345 :     use constant UPD_CONTIG => 3;
2346 :     use constant UPD_START => 4;
2347 :     use constant UPD_STOP => 5;
2348 :     use constant UPD_LEFT => 6;
2349 :     use constant UPD_RIGHT => 7;
2350 :     use constant UPD_STRAND => 8;
2351 :     use constant UPD_TYPE => 9;
2352 : gdpusch 1.74
2353 :     if (not $calling_method) {
2354 :     $calling_method = "unspecified means";
2355 :     }
2356 : overbeek 1.19
2357 :     my $keep = {};
2358 : gdpusch 1.63 foreach my $entry (@$tbl) {
2359 : overbeek 1.25 my ($fid, $loc, $aliases) = @$entry;
2360 : overbeek 1.31 ($fid =~ /^fig\|\d+\.\d+\.(\w+)\.\d+$/o) || confess "Could not parse FID=$fid";
2361 : overbeek 1.19 my $type = $1;
2362 :    
2363 : overbeek 1.31 ($loc =~ m/^(\S+)_(\d+)_(\d+)$/o) || confess "Could not parse loc=$loc";
2364 : overbeek 1.25 my ($contig, $beg, $end) = ($1, $2, $3);
2365 : overbeek 1.19
2366 :     push(@{$keep->{$contig}}, [$fid,
2367 :     $loc,
2368 :     $aliases,
2369 :     $contig,
2370 :     $beg,
2371 :     $end,
2372 :     &FIG::min($beg,$end),
2373 :     &FIG::max($beg,$end),
2374 :     ($beg < $end) ? "+" : "-",
2375 :     $type,
2376 :     ]);
2377 :     }
2378 :    
2379 : gdpusch 1.63 foreach my $contig (keys %$keep) {
2380 : overbeek 1.19 my $x = $keep->{$contig};
2381 :     $keep->{$contig} = [sort { ($a->[UPD_LEFT] <=> $b->[UPD_LEFT]) ||
2382 :     ($b->[UPD_RIGHT] <=> $a->[UPD_RIGHT]) } @$x];
2383 :     }
2384 :    
2385 : overbeek 1.25 if ($ENV{VERBOSE}) {
2386 :     print STDERR "\nkeep:\n";
2387 :     foreach my $contig (sort keys %$keep) {
2388 :     print STDERR ">$contig\n";
2389 :     my $x = $keep->{$contig};
2390 :     for (my $i=0; $i < @$x; ++$i) {
2391 : overbeek 1.26 my $y = &flatten_dumper($x->[$i]);
2392 : overbeek 1.25 print STDERR "keep $i:\t$y\n";
2393 :     }
2394 :     print STDERR "\n";
2395 :     }
2396 :     }
2397 :    
2398 : overbeek 1.19 my $entry;
2399 : overbeek 1.21 print STDERR "\n" if $ENV{VERBOSE};
2400 : gdpusch 1.63 foreach my $new_entry (@$new_entries) {
2401 : overbeek 1.19 my $orf = $new_entry->[UPD_FID];
2402 :     my $loc = $new_entry->[UPD_LOC];
2403 : overbeek 1.25 my ($contig, $beg, $end) = $fig->boundaries_of($loc);
2404 :    
2405 :     my $left = $new_entry->[UPD_LEFT];
2406 :     my $right = $new_entry->[UPD_RIGHT];
2407 :     my $strand = $new_entry->[UPD_STRAND];
2408 : overbeek 1.19
2409 :     #...Unconditionally delete ORF object, regardless of whether it gets promoted...
2410 :     unless ($self->delete_feature($orf)) {
2411 :     confess "Could not delete ORF $orf:\n"
2412 :     , Dumper($self->get_feature_object($orf));
2413 :     }
2414 :    
2415 :     my $peg;
2416 :     my $where;
2417 : gdpusch 1.63 if (defined($where = &keep_this_one($parms, $keep, $contig, $beg, $end, $orf))) {
2418 : overbeek 1.25 print STDERR "Inserting $contig\_$beg\_$end before keep[$where]\n" if $ENV{VERBOSE};
2419 : overbeek 1.21
2420 : gdpusch 1.74 if ($peg = $self->add_feature( -type => 'peg',
2421 :     -loc => $self->from_seed_loc($loc),
2422 :     -annot => [qq(RAST), qq(Called by $calling_method.)],
2423 :     )
2424 :     )
2425 :     {
2426 :     if (defined($called_by_fh)) {
2427 :     print $called_by_fh "$peg\t$calling_method\n";
2428 :     }
2429 :     }
2430 :     else {
2431 : overbeek 1.19 confess "Could not promote $orf to a PEG:\n"
2432 : overbeek 1.21 , Dumper($self);
2433 : overbeek 1.19 }
2434 :    
2435 :     $new_entry->[UPD_FID] = $peg;
2436 : overbeek 1.25 splice(@ { $keep->{$contig} }, $where, 0, $new_entry);
2437 :     }
2438 :     else {
2439 :     print STDERR "Not inserting feature $left$strand$right\n" if $ENV{VERBOSE};
2440 : overbeek 1.19 }
2441 : overbeek 1.21 print STDERR "\n" if $ENV{VERBOSE};
2442 : overbeek 1.22
2443 : overbeek 1.25 if (my $x = $keep->{$contig}) {
2444 :     for (my $i=0; $i < @$x; ++$i) {
2445 : overbeek 1.26 my $y = &flatten_dumper($x->[$i]);
2446 : overbeek 1.25 print STDERR "keep $i:\t$y\n";
2447 :     }
2448 : overbeek 1.22 }
2449 : overbeek 1.25 print STDERR "\n----------\n\n" if $ENV{VERBOSE};
2450 : overbeek 1.19 }
2451 :     # print STDERR Dumper($keep);
2452 :    
2453 : overbeek 1.42 return (scalar keys %$keep);
2454 : overbeek 1.19 }
2455 :    
2456 : gdpusch 1.63 sub keep_this_one {
2457 : overbeek 1.19 my ($parms,$keep,$contig,$beg,$end,$fid) = @_;
2458 :    
2459 :     my ($ln, $x, @overlaps, $left, $right, $strand, $i);
2460 : overbeek 1.25 my $where = undef;
2461 : overbeek 1.19
2462 :     my $min_peg_ln = $parms->{min_peg_ln};
2463 :    
2464 : gdpusch 1.63 if (($ln = (abs($end-$beg)+1)) < $min_peg_ln) {
2465 : overbeek 1.19 print STDERR "FID failed length test: FID=$fid, beg=$beg, end=$end, len=$ln, min=$min_peg_ln\n"
2466 :     , &Dumper($fid,$keep) if $ENV{VERBOSE};
2467 :     return undef;
2468 :     }
2469 :    
2470 : overbeek 1.25 $left = &FIG::min($beg,$end);
2471 :     $right = &FIG::max($beg,$end);
2472 :     $strand = ($beg < $end) ? "+" : "-";
2473 :     print STDERR "Processing $fid: $left$strand$right\n" if $ENV{VERBOSE};
2474 :    
2475 : gdpusch 1.63 if ($x = $keep->{$contig}) {
2476 : overbeek 1.19 for ($i=0; ($i < @$x) && ($left > $x->[$i]->[UPD_RIGHT]); ++$i) {}
2477 : overbeek 1.25 print STDERR "Possible insertion before keep[$i] --- checking for overlaps\n" if $ENV{VERBOSE};
2478 : overbeek 1.19
2479 :     @overlaps = ();
2480 : gdpusch 1.63 while (($i < @$x) && ($right >= $x->[$i]->[UPD_LEFT])) {
2481 : overbeek 1.26 print STDERR "keep[$i]: " , &flatten_dumper($x->[$i]), "\n" if $ENV{VERBOSE};
2482 : overbeek 1.25
2483 :     if ($left <= $x->[$i]->[UPD_LEFT]) {
2484 :     print STDERR "Setting insertion point to before keep[$i]\n" if $ENV{VERBOSE};
2485 :     $where = $i;
2486 :     }
2487 :    
2488 :     if ($ENV{VERBOSE}) {
2489 : overbeek 1.26 my $y = &flatten_dumper($x->[$i]);
2490 : overbeek 1.25 print STDERR " overlap = "
2491 :     , &overlap($left, $right, $x->[$i]->[UPD_LEFT], $x->[$i]->[UPD_RIGHT])
2492 :     , ",\n pushing feature\t$y\n";
2493 :     }
2494 :    
2495 : overbeek 1.19 push(@overlaps,$x->[$i]);
2496 :     ++$i;
2497 :     }
2498 :    
2499 : overbeek 1.25 if (not defined $where) {
2500 :     $where = $i;
2501 :     print STDERR "Insertion point defaults to before keep[$where]\n" if $ENV{VERBOSE};
2502 :     }
2503 :    
2504 : overbeek 1.19 my $serious = 0;
2505 : gdpusch 1.63 for ($i=0; ($i < @overlaps); ++$i) {
2506 : overbeek 1.26 my $y = &flatten_dumper($overlaps[$i]);
2507 : overbeek 1.25 print STDERR " overlap list[$i]:\t$y\n" if $ENV{VERBOSE};
2508 : overbeek 1.19 if (&serious_overlap($parms,$left,$right,$strand,$overlaps[$i],$fid)) { ++$serious; }
2509 : overbeek 1.21 # print STDERR "\n" if $ENV{VERBOSE};
2510 : overbeek 1.19 }
2511 : overbeek 1.21 # print STDERR "\n" if $ENV{VERBOSE};
2512 : overbeek 1.19 return undef if $serious;
2513 :     }
2514 :    
2515 :     return $where;
2516 :     }
2517 :    
2518 : gdpusch 1.63 sub overlap {
2519 : overbeek 1.19 my ($min, $max, $minO, $maxO) = @_;
2520 :     my $olap = &FIG::max(0, (&FIG::min($max,$maxO) - &FIG::max($min,$minO) + 1));
2521 :     # warn "min=$min,\tmax=$max,\tminO=$minO,\tmaxO=$maxO,\tolap=$olap\n" if $ENV{VERBOSE};
2522 :     return $olap;
2523 :     }
2524 :    
2525 :     sub serious_overlap {
2526 :     my ($parms, $min, $max, $strand, $overlap, $fid) = @_;
2527 :     my $minO = $overlap->[UPD_LEFT];
2528 :     my $maxO = $overlap->[UPD_RIGHT];
2529 :     my $strandO = $overlap->[UPD_STRAND];
2530 :     my $typeO = $overlap->[UPD_TYPE];
2531 :     my $fidO = $overlap->[UPD_FID];
2532 :    
2533 :     my $olap = &overlap($min, $max, $minO, $maxO);
2534 :     print STDERR "olap=$olap\n" if $ENV{VERBOSE};
2535 :    
2536 :     if (&embedded($min,$max,$minO,$maxO))
2537 :     {
2538 : overbeek 1.25 print STDERR "$minO$strandO$maxO is embedded in $min$strand$max [kept]\n" if $ENV{VERBOSE};
2539 : overbeek 1.19 # die Dumper($overlap);
2540 :     return 1;
2541 :     }
2542 :    
2543 : gdpusch 1.63 if (&embedded($minO,$maxO,$min,$max)) {
2544 : overbeek 1.19 print STDERR "$min$strand$max is embedded in $minO$strandO$maxO [kept]\n" if $ENV{VERBOSE};
2545 :     # die Dumper($overlap);
2546 :     return 1;
2547 :     }
2548 :    
2549 : gdpusch 1.63 if (($typeO eq "rna") && ($olap > $parms->{max_RNA_overlap})) {
2550 : overbeek 1.19 print STDERR "too much RNA overlap: $min$strand$max overlaps $minO$strandO$maxO ($olap)\n" if $ENV{VERBOSE};
2551 :     return 1;
2552 :     }
2553 :    
2554 : overbeek 1.31 if (($typeO !~ /^rna/o) && ($olap > $parms->{max_convergent_overlap})
2555 : overbeek 1.19 && &convergent($min,$max,$strand,$minO,$maxO,$strandO))
2556 :     {
2557 :     print STDERR "too much convergent overlap: $min$strand$max overlaps $minO$strandO$maxO ($olap)\n" if $ENV{VERBOSE};
2558 :     return 1;
2559 :     }
2560 :    
2561 : overbeek 1.31 if (($typeO !~ /^rna/o) && ($olap > $parms->{max_divergent_overlap})
2562 : overbeek 1.19 && &divergent($min,$max,$strand,$minO,$maxO,$strandO))
2563 :     {
2564 :     print STDERR "too much divergent overlap: $min$strand$max overlaps $minO$strandO$maxO ($olap)\n" if $ENV{VERBOSE};
2565 :     return 1;
2566 :     }
2567 :    
2568 : overbeek 1.31 if (($typeO !~ /^rna/o) && ($strand eq $strandO)
2569 : overbeek 1.19 && ($olap > $parms->{max_same_strand_overlap}))
2570 :     {
2571 :     print STDERR "too much same_strand overlap: $min$strand$max overlaps $minO$strandO$maxO ($olap)\n" if $ENV{VERBOSE};
2572 :     return 1;
2573 :     }
2574 :    
2575 :     return 0;
2576 :     }
2577 :    
2578 :     sub embedded {
2579 :     my($min,$max, $minO,$maxO) = @_;
2580 :    
2581 : gdpusch 1.63 if (($min <= $minO) && ($maxO <= $max)) {
2582 : overbeek 1.19 return 1;
2583 :     }
2584 :     return 0;
2585 :     }
2586 :    
2587 :     sub convergent {
2588 :     my($min,$max,$strand,$minO,$maxO,$strandO) = @_;
2589 :    
2590 :     if (($strand ne $strandO) &&
2591 :     ((($min < $minO) && ($strand eq "+")) ||
2592 :     (($minO < $min) && ($strandO eq "+"))))
2593 :     {
2594 :     return 1;
2595 :     }
2596 :     return 0;
2597 :     }
2598 :    
2599 :     sub divergent {
2600 :     my($min,$max,$strand,$minO,$maxO,$strandO) = @_;
2601 :    
2602 :     if (($strand ne $strandO) &&
2603 :     ((($min < $minO) && ($strand eq "-")) ||
2604 :     (($minO < $min) && ($strandO eq "-"))))
2605 :     {
2606 :     return 1;
2607 :     }
2608 :     return 0;
2609 :     }
2610 :    
2611 :    
2612 :    
2613 : gdpusch 1.63 sub call_start {
2614 : overbeek 1.1 my ($self, $fid, $sims) = @_;
2615 : overbeek 1.26 my ($id2, $ln2, $b1, $b2, $iden, $bsc);
2616 :     my ($tmp_loc, $tmp_len, $loc_str);
2617 : overbeek 1.34 my (%votes_for, %codon_for, $codon);
2618 : overbeek 1.26
2619 :     my $num_sims = @$sims;
2620 :     if (defined($ENV{VERBOSE})) {
2621 :     my $s = ($num_sims == 1) ? qq() : qq(s);
2622 :     print STDERR "\nRe-calling START of $fid using $num_sims sim$s";
2623 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2)) {
2624 : gdpusch 1.59 print STDERR ":\n";
2625 : overbeek 1.26 foreach my $sim (@$sims) {
2626 :     print STDERR (&flatten_dumper($sim), "\n");
2627 :     }
2628 :     }
2629 :     print STDERR "\n";
2630 :     }
2631 :    
2632 :     my $original_loc = $self->get_feature_loc($fid);
2633 : gdpusch 1.63 if (@$original_loc > 1) {
2634 :     warn "Cannot recall START of multi-exon feature $fid";
2635 :     return $original_loc;
2636 :     }
2637 : overbeek 1.26 my $original_len = $self->get_exon_length($original_loc->[-1]);
2638 :     my $orf_len = $self->get_orf_length($original_loc);
2639 :     my $original_start_codon = substr($self->get_dna_subseq($original_loc), 0, 3);
2640 :     if (defined $original_start_codon) {
2641 : overbeek 1.31 if ($original_start_codon =~ m/^[agt]tg$/io) {
2642 : overbeek 1.26 print STDERR "FID $fid has original START=$original_start_codon,"
2643 :     , " original_len=$original_len,"
2644 :     , " orf_len=$orf_len,"
2645 :     , " original_loc = ", &flatten_dumper($original_loc), "\n"
2646 :     if $ENV{VERBOSE};
2647 :     }
2648 :     else {
2649 :     print STDERR "Original START codon $original_start_codon is invalid --- setting undef\n"
2650 :     if $ENV{VERBOSE};
2651 :     $original_start_codon = undef;
2652 :     }
2653 :     }
2654 :     else {
2655 :     confess "Could not get called START codon for FID $fid, loc="
2656 :     , &flatten_dumper($original_loc);
2657 :     }
2658 :    
2659 : gdpusch 1.63 if (not defined($orf_len)) {
2660 : overbeek 1.13 warn "FID $fid has an undefined ORF-length; recomputing now...\n" if $ENV{VERBOSE};
2661 : overbeek 1.26 $original_loc = $self->search_for_upstream_stop($original_loc);
2662 :     $original_loc = $self->set_feature_loc($fid, $original_loc)
2663 : overbeek 1.13 || confess "Could not set ORF-length for $fid";
2664 : overbeek 1.26 ($orf_len = $self->get_orf_length($original_loc))
2665 :     || confess "Setting recomputed ORF-length failed;"
2666 :     , " original_loc=", &flatten_dumper($original_loc), "\n"
2667 :     , " feature=", &flatten_dumper($self->get_feature($fid));
2668 :     }
2669 :     my $orf_loc = $self->copy_loc($original_loc->[-1]);
2670 :    
2671 :    
2672 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2673 :     # ... Compute mean and variance of hit lengths ...
2674 :     #-----------------------------------------------------------------------
2675 : overbeek 1.30 my $mean_len = $original_len;
2676 : overbeek 1.29
2677 :     #...Cheat to suppress divide-by-zero errors
2678 :     my $variance = $original_len / ($num_sims+1);
2679 :     my $std_dev = sqrt($variance);
2680 :    
2681 : overbeek 1.26 my $begin_sim_region = 0;
2682 : gdpusch 1.63 foreach my $sim (@$sims) {
2683 : overbeek 1.29 print STDERR &flatten_dumper($sim), "\n"
2684 : gdpusch 1.60 if (defined($ENV{VERBOSE}) && (($num_sims <= 10) || ($ENV{VERBOSE} > 2)));
2685 : overbeek 1.26
2686 :     $b1 = $sim->b1;
2687 :     $ln2 = 3 * ( 1 + $sim->ln2 );
2688 :    
2689 :     $mean_len += $ln2;
2690 :     $variance += ($ln2 - $original_len)**2;
2691 :    
2692 :     if ($ENV{VERBOSE}) {
2693 :     if ($ENV{VERBOSE} > 2) {
2694 :     print STDERR " b1=$b1,\tbegin=$begin_sim_region --> ";
2695 :     $begin_sim_region = &FIG::max( $begin_sim_region, ($original_len - 3*($b1-1)) );
2696 :     print STDERR "begin=$begin_sim_region\n";
2697 :     }
2698 :     else {
2699 :     $tmp_len = ($original_len - 3*($b1-1));
2700 :     if ($tmp_len > $begin_sim_region) {
2701 :     print STDERR " b1=$b1,\tbegin=$begin_sim_region --> ";
2702 :     $begin_sim_region = $tmp_len;
2703 :     print STDERR "begin=$begin_sim_region\n";
2704 :     }
2705 :     }
2706 :     }
2707 : overbeek 1.23 }
2708 :    
2709 : overbeek 1.29 if ($num_sims > 1) {
2710 : overbeek 1.30 $mean_len /= ($num_sims + 1); #...Add one to acct for including the original length...
2711 :     $variance -= $num_sims * ($mean_len - $original_len)**2; #...Correct for offset...
2712 :     $variance /= $num_sims; #...Again, add one to "sample-variance" denominator...
2713 : overbeek 1.23 $std_dev = sqrt($variance);
2714 : overbeek 1.26 print STDERR "Setting num_sims=$num_sims, begin=$begin_sim_region, mean_len=$mean_len, std_dev=$std_dev\n"
2715 :     if $ENV{VERBOSE};
2716 : overbeek 1.23 }
2717 :     else {
2718 : overbeek 1.26 $mean_len = $original_len;
2719 :     $std_dev = sqrt($original_len);
2720 :     print STDERR "Default num_sims=$num_sims, begin=$begin_sim_region, mean_len=$mean_len, std_dev=$std_dev\n"
2721 :     if $ENV{VERBOSE};
2722 : overbeek 1.23 }
2723 : overbeek 1.12
2724 : overbeek 1.26 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2725 :     # ... Vote on possible START positions ...
2726 :     #-----------------------------------------------------------------------
2727 : overbeek 1.34
2728 :     my @tmp = $self->possible_starts($orf_loc);
2729 :     my @possible_starts = map { $_->[0] } @tmp;
2730 :     %codon_for = map { $_->[0] => $_->[1] } @tmp;
2731 :     my $first_start_len = $possible_starts[0];
2732 :     my $first_start_codon = $codon_for{$first_start_len};
2733 :     confess "Could not extract first START candidate for FID $fid, loc="
2734 :     , &flatten_dumper($orf_loc)
2735 :     unless ($first_start_len && $first_start_codon);
2736 :    
2737 : overbeek 1.26 print STDERR "Possible STARTs: "
2738 : overbeek 1.34 , join(", ", @possible_starts)
2739 :     , "\n"
2740 : overbeek 1.26 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2741 :    
2742 :     #...cast a nominal vote for the original call...
2743 :     $votes_for{$original_len} += &_score_len(10, 80, $original_start_codon, $original_len, $mean_len, $std_dev);
2744 :    
2745 : overbeek 1.34 my $estimate;
2746 : overbeek 1.36 for (my $i=0; $i < @$sims; ++$i) {
2747 :     last if ($i >= 10);
2748 :    
2749 :     my $sim = $sims->[$i];
2750 :    
2751 : overbeek 1.26 $id2 = $sim->id2;
2752 :     $ln2 = 3 * ( 1 + $sim->ln2 );
2753 :    
2754 :     $b1 = $sim->b1;
2755 :     $b2 = $sim->b2;
2756 :     $iden = $sim->iden;
2757 :     $bsc = $sim->bsc;
2758 :    
2759 : overbeek 1.41 $estimate = 3 * ( int($orf_len/3) - $b1 + $b2 );
2760 : overbeek 1.26
2761 : gdpusch 1.63 if ($estimate > $orf_len) {
2762 : overbeek 1.29 print STDERR "For sim \@ $id2: ln2=$ln2, b1=$b1, b2=$b2, estimated len=$estimate > orf_len=$orf_len"
2763 : overbeek 1.26 , " --- using len=$first_start_len\n"
2764 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2765 :     $votes_for{$first_start_len} += &_score_len($bsc, $iden, $first_start_codon, $first_start_len, $mean_len, $std_dev);
2766 :     }
2767 :     else {
2768 : overbeek 1.29 print STDERR "For sim \@ $id2: ln2=$ln2, b1=$b1, b2=$b2, estimated len=$estimate\n"
2769 : overbeek 1.26 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2770 :    
2771 : overbeek 1.34 my $best_len = $first_start_len;
2772 :     $codon = $first_start_codon;
2773 :     foreach my $start (@possible_starts) {
2774 :     if (abs($start - $estimate) < abs($best_len - $estimate)) {
2775 :     $best_len = $start;
2776 :     $codon = $codon_for{$best_len};
2777 : overbeek 1.26 }
2778 :     }
2779 :    
2780 : overbeek 1.34 #...penalize best_len candidate if $best_len < $begin_sim_region...
2781 :     my $weight = ($best_len <= $begin_sim_region) ? (0.5) : (1.0) ;
2782 :     $votes_for{$best_len} += $weight * &_score_len($bsc, $iden, $codon, $best_len, $mean_len, $std_dev);
2783 :     $codon_for{$best_len} = $codon;
2784 :    
2785 :     print STDERR " Candidate is codon=$codon, len=$best_len,"
2786 :     , "\n"
2787 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2788 : overbeek 1.12 }
2789 : overbeek 1.26 }
2790 :     print STDERR "Votes: ", &flatten_dumper(\%votes_for), "\n" if $ENV{VERBOSE};
2791 :    
2792 :     my $best_len = $original_len;
2793 :     my $best_vote = $votes_for{$original_len};
2794 :     print STDERR "begin=$begin_sim_region\n" if $ENV{VERBOSE};
2795 : gdpusch 1.63 foreach my $len (sort { $a <=> $b } keys %votes_for) {
2796 : overbeek 1.26 if ($votes_for{$len} > $best_vote) {
2797 :     $best_len = $len;
2798 :     $best_vote = $votes_for{$len};
2799 :     print STDERR "Best --> len=$best_len, codon=$codon_for{$best_len}, vote=$best_vote\n" if $ENV{VERBOSE};
2800 :     }
2801 :     else {
2802 :     print STDERR " len=$len, codon=$codon_for{$len}, vote=$votes_for{$len}"
2803 :     . " --- keeping best_len=$best_len, best_codon=$codon_for{$best_len}, best_vote=$best_vote\n"
2804 :     if $ENV{VERBOSE};
2805 :     }
2806 :     }
2807 :     print STDERR "Selected START codon=$codon_for{$best_len}, new_len=$best_len, original_len=$original_len, orf_len=$orf_len,"
2808 :     , " given mean_len=$mean_len, std_dev=$std_dev\n"
2809 :     if $ENV{VERBOSE};
2810 :    
2811 : overbeek 1.31 if ($codon_for{$best_len} !~ m/^[atg]tg$/io) {
2812 : overbeek 1.26 if ($ENV{VERBOSE}) {
2813 :     print STDERR "Invalid START codon=$codon_for{$best_len} selected, len=$best_len, orf_len=$orf_len\n";
2814 :     print STDERR "Attempting to default to first_start_len=$first_start_len, first_start_codon=$first_start_codon\n";
2815 :     }
2816 :    
2817 :     if (defined($first_start_codon)) {
2818 :     $best_len = $first_start_len;
2819 :     print STDERR "Default succeded\n" if $ENV{VERBOSE};
2820 :     }
2821 :     else {
2822 :     $loc_str = &flatten_dumper($original_loc);
2823 :     confess "Aborting ---something is seriously wrong with this ORF: $loc_str";
2824 : overbeek 1.12 }
2825 :     }
2826 :    
2827 : overbeek 1.26 $tmp_loc = $self->copy_loc($orf_loc);
2828 :     $self->set_exon_length($tmp_loc, $best_len)
2829 : overbeek 1.29 || confess "Could not set recalled ORF length to loc=" . &flatten_dumper($tmp_loc);
2830 : overbeek 1.26
2831 :     return $tmp_loc;
2832 :     }
2833 :    
2834 :     sub _score_len {
2835 :     my ($bsc, $iden, $codon, $len, $mean_len, $std_dev) = @_;
2836 :    
2837 : overbeek 1.36 return 0 if ($iden < 75.0);
2838 :    
2839 :     return ($bsc * (($iden - 75.0) / 100.0)
2840 : overbeek 1.26 * (($codon eq 'atg') ? 5 : 1)
2841 : overbeek 1.36 * exp(-(0.5) * (($len - $mean_len) / $std_dev)**2)
2842 :     / $std_dev
2843 :     / sqrt(2.0 * 3.14159265359)
2844 :     );
2845 : overbeek 1.26 }
2846 :    
2847 :    
2848 : gdpusch 1.63 sub possible_starts {
2849 : overbeek 1.26 my ($self, $loc, $maxlen) = @_;
2850 :     my @starts = ();
2851 :    
2852 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
2853 :     my ($contig_id, $end, $strand, $len, $orf_len) = @ { $loc->[0] };
2854 :     my $contig_len = $self->get_contig_length($contig_id);
2855 :    
2856 : gdpusch 1.63 if (not defined($orf_len)) {
2857 : overbeek 1.26 print STDERR "loc=", &flatten_dumper($loc), " has an undefined ORF-length; recomputing now...\n"
2858 :     if $ENV{VERBOSE};
2859 :     $loc = $self->search_for_upstream_stop($loc)
2860 :     || confess "Could not find upstrem STOP for ", &flatten_dumper($loc);
2861 :     ($orf_len = $self->get_orf_length($loc))
2862 :     || confess "Could not find ORF-length for loc=", &flatten_dumper($loc);
2863 :     }
2864 :    
2865 :     my $sign = ($strand eq qq(+)) ? +1 : -1 ;
2866 :    
2867 :     my $beg = $end - $sign*($orf_len-1);
2868 :    
2869 :     my $seq = $self->get_dna_subseq($loc);
2870 :     print STDERR "Searching for STARTs, contig=$contig_id ($contig_len),"
2871 : overbeek 1.29 . " beg=$beg, end=$end, strand=$strand, len=$len, orf_len=$orf_len\n"
2872 : overbeek 1.26 . (join(" ", unpack(("A3" x ((2+length($seq))/3)), $seq))) . "\n"
2873 : overbeek 1.37 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2874 : overbeek 1.26 # print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2875 :    
2876 :     my $e = $end - $sign*($orf_len-3);
2877 :     my $codon = "";
2878 : gdpusch 1.63 for (my $i=$orf_len; $i > 3; $i -= 3) {
2879 : overbeek 1.26 $e = $end - $sign*($i-3);
2880 :     $codon = $self->get_dna_subseq([$contig_id, $e, $strand, 3]);
2881 :    
2882 :     print STDERR "PS\t--- $i\t$e\t$codon"
2883 : overbeek 1.37 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2884 : overbeek 1.26
2885 : gdpusch 1.63 if ($codon =~ m/atg|gtg|ttg/io) {
2886 : overbeek 1.26 $len = $i;
2887 :     push @starts, [$len, $codon];
2888 : overbeek 1.37 print STDERR " <---" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2889 : overbeek 1.26 }
2890 :    
2891 : overbeek 1.37 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2892 : overbeek 1.26
2893 : gdpusch 1.63 if ((not $codon) || (&FIG::translate($codon, $self->get_translation_table) eq '*')) {
2894 : gdpusch 1.55 confess "Found internal STOP at strand=$strand, len=$i, e=$e, end=$end, codon=$codon\n\n"
2895 : overbeek 1.12 }
2896 :     }
2897 : overbeek 1.37 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2898 : overbeek 1.1
2899 : overbeek 1.26 return @starts;
2900 : overbeek 1.1 }
2901 :    
2902 : overbeek 1.29 my %cache_downstream_start_loc;
2903 :     my %cache_downstream_start_codon;
2904 : gdpusch 1.63 sub search_for_downstream_start {
2905 : overbeek 1.12 my ($self, $loc, $minlen) = @_;
2906 : overbeek 1.1 unless (defined($minlen)) { $minlen = 90; }
2907 : overbeek 1.29 my ($orf_loc, $codon);
2908 : overbeek 1.1
2909 : overbeek 1.12 if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
2910 : overbeek 1.29 my $x = &flatten_dumper($loc);
2911 :     if (defined($orf_loc = $cache_downstream_start_loc{$x})) {
2912 :     $codon = $cache_downstream_start_codon{$x};
2913 : overbeek 1.31 print STDERR " Returning cached downstream START codon=$codon, loc="
2914 : overbeek 1.29 , &flatten_dumper($orf_loc), "\n"
2915 :     if $ENV{VERBOSE};
2916 :     if (wantarray) { return ($orf_loc, $codon); }
2917 :     else { return $orf_loc; }
2918 :     }
2919 :    
2920 : overbeek 1.12 my ($contig_id, $end, $strand, $len, $orf_len) = @ { $loc->[0] };
2921 : overbeek 1.29 confess "Bad loc=", &flatten_dumper($loc) unless ($contig_id && $end && $strand && $len);
2922 : overbeek 1.12
2923 : overbeek 1.26 my $sign = ($strand eq qq(+)) ? +1 : -1 ;
2924 : overbeek 1.29 my $beg = $end - $sign*($len-1);
2925 : overbeek 1.12
2926 : overbeek 1.26 print STDERR "Searching for downstream START, contig=$contig_id,"
2927 :     . " beg=$beg, end=$end, strand=$strand, len=$len, orf_len=$orf_len"
2928 : overbeek 1.12 if $ENV{VERBOSE};
2929 : overbeek 1.26 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2930 : overbeek 1.12
2931 : overbeek 1.29 my $e;
2932 :     $codon = "";
2933 : gdpusch 1.63 for (my $i=$len; $i >= $minlen; $i -= 3) {
2934 : overbeek 1.26 $e = $end - $sign*($i-3);
2935 :     $codon = $self->get_dna_subseq([$contig_id, $e, $strand, 3]);
2936 :     print STDERR "SD^\t--- $i\t$e\t$codon\n"
2937 : overbeek 1.41 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2938 : overbeek 1.12
2939 : gdpusch 1.63 if ($codon =~ m/atg|gtg|ttg/io) {
2940 : overbeek 1.12 $len = $i;
2941 :     last;
2942 :     }
2943 :    
2944 : gdpusch 1.63 if ((not $codon) || (&FIG::translate($codon, $self->get_translation_table) eq '*')) {
2945 : overbeek 1.26 cluck "Search failed ($i, $e, $end, $codon) --- going with default length\n" if $ENV{VERBOSE};
2946 : overbeek 1.12 print STDERR "\n" if $ENV{VERBOSE};
2947 :    
2948 :     $len = $orf_len;
2949 :     last;
2950 :     }
2951 :     }
2952 : overbeek 1.26 print STDERR " --> len=$len, codon=$codon\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} == 1));
2953 : overbeek 1.41 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
2954 : overbeek 1.12
2955 : overbeek 1.29 $orf_loc = [ $contig_id, $end, $strand, $len, $orf_len ];
2956 :     $cache_downstream_start_loc{$x} = $orf_loc;
2957 :     $cache_downstream_start_codon{$x} = $codon;
2958 : gdpusch 1.74
2959 : overbeek 1.29 if (wantarray) { return ($orf_loc, $codon); }
2960 :     else { return $orf_loc; }
2961 : overbeek 1.12 }
2962 :    
2963 : overbeek 1.29 my %cache_upstream_start_loc;
2964 :     my %cache_upstream_start_codon;
2965 : gdpusch 1.63 sub search_for_upstream_start {
2966 : overbeek 1.12 my ($self, $loc, $maxlen) = @_;
2967 : overbeek 1.29 my ($orf_loc, $codon);
2968 : overbeek 1.12
2969 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
2970 : overbeek 1.29 my $x = &flatten_dumper($loc);
2971 :     if (defined($orf_loc = $cache_upstream_start_loc{$x})) {
2972 :     $codon = $cache_upstream_start_codon{$x};
2973 : overbeek 1.31 print STDERR " Returning cached upstream START codon=$codon, loc="
2974 : overbeek 1.29 , &flatten_dumper($orf_loc), "\n"
2975 :     if $ENV{VERBOSE};
2976 :     if (wantarray) { return ($orf_loc, $codon); }
2977 :     else { return $orf_loc; }
2978 :     }
2979 :    
2980 : overbeek 1.26 my ($contig_id, $end, $strand, $len, $orf_len) = @ { $loc->[-1] };
2981 :     if (not defined($orf_len)) {
2982 : overbeek 1.12 confess "Undefined ORF-length for loc:\n"
2983 : overbeek 1.26 , &flatten_dumper($loc);
2984 : overbeek 1.1 }
2985 :    
2986 : overbeek 1.26 my $sign = ($strand eq qq(+)) ? +1 : -1 ;
2987 :     my $beg = $end - $sign*($len-1);
2988 : overbeek 1.29 my $contig_len = $self->get_contig_length($contig_id);
2989 : overbeek 1.12
2990 : overbeek 1.26 print STDERR "Searching for upstream START, contig=$contig_id ($contig_len),"
2991 :     . " beg=$beg, end=$end, strand=$strand, len=$len, orf_len=$orf_len"
2992 : overbeek 1.12 if $ENV{VERBOSE};
2993 : overbeek 1.26 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
2994 : overbeek 1.1
2995 : overbeek 1.29 my $e;
2996 :     $codon = "";
2997 : gdpusch 1.63 for (my $i=$len; $i <= ($orf_len-3); $i += 3) {
2998 : overbeek 1.12 $e = $end - $sign*($i-3);
2999 :     $codon = $self->get_dna_subseq([$contig_id, $e, $strand, 3]);
3000 : overbeek 1.26 print STDERR "SU^\t--- $i\t$e\t$codon\n"
3001 : overbeek 1.41 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3002 : overbeek 1.1
3003 : gdpusch 1.63 if ($codon =~ m/atg|gtg|ttg/io) {
3004 : overbeek 1.1 $len = $i;
3005 : overbeek 1.26 last;
3006 : overbeek 1.1 }
3007 :    
3008 : gdpusch 1.63 if ((not $codon) || (&FIG::translate($codon, $self->get_translation_table) eq '*')) {
3009 : overbeek 1.12 cluck "Search failed ($i, $e, $end, $codon) --- going with default length\n" if $ENV{VERBOSE};
3010 :     print STDERR "\n" if $ENV{VERBOSE};
3011 :    
3012 : overbeek 1.1 $len = $orf_len;
3013 :     last;
3014 :     }
3015 :     }
3016 : overbeek 1.26 print STDERR " --> len=$len, codon=$codon\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} == 1));
3017 : overbeek 1.41 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3018 : overbeek 1.1
3019 : overbeek 1.29 $orf_loc = [ $contig_id, $end, $strand, $len, $orf_len ];
3020 :     $cache_upstream_start_loc{$x} = $orf_loc;
3021 :     $cache_upstream_start_codon{$x} = $codon;
3022 : gdpusch 1.74
3023 : overbeek 1.29 if (wantarray) { return ($orf_loc, $codon); }
3024 :     else { return $orf_loc; }
3025 :     }
3026 : overbeek 1.1
3027 : overbeek 1.29 my %cache_upstream_stop_loc;
3028 :     my %cache_upstream_stop_codon;
3029 : gdpusch 1.63 sub search_for_upstream_stop {
3030 : overbeek 1.12 my ($self, $loc, $maxlen) = @_;
3031 : overbeek 1.29 my ($orf_loc, $codon);
3032 :    
3033 : overbeek 1.26 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3034 :     # ... NOTE: Returns an ORF and a codon, not an exon ptr and a codon ...
3035 :     #-----------------------------------------------------------------------
3036 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
3037 : overbeek 1.29 my $x = &flatten_dumper($loc);
3038 :     if (defined($orf_loc = $cache_upstream_stop_loc{$x})) {
3039 :     $codon = $cache_upstream_stop_codon{$x};
3040 : overbeek 1.31 print STDERR " Returning cached upstream STOP codon=$codon, loc="
3041 : overbeek 1.29 , &flatten_dumper($orf_loc), "\n"
3042 :     if $ENV{VERBOSE};
3043 :     if (wantarray) { return ($orf_loc, $codon); }
3044 :     else { return $orf_loc; }
3045 :     }
3046 :    
3047 : overbeek 1.26 my ($contig_id, $end, $strand, $len, $orf_len) = @ { $loc->[-1] };
3048 :     confess "Bad loc=", &flatten_dumper($loc) unless ($contig_id && $end && $strand && $len);
3049 : overbeek 1.12
3050 : overbeek 1.26 my $sign = ($strand eq qq(+)) ? +1 : -1 ;
3051 : overbeek 1.12 my $contig_len = $self->get_contig_length($contig_id);
3052 :    
3053 : gdpusch 1.63 if (defined($orf_len)) {
3054 : overbeek 1.26 print STDERR "search_for_upstream_stop: ORF-length already defined for loc=\n"
3055 :     , &flatten_dumper($loc)
3056 :     if $ENV{VERBOSE};
3057 :    
3058 :     $codon = $self->get_dna_subseq([$contig_id, ($end - $sign*$orf_len), $strand, 3]);
3059 : overbeek 1.12 }
3060 : gdpusch 1.63 else {
3061 : overbeek 1.12 $orf_len = $len; #...default length...
3062 : overbeek 1.26
3063 : gdpusch 1.63 unless (defined($maxlen)) {
3064 : overbeek 1.26 #...Check this for correctness...
3065 :     if ($strand eq qq(+)) {
3066 :     $maxlen = $end;
3067 :     } else {
3068 :     $maxlen = $contig_len - $end + 1;
3069 :     }
3070 :     }
3071 :    
3072 :     my ($beg, $minlen);
3073 :     if (defined($len)) {
3074 :     $minlen = $len;
3075 :     $beg = $end - $sign*($len-1);
3076 : overbeek 1.12 } else {
3077 : overbeek 1.26 $minlen = 6;
3078 :     $beg = $end - $sign*($maxlen-1);
3079 : overbeek 1.12 }
3080 :    
3081 : overbeek 1.26 print STDERR "Searching for upstream STOP, contig=$contig_id ($contig_len),"
3082 :     , " beg=$beg, end=$end, strand=$strand, len=$len, minlin=$minlen, maxlen=$maxlen\n"
3083 :     if $ENV{VERBOSE};
3084 :    
3085 : overbeek 1.29 $codon = "";
3086 : overbeek 1.26 my ($e, $failed, $succeeded);
3087 : gdpusch 1.63 for (my $i=$minlen; $i <= $maxlen; $i += 3) {
3088 : overbeek 1.26 $e = $end - $sign*($i-3);
3089 :     $codon = $self->get_dna_subseq([$contig_id, $e, $strand, 3]);
3090 :     print STDERR "SU*\t--- $i\t$e\t$codon\n"
3091 : overbeek 1.41 if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3092 : overbeek 1.26
3093 : gdpusch 1.63 if ($codon) {
3094 : overbeek 1.26 $orf_len = $i;
3095 : gdpusch 1.63 if (&FIG::translate($codon, $self->get_translation_table) eq '*') {
3096 : overbeek 1.26 $succeeded = 1;
3097 :     $orf_len = ($i-3);
3098 :     last;
3099 :     }
3100 :     }
3101 : gdpusch 1.63 else {
3102 : overbeek 1.26 $failed = 1;
3103 :     print STDERR "Search failed ($i, $e, $end) --- going with last successful length, $orf_len\n";
3104 :     last;
3105 : overbeek 1.12 }
3106 :     }
3107 : overbeek 1.41 print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3108 : overbeek 1.26
3109 :     $codon = $self->get_dna_subseq([$contig_id, ($end - $sign*$orf_len), $strand, 3]);
3110 : gdpusch 1.63 if ($codon
3111 :     && (&FIG::translate($codon, $self->get_translation_table) ne '*')
3112 :     && ($succeeded || (not $failed))
3113 :     )
3114 : overbeek 1.12 {
3115 : overbeek 1.26 confess "Something is wrong --- non-failed search yielded preceeding codon=$codon";
3116 : overbeek 1.12 }
3117 :     }
3118 :    
3119 : overbeek 1.29 $orf_loc = [ $contig_id, $end, $strand, $len, $orf_len ];
3120 :     $cache_upstream_stop_loc{$x} = $orf_loc;
3121 :     $cache_upstream_stop_codon{$x} = $codon;
3122 : gdpusch 1.74
3123 : overbeek 1.29 if (wantarray) { return ($orf_loc, $codon); }
3124 :     else { return $orf_loc; }
3125 : overbeek 1.12 }
3126 : gdpusch 1.52
3127 : gdpusch 1.74 my %cache_extend_to_downstream_stop_loc;
3128 :     my %cache_extend_to_downstream_stop_codon;
3129 :     sub extend_to_downstream_stop {
3130 :     my ($self, $loc, $minlen) = @_;
3131 :     unless (defined($minlen)) { $minlen = 90; }
3132 :     my ($orf_loc, $codon);
3133 :    
3134 :     if (ref($loc->[0]) ne 'ARRAY') { $loc = [$loc]; }
3135 :     my $x = &flatten_dumper($loc);
3136 :     if (defined($orf_loc = $cache_extend_to_downstream_stop_loc{$x})) {
3137 :     $codon = $cache_extend_to_downstream_stop_codon{$x};
3138 :     print STDERR " Returning cached downstream START codon=$codon, loc="
3139 :     , &flatten_dumper($orf_loc), "\n"
3140 :     if $ENV{VERBOSE};
3141 :     if (wantarray) { return ($orf_loc, $codon); }
3142 :     else { return $orf_loc; }
3143 :     }
3144 :    
3145 :     my ($contig_id, $end, $strand, $len, $orf_len) = @ { $loc->[0] };
3146 :     confess "Bad loc=", &flatten_dumper($loc) unless ($contig_id && $end && $strand && $len);
3147 :    
3148 :     my $sign = ($strand eq qq(+)) ? +1 : -1 ;
3149 :     my $beg = $end - $sign*($len-1);
3150 :    
3151 :     print STDERR "Searching for downstream STOP, contig=$contig_id,"
3152 :     . " beg=$beg, end=$end, strand=$strand, len=$len, orf_len=$orf_len"
3153 :     if $ENV{VERBOSE};
3154 :     print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
3155 :    
3156 :     my $e;
3157 :     $codon = "";
3158 :     for (my $i=0; $i < 9000; $i += 3) {
3159 :     $e = $end + $sign * $i;
3160 :     last unless $self->check_bounds([$contig_id, $e, $strand, 3]);
3161 :    
3162 :     $codon = $self->get_dna_subseq([$contig_id, $e, $strand, 3]);
3163 :     print STDERR "ED\*\t--- $i\t$e\t$codon\n"
3164 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3165 :    
3166 :     if (&FIG::translate($codon, $self->get_translation_table) eq '*') {
3167 :     last;
3168 :     }
3169 :    
3170 :     $end = $e;
3171 :     $len += 3;
3172 :     if (defined($orf_len)) { $orf_len += 3; }
3173 :     }
3174 :     print STDERR " --> len=$len, codon=$codon\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} == 1));
3175 :     print STDERR "\n" if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 2));
3176 :    
3177 :     $orf_loc = [ $contig_id, $end, $strand, $len, $orf_len ];
3178 :     $cache_extend_to_downstream_stop_loc{$x} = $orf_loc;
3179 :     $cache_extend_to_downstream_stop_codon{$x} = $codon;
3180 :    
3181 :     if (wantarray) { return ($orf_loc, $codon); }
3182 :     else { return $orf_loc; }
3183 :     }
3184 : gdpusch 1.52
3185 :    
3186 :     sub kmer_counts {
3187 :     my ($self, $kmer_lens, $seqs ) = @_ ;
3188 :    
3189 :     $kmer_lens = [ sort {$b <=> $a} eval($kmer_lens) ];
3190 :     my $max_kmer_len = $kmer_lens->[0];
3191 :    
3192 :     my $table = {};
3193 :     foreach my $kmer_len ( @$kmer_lens ) {
3194 :     my $counts = $table->{$kmer_len} = {};
3195 :    
3196 :     #...initialize subtable 'totals' counters: [ grand, frame1, frame2, frame3 ]
3197 :     $counts->{0} = [ 0, 0, 0, 0 ];
3198 :    
3199 :     foreach my $kmer ( &_generate_kmers($kmer_len) ) {
3200 :     $counts->{$kmer} = [0,0,0,0];
3201 :     }
3202 :     }
3203 :    
3204 :     foreach my $seq ( @$seqs ) {
3205 : gdpusch 1.53 $seq = lc( $seq );
3206 : gdpusch 1.52 for ( my $i = 0, my $frame = 1;
3207 :     ( $i < length($seq) );
3208 :     ++$i, $frame = (($frame % 3)+1) )
3209 :     {
3210 :     foreach my $kmer_len ( @$kmer_lens ) {
3211 : gdpusch 1.53 my $kmer = substr( $seq, $i, $kmer_len );
3212 : gdpusch 1.52 next if ($kmer_len > length($kmer)); #...skip too-short entries at end of seq
3213 :    
3214 :     my $counts = $table->{$kmer_len};
3215 :     my $totals = $counts->{"0"};
3216 :    
3217 :     ++ $totals->[ 0 ];
3218 :     ++ $totals->[ $frame ];
3219 :    
3220 :     my @kmers = &_expand_ambigs( $kmer );
3221 : gdpusch 1.55 # print STDERR "Expanded $kmer to ", (scalar @kmers), " kmers\n" if ($ENV{VERBOSE} && (@kmers > 1));
3222 : gdpusch 1.52 my $frac = 1 / (scalar @kmers);
3223 :    
3224 :     foreach my $kmer ( @kmers ) {
3225 :     my $count = $counts->{$kmer};
3226 :     $count->[ 0 ] += $frac;
3227 :     $count->[ $frame ] += $frac;
3228 :     }
3229 :     }
3230 :     }
3231 :     }
3232 : gdpusch 1.55 # print STDERR "\n\n" if $ENV{VERBOSE};
3233 : gdpusch 1.52
3234 :     return $table;
3235 :     }
3236 :    
3237 :     sub _generate_kmers {
3238 :     my ($k) = @_;
3239 :    
3240 :     if ($k == 0) {
3241 :     return ("");
3242 :     }
3243 :     else {
3244 :     my @kmers = ();
3245 :     foreach my $kmer ( &_generate_kmers($k-1) ) {
3246 : gdpusch 1.53 foreach my $base (qw(a c g t)) {
3247 : gdpusch 1.52 push( @kmers, $base.$kmer);
3248 :     }
3249 :     }
3250 :     return @kmers;
3251 :     }
3252 :     }
3253 :    
3254 :     sub _expand_ambigs {
3255 :     my (@stack) = @_;
3256 : gdpusch 1.55 print STDERR "Expanding ", join(", ", @stack), "\n"
3257 :     if ($ENV{VERBOSE} && (@stack > 1));
3258 : gdpusch 1.52
3259 :     my @out;
3260 :     while (@stack > 0) {
3261 : gdpusch 1.53 # m = (a|c)
3262 :     if ($stack[0] =~ m/^([^m]*)m(.*)$/) {
3263 : gdpusch 1.52 shift( @stack );
3264 : gdpusch 1.53 unshift( @stack, ( "$1a$2", "$1c$2") );
3265 : gdpusch 1.52 }
3266 :    
3267 : gdpusch 1.53 # r = (a|g)
3268 :     if ($stack[0] =~ m/^([^r]*)r(.*)$/) {
3269 : gdpusch 1.52 shift( @stack );
3270 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1g$2") );
3271 : gdpusch 1.52 }
3272 :    
3273 : gdpusch 1.53 # w = (a|t)
3274 :     if ($stack[0] =~ m/^([^w]*)w(.*)$/) {
3275 : gdpusch 1.52 shift( @stack );
3276 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1t$2") );
3277 : gdpusch 1.52 }
3278 :    
3279 : gdpusch 1.53 # s = (c|g)
3280 :     if ($stack[0] =~ m/^([^s]*)s(.*)$/) {
3281 : gdpusch 1.52 shift( @stack );
3282 : gdpusch 1.53 unshift( @stack, ("$1c$2", "$1g$2") );
3283 : gdpusch 1.52 }
3284 :    
3285 : gdpusch 1.53 # y = (c|t)
3286 :     if ($stack[0] =~ m/^([^y]*)y(.*)$/) {
3287 : gdpusch 1.52 shift( @stack );
3288 : gdpusch 1.53 unshift( @stack, ("$1c$2", "$1t$2") );
3289 : gdpusch 1.52 }
3290 :    
3291 : gdpusch 1.53 # k = (g|t)
3292 :     if ($stack[0] =~ m/^([^k]*)k(.*)$/) {
3293 : gdpusch 1.52 shift( @stack );
3294 : gdpusch 1.53 unshift( @stack, ("$1g$2", "$1t$2") );
3295 : gdpusch 1.52 }
3296 :    
3297 : gdpusch 1.53 # v = (a|c|g)
3298 :     if ($stack[0] =~ m/^([^v]*)v(.*)$/) {
3299 : gdpusch 1.52 shift( @stack );
3300 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1c$2", "$1g$2") );
3301 : gdpusch 1.52 }
3302 :    
3303 : gdpusch 1.53 # h = (a|c|t)
3304 :     if ($stack[0] =~ m/^([^h]*)h(.*)$/) {
3305 : gdpusch 1.52 shift( @stack );
3306 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1c$2", "$1t$2") );
3307 : gdpusch 1.52 }
3308 :    
3309 : gdpusch 1.53 # d = (a|g|t)
3310 :     if ($stack[0] =~ m/^([^d]*)d(.*)$/) {
3311 : gdpusch 1.52 shift( @stack );
3312 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1g$2", "$1t$2") );
3313 : gdpusch 1.52 }
3314 :    
3315 : gdpusch 1.53 # b = (c|g|t)
3316 :     if ($stack[0] =~ m/^([^b]*)b(.*)$/) {
3317 : gdpusch 1.52 shift( @stack );
3318 : gdpusch 1.53 unshift( @stack, ("$1c$2", "$1g$2", "$1t$2") );
3319 : gdpusch 1.52 }
3320 :    
3321 : gdpusch 1.53 # (n|x) = (a|c|g|t)
3322 :     if ($stack[0] =~ m/^([^xn]*)[xn](.*)$/) {
3323 : gdpusch 1.52 shift( @stack );
3324 : gdpusch 1.53 unshift( @stack, ("$1a$2", "$1c$2", "$1g$2", "$1t$2") );
3325 : gdpusch 1.52 }
3326 :    
3327 : gdpusch 1.53 while ( (@stack > 0) && ($stack[0] !~ m/[mrwsykvhdbxn]/)) {
3328 : gdpusch 1.52 push( @out, shift(@stack) );
3329 :     }
3330 :    
3331 :     last if (@stack == 0);
3332 :     }
3333 :    
3334 :     return @out;
3335 :     }
3336 : gdpusch 1.55
3337 :     sub print_kmer_table
3338 :     {
3339 :     my ($self, $table) = @_;
3340 :    
3341 :     my @kmer_lens = ( sort {$a <=> $b} keys %$table );
3342 :    
3343 :     my $s = (@kmer_lens > 1) ? qq(s) : qq();
3344 :     print STDOUT "# kmer tables for length$s: ", join(qq(,), @kmer_lens), "\n";
3345 :     foreach my $kmer_len (@kmer_lens) {
3346 :     my $counts = $table->{$kmer_len};
3347 :     foreach my $kmer (sort keys %$counts) {
3348 :     my $entry = $counts->{$kmer};
3349 :     print STDOUT "$kmer\t", join("\t", @$entry), "\n";
3350 :     }
3351 :     print STDOUT "//\n";
3352 :     }
3353 :     }
3354 : overbeek 1.26 1;
3355 : overbeek 1.12
3356 : overbeek 1.2
3357 :    
3358 : overbeek 1.26 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3359 : overbeek 1.2 package NewGenome::ORF;
3360 : overbeek 1.26 #-----------------------------------------------------------------------
3361 : overbeek 1.1
3362 : overbeek 1.46 use ToCall;
3363 : overbeek 1.1 use Data::Dumper;
3364 : overbeek 1.2 use Carp qw(:DEFAULT cluck);
3365 : overbeek 1.1
3366 :     sub new {
3367 : overbeek 1.2 my($class, $newG, $fid) = @_;
3368 : overbeek 1.1
3369 : overbeek 1.2 return bless [$newG, $fid], $class;
3370 : overbeek 1.1 }
3371 :    
3372 :     sub seq {
3373 : overbeek 1.2 my ($self) = @_;
3374 :    
3375 :     # confess Dumper($self);
3376 :     my ($newG, $fid) = @$self;
3377 :    
3378 : overbeek 1.5 my $seq;
3379 :     if (defined($seq = $newG->get_feature_sequence($fid)))
3380 : overbeek 1.2 {
3381 : overbeek 1.5 return $seq;
3382 : overbeek 1.2 } else {
3383 :     cluck "!!! Could not get sequence for $fid";
3384 : overbeek 1.26 return undef;
3385 : overbeek 1.2 }
3386 :     }
3387 : overbeek 1.1
3388 : gdpusch 1.49 sub orf_seq {
3389 :     my ($self) = @_;
3390 :    
3391 :     # confess Dumper($self);
3392 :     my ($newG, $fid) = @$self;
3393 :    
3394 :     if ($newG->get_feature_type($fid) ne 'orf') {
3395 :     cluck "!!! Not an ORF: $fid";
3396 :     return undef;
3397 :     }
3398 :    
3399 :     my $seq;
3400 :     if (defined($seq = $newG->get_orf_sequence($fid)))
3401 :     {
3402 :     return $seq;
3403 :     } else {
3404 :     cluck "!!! Could not get ORF sequence for $fid";
3405 :     return undef;
3406 :     }
3407 :     }
3408 :    
3409 : overbeek 1.2 sub get_fid
3410 :     {
3411 :     my ($self) = @_;
3412 :    
3413 :     return $self->[1];
3414 :     }
3415 :    
3416 :     sub set_fid
3417 :     {
3418 :     my ($self, $fid) = @_;
3419 :    
3420 :     return ($self->[1] = $fid);
3421 : overbeek 1.1 }
3422 :    
3423 : overbeek 1.12 sub call_start
3424 :     {
3425 :     my ($self, $sims) = @_;
3426 :     my ($newG, $orf_id) = @$self;
3427 : overbeek 1.46
3428 : overbeek 1.12 my $loc = $newG->call_start($orf_id, $sims);
3429 :     confess "Could not call START for $orf_id", &Dumper($self, $sims) unless defined($loc);
3430 :    
3431 :     return $loc
3432 :     }
3433 :    
3434 : overbeek 1.1 sub promote_to_peg {
3435 : gdpusch 1.74 my ($self, $sims, $func, $annotation) = @_;
3436 : overbeek 1.2 my ($newG, $orf_id) = @$self;
3437 :     my ($loc, $new_peg);
3438 : overbeek 1.1
3439 : overbeek 1.13 if (defined($sims)) {
3440 :     $loc = $newG->call_start($orf_id, $sims);
3441 :     } else {
3442 :     $loc = $newG->get_feature_loc($orf_id);
3443 :     }
3444 : overbeek 1.1
3445 : overbeek 1.26 print STDERR "Attempting to promote $orf_id to a PEG\n" if $ENV{VERBOSE};
3446 :    
3447 : overbeek 1.10 unless ($newG->delete_feature($orf_id))
3448 :     {
3449 :     confess "Could not delete ORF $orf_id:\n"
3450 :     , Dumper($newG->get_feature_object($orf_id));
3451 :     }
3452 :    
3453 : gdpusch 1.74 unless ($new_peg = $newG->add_feature( -type => 'peg', -loc => $loc, -func => $func, -annot => $annotation ) )
3454 : overbeek 1.10 {
3455 : overbeek 1.2 confess "Could not promote $orf_id to a PEG:\n"
3456 : overbeek 1.10 , Dumper($newG, $new_peg, $self);
3457 : overbeek 1.1 }
3458 :    
3459 : overbeek 1.13 unless ($self->set_fid($new_peg))
3460 : overbeek 1.10 {
3461 :     confess "Could not reset $orf_id to $new_peg";
3462 : overbeek 1.1 }
3463 :    
3464 : overbeek 1.12 print STDERR "Promoted $orf_id to $new_peg\n" if $ENV{VERBOSE};
3465 :    
3466 : overbeek 1.2 return $new_peg;
3467 : overbeek 1.1 }
3468 :    
3469 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3