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

Annotation of /FigKernelPackages/SeedDas.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.20 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : olson 1.1 package SeedDas;
19 : olson 1.3 use Data::Dumper;
20 :     #$Data::Dumper::Deparse = 1;
21 :    
22 :     use strict;
23 :    
24 : olson 1.4 use Text::ParseWords;
25 :    
26 : olson 1.3 use FIG;
27 :     use FIG_Config;
28 : olson 1.12 use SproutSearch;
29 : olson 1.1
30 :     #
31 :     # DAS Support
32 :     #
33 :    
34 :     sub new
35 :     {
36 :     my($class, $fig, $data_dir, $url, $dsn) = @_;
37 :    
38 :     my $genome;
39 :    
40 : olson 1.3 if ($dsn =~ /^(\d+)[_.](\d+)$/)
41 : olson 1.1 {
42 : parrello 1.19 $genome = "$1.$2";
43 : olson 1.1 }
44 :     else
45 :     {
46 : parrello 1.19 warn "Genome not found for dsn $dsn\n";
47 : olson 1.1 }
48 :    
49 :     my $self = {
50 : parrello 1.19 fig => $fig,
51 :     dir => $data_dir,
52 :     url => $url,
53 :     dsn => $dsn,
54 :     genome => $genome,
55 :     peg_index => {},
56 :     cache_peg_name => {},
57 :     cache_peg_func => {},
58 : olson 1.1 };
59 :    
60 : olson 1.3 bless $self, $class;
61 :    
62 :     $self->load_track_definitions();
63 :     $self->load_track_bindings();
64 :    
65 : olson 1.4 # warn Dumper($self);
66 :    
67 : olson 1.3 return $self;
68 :     }
69 :    
70 :    
71 :     =pod
72 :    
73 :     =head1 load_tracks
74 :    
75 :     Load the track definitions (binding between track name and code that
76 :     implements the track), and the track bindings (binding between a
77 :     genome and the tracks that are available on that genome).
78 :    
79 :    
80 :     Track definitions are found in FIG/Data/Global/DAS/Track.defs.
81 :    
82 :     Track bindings are found in FIG/Data/Global/DAS/Default.bindings (bindings
83 :     for all genomes), with additional per-genome bindings in
84 :     DAS/<GenomeID>.bindings.
85 : olson 1.4
86 :     After loading a binding, we perform variable expansion on the fields in the
87 :     track definition that maps to this binding, and store a copy in the binding.
88 : olson 1.3
89 :     =cut
90 :    
91 :     sub load_track_definitions
92 :     {
93 :     my($self) = @_;
94 :    
95 :     my $trackdefs = "$self->{dir}/Track.defs";
96 :     my $tfh;
97 :    
98 :     if (!open($tfh, "<$trackdefs"))
99 :     {
100 : parrello 1.19 warn "Could not open $trackdefs: $!\n";
101 :     return undef;
102 : olson 1.3 }
103 :    
104 :     local $/ = "//\n";
105 :    
106 :     while (<$tfh>)
107 :     {
108 : parrello 1.19 chomp;
109 : olson 1.3
110 : parrello 1.19 next if $_ eq "";
111 : olson 1.3
112 : parrello 1.19 my($name, @defs) = split(/\n/, $_);
113 : olson 1.3
114 : parrello 1.19 my $track = {};
115 : olson 1.3
116 : parrello 1.19 for my $defline (@defs)
117 :     {
118 :     my($k, $v) = split(/\t/, $defline, 2);
119 :     $track->{$k} = $v;
120 :     }
121 :    
122 :     $self->{tracks}->{$name} = $track;
123 :     push(@{$self->{track_list}}, $track);
124 :    
125 :     #
126 :     # map perl binding to a code ref.
127 :     #
128 :    
129 :     if (my $func = $track->{perl_binding})
130 :     {
131 :     my $ref = \&$func;
132 :    
133 :     if (defined(&$ref))
134 :     {
135 :     $track->{perl_binding_ref} = $ref;
136 :     }
137 :     else
138 :     {
139 :     warn "No function found for perl binding $func\n";
140 :     }
141 :     }
142 : olson 1.3 }
143 :     }
144 :    
145 :     sub load_track_bindings
146 :     {
147 :     my($self) = @_;
148 :    
149 :     #
150 :     # Load defaults.
151 :     #
152 :    
153 :     my $bindings = "$self->{dir}/Default.bindings";
154 :     my $bfh;
155 :    
156 :     if (!open($bfh, "<$bindings"))
157 :     {
158 : parrello 1.19 warn "Could not open $bindings: $!\n";
159 :     return undef;
160 : olson 1.3 }
161 :    
162 :     local $/ = "//\n";
163 :    
164 :     while (<$bfh>)
165 :     {
166 : parrello 1.19 chomp;
167 : olson 1.3
168 : parrello 1.19 next if $_ eq "";
169 : olson 1.3
170 : parrello 1.19 my($name, @rest) = split(/\n/, $_);
171 :    
172 :     #
173 :     # Read @rest into a binding hash.
174 :     #
175 :     # The associated track is binding_hash{track} if present,
176 :     # otherwise it is assumed to be the same as the name of the binding.
177 :     #
178 :    
179 :     my $binding = {};
180 :     for my $ent (@rest)
181 :     {
182 :     my($k, $v) = split(/\t/, $ent, 2);
183 :     $binding->{$k} = $v;
184 :     }
185 : olson 1.3
186 : parrello 1.19 my $track_name = $binding->{track};
187 :     if (!$track_name)
188 :     {
189 :     $track_name = $name;
190 :     }
191 :    
192 :     $self->{bindings}->{$name} = [$track_name, undef, $name, $binding];
193 : olson 1.4 }
194 :    
195 :     close($bfh);
196 :    
197 :     #
198 :     # Now load the genome-specific binding information, if any.
199 :     #
200 :    
201 :     if (open($bfh, "$self->{dir}/$self->{genome}.bindings"))
202 :     {
203 : parrello 1.19 while (<$bfh>)
204 :     {
205 :     chomp;
206 :    
207 :     next if $_ eq "";
208 :    
209 :     my($name, @rest) = split(/\n/, $_);
210 :    
211 :     #
212 :     # Read @rest into a binding hash.
213 :     #
214 :     # The associated track is binding_hash{track} if present,
215 :     # otherwise it is assumed to be the same as the name of the binding.
216 :     #
217 :    
218 :     my $binding_info = $self->{bindings}->{$name};
219 :     my $binding;
220 :    
221 :     #
222 :     # See if this is a new binding. If it is, create a new hash.
223 :     # If not, use the existing one.
224 :     #
225 :     if (!defined($binding_info))
226 :     {
227 :     $binding = {};
228 :     }
229 :     else
230 :     {
231 :     $binding = $binding_info->[3];
232 :     }
233 :    
234 :     for my $ent (@rest)
235 :     {
236 :     my($k, $v) = split(/\t/, $ent, 2);
237 :     $binding->{$k} = $v;
238 :     }
239 :    
240 :     #
241 :     # See if we are setting or overriding the track.
242 :     #
243 :    
244 :     my $track_name;
245 :     if (!defined($binding_info))
246 :     {
247 :     $track_name = $binding->{track};
248 :     if (!$track_name)
249 :     {
250 :     $track_name = $name;
251 :     }
252 :     }
253 :     else
254 :     {
255 :     $track_name = $binding_info->[0];
256 :     if (defined($binding->{track}))
257 :     {
258 :     $track_name = $binding->{track};
259 :     }
260 :     }
261 :    
262 :     if (!$track_name)
263 :     {
264 :     warn "Track not found for binding '$name'\n";
265 :     next;
266 :     }
267 :    
268 :     $self->{bindings}->{$name} = [$track_name, undef, $name, $binding];
269 :     }
270 :     close($bfh);
271 : olson 1.4 }
272 :    
273 :     #
274 :     # We have now loaded all the binding information.
275 :     #
276 :     # Walk them, filling out the tracks based on the bindings.
277 :     #
278 :    
279 :     while (my($bname, $binding_info) = each(%{$self->{bindings}}))
280 :     {
281 : parrello 1.19 my($track_name, undef, $name, $binding) = @$binding_info;
282 :    
283 :     my $orig_track = $self->{tracks}->{$track_name};
284 : olson 1.4
285 : parrello 1.19 my $track = {};
286 : olson 1.4
287 : parrello 1.19 for my $k (keys(%$orig_track))
288 :     {
289 :     #
290 :     # Do magic variable replacement on values.
291 :     #
292 :    
293 :     my $v = $orig_track->{$k};
294 :     $v =~ s,\${(\w+)},defined($binding->{$1}) ? $binding->{$1} : $1,eg;
295 :     $track->{$k} = $v;
296 :     }
297 :     $binding_info->[1] = $track;
298 :    
299 :     #
300 :     # If we have any perl function args, split them up here using shellwords.
301 :     #
302 :     if (defined($track->{perl_binding_args}))
303 :     {
304 :     $track->{perl_binding_args} = [shellwords($track->{perl_binding_args})];
305 :     }
306 :     else
307 :     {
308 :     $track->{perl_binding_args} = [];
309 :     }
310 : olson 1.4
311 : parrello 1.19 #
312 :     # The track type is updated to be the specific track type we've
313 :     # determined here.
314 :     #
315 : olson 1.4
316 : parrello 1.19 $track->{type} = $name;
317 : olson 1.3 }
318 : olson 1.1 }
319 :    
320 :     sub types_header
321 :     {
322 :     my($self) = @_;
323 :    
324 :     return <<END;
325 :     <?xml version="1.0" standalone="yes"?>
326 :     <!DOCTYPE DASTYPES SYSTEM "http://www.biodas.org/dtd/dastypes.dtd">
327 :     <DASTYPES>
328 : olson 1.3 <GFF version="1.2" summary="yes" href="$self->{url}">
329 : olson 1.1 <SEGMENT>
330 :     END
331 :     ;
332 :    
333 :     }
334 :    
335 : olson 1.3 sub types_footer
336 :     {
337 :     my($self) = @_;
338 :    
339 :     return <<END;
340 :     </SEGMENT>
341 :     </GFF>
342 :     </DASTYPES>
343 :     END
344 :     }
345 :    
346 :     =pod
347 :    
348 :     =head1 all_types
349 :    
350 :     Return XML for the list of all types we support.
351 :    
352 :     This comes from the list of bindings for the current genome (dsn).
353 :    
354 :     =cut
355 :    
356 :     sub all_types
357 :     {
358 :     my($self) = @_;
359 :    
360 :     my $out = "";
361 :    
362 :     #
363 :     # First entry in the bindings list entry is the track.
364 :     #
365 : olson 1.5
366 :     for my $ent (values(%{$self->{bindings}}))
367 : olson 1.3 {
368 : parrello 1.19 my($track_name, $track, $name, $other_info) = @$ent;
369 :     $out .= $self->get_track_xml($track);
370 : olson 1.3 }
371 :     return $out;
372 :     }
373 :    
374 :     sub get_track_xml
375 :     {
376 :     my($self, $track) = @_;
377 :    
378 :     return sprintf(qq(<TYPE id="%s" category="%s" method="%s" source="%s"/>\n),
379 : parrello 1.19 $track->{type},
380 :     $track->{category},
381 :     $track->{method},
382 :     $track->{source});
383 : olson 1.3 }
384 :    
385 :    
386 :     #
387 :     # Features.
388 :     #
389 :    
390 :     sub features_header
391 : olson 1.1 {
392 :     my($self) = @_;
393 :    
394 : olson 1.3 return <<END
395 :     <?xml version="1.0" standalone="yes"?>
396 :     <!DOCTYPE DASGFF SYSTEM "http://www.biodas.org/dtd/dasgff.dtd">
397 :     <DASGFF>
398 :     <GFF version="1.01" href="$self->{url}">
399 :     END
400 :    
401 :     }
402 :    
403 :     sub features_footer
404 :     {
405 :     my($self) = @_;
406 :     return <<END
407 :     </GFF>
408 :     </DASGFF>
409 :     END
410 :    
411 :     }
412 :    
413 :     sub features
414 :     {
415 : olson 1.12 my($self, $segments, $types, $feature_ids) = @_;
416 :    
417 :     #
418 :     # Determine if this is a normal get-the-features-from-a-segment
419 :     # call, or if it was likely performed as a result of a gbrowse search.
420 :     #
421 :     # We detect the latter by noticing if @$segments is empty and @$feature_ids is not.
422 :     #
423 :    
424 :     if (@$segments == () and @$feature_ids > 0)
425 :     {
426 : parrello 1.19 return $self->features_search($types, $feature_ids);
427 : olson 1.12 }
428 :    
429 : olson 1.3
430 :     my %types;
431 :    
432 :     #
433 :     # Build a hash of the types we're looking for.
434 :     #
435 :    
436 :     map { $types{$_}++ } @$types;
437 :    
438 :     my $fig = $self->{fig};
439 :    
440 : olson 1.12 #warn "features: segments=", Dumper($segments);
441 :     #warn " types=", Dumper($types);
442 :     #warn " feature_ids=", Dumper($feature_ids);
443 : olson 1.3
444 :     #
445 :     # Determine what tracks we need to extract by matching against
446 :     # the list of types requested. Also determine, for the tracks
447 :     # we are keeping, if we need to map segments to features.
448 :     #
449 :    
450 :     my $need_features = 0;
451 : olson 1.1
452 : olson 1.4 my @bindings = values(%{$self->{bindings}});
453 :     my @all_tracks = map { $_->[1] } @bindings;
454 : olson 1.3 my @tracks = ();
455 :    
456 : olson 1.4 # warn "Types: ", Dumper(\%types);
457 : olson 1.3
458 : olson 1.5 #
459 :     # We keep track of all the tracks of method Related so that
460 :     # they can be processed first. This is a little ugly, but
461 :     # there is some cross-track processing that needs to be supported.
462 :     #
463 :    
464 :     my @related_tracks;
465 : olson 1.3 for my $track (@all_tracks)
466 :     {
467 : parrello 1.19 # warn "Testing track ", Dumper($track);
468 :     if ($types{$track->{type}})
469 :     {
470 :     # warn "Adding track $track->{type}\n";
471 :    
472 :     if ($track->{method} eq "Related")
473 :     {
474 :     push(@related_tracks, $track);
475 :     }
476 :     else
477 :     {
478 :     push(@tracks, $track);
479 :     }
480 :     $need_features++ if $track->{need_pegs};
481 :     }
482 : olson 1.3 }
483 :    
484 :     #
485 : olson 1.5 # Stuff the related tracks onto the front of the track list.
486 :     #
487 :    
488 :     unshift(@tracks, @related_tracks);
489 :    
490 :     #
491 : olson 1.3 # Walk the segments. For each segment, we emit the segment header,
492 :     # and determine if we have any active tracks. If we do,
493 :     # invoke that code to generate features. Then emit the segment footer.
494 :     #
495 :     # At the start of each segment we determine if we have any peg mapping to
496 :     # determine, and do it once.
497 :     #
498 :    
499 : olson 1.5 # warn "have segs ", Dumper($segments);
500 : olson 1.3 my $output = [];
501 :     for my $segment (@$segments)
502 :     {
503 : parrello 1.19 #
504 :     # Each segment is a list [$reference, $class, $start, $end]
505 :     #
506 :     # $reference is the contig.
507 :     # $class is the class of this segment.
508 :     # $start and $end are the starting and ending points on the contig.
509 :     #
510 :    
511 :     my($seg_id, $class, $start, $end) = @$segment;
512 :    
513 :     my $contig = $seg_id;
514 :     $contig =~ s/--/:/g;
515 :    
516 :     #
517 :     # Determine the start and stop if it was not passed in
518 :     # (an implicit reference to the entire contig).
519 :     #
520 :     if (!defined($start))
521 :     {
522 :     $start = 1;
523 :     $end = $fig->contig_ln($self->{genome}, $contig);
524 :    
525 :     $segment->[2] = $start;
526 :     $segment->[3] = $end;
527 :     }
528 :    
529 :     push(@$output, qq(<SEGMENT id="$seg_id" start="$start" stop="$end" version="1.0">));
530 :    
531 :     my($pegs, $peg_start, $peg_end);
532 :    
533 :     if ($need_features)
534 :     {
535 :     ($pegs, $peg_start, $peg_end) = $self->{fig}->genes_in_region($self->{genome},
536 :     $contig,
537 :     $start,
538 :     $end);
539 :    
540 :     # warn "Checking segment. start=$start end=$end contig=$contig\n";
541 :     # warn " pegs $pegs $peg_start $peg_end\n";
542 :    
543 :     # warn "$fig have pegs for $contig $start $end ", Dumper($pegs);
544 :     #
545 :     # We assign a peg index to each peg, and store it in $self->{peg_index}.
546 :     #
547 :     # This is used for coloring related pegs when enabled.
548 :     #
549 :    
550 :     my $idx = 1;
551 :     for my $peg (@$pegs)
552 :     {
553 :     $self->{peg_index}->{$peg} = $idx++;
554 :     }
555 :     }
556 :    
557 :     for my $track (@tracks)
558 :     {
559 :     my $func = $track->{perl_binding_ref};
560 :    
561 :     if ($func)
562 :     {
563 :     my $track_output = $self->$func($track, $segment, $pegs);
564 :     push(@$output, @$track_output);
565 :     }
566 :     else
567 :     {
568 :     die "No function found for ", Dumper($track);
569 :     }
570 :     }
571 :     push(@$output, qq(</SEGMENT>\n));
572 : olson 1.3 }
573 :     return $output;
574 :     }
575 :    
576 : olson 1.12 #
577 :     # This code is a close replica of the features code. In fact, at some point,
578 :     # it might be possible to factor the common code into two chunks, one that determines
579 :     # a set of segments and a set of pegs and tracks in those segments to show.
580 :     #
581 :     # Until then, we'll keep them separate.
582 :     #
583 :    
584 :     sub features_search
585 :     {
586 :     my($self, $types, $feature_ids) = @_;
587 :    
588 :     #
589 :     # Create a search object and perform the search.
590 :     #
591 :    
592 :     my $fig = $self->{fig};
593 :     my $search = new SproutSearch($fig);
594 :    
595 :     my $search_str = join(" ", @$feature_ids);
596 :     #
597 :     # Elide the trailing * that might be on there, gbrowse gives it meaning..
598 :     #
599 :     $search_str =~ s/\*$//;
600 :    
601 :     my($features, $genomes, $subsystems) = $search->search($search_str);
602 :    
603 :     # warn "Search returns features ", Dumper($features);
604 :    
605 :     #
606 :     # Walk the list of features, determining which contigs the features are
607 :     # on. This will create our list of segments.
608 :     #
609 :    
610 :     my(%contig_pegs, %contig_min, %contig_max);
611 :    
612 :     for my $feature (@$features)
613 :     {
614 : parrello 1.19 my($peg, $anno, $alias) = @$feature;
615 : olson 1.12
616 : parrello 1.19 my $loc = $fig->feature_location($peg);
617 :     my $contig = $fig->contig_of($loc);
618 : olson 1.12
619 : parrello 1.19 my $beg = $fig->beg_of($loc);
620 :     my $end = $fig->end_of($loc);
621 : olson 1.12
622 : parrello 1.19 my $min = \$contig_min{$contig};
623 :     my $max = \$contig_max{$contig};
624 : olson 1.12
625 : parrello 1.19 $$min = $beg if !defined($$min) or $beg < $$min;
626 :     $$max = $end if !defined($$max) or $end > $$max;
627 : olson 1.12
628 : parrello 1.19 push(@{$contig_pegs{$contig}}, [$peg, $loc, $beg, $end]);
629 : olson 1.12 }
630 :    
631 :    
632 :     #
633 :     # Now we can create our segment list from the contigs.
634 :     #
635 :    
636 :     my $segments = [];
637 :    
638 :     for my $contig (keys(%contig_pegs))
639 :     {
640 : parrello 1.19 my $pegs = $contig_pegs{$contig};
641 :     my $min = $contig_min{$contig};
642 :     my $max = $contig_max{$contig};
643 :    
644 :     #
645 :     # The undef here makes this code Sprout-specific.
646 :     #
647 :     my $len = $fig->contig_ln(undef, $contig);
648 :    
649 :     # warn "contig $contig: min=$min max=$max\n";
650 :    
651 :     #
652 :     # Ensure we have at least a 10k wide context.
653 :     #
654 :     if ($max - $min < 10000)
655 :     {
656 :     my $pad = int((10000 - ($max - $min)) / 2);
657 :    
658 :     $min -= $pad;
659 :     $min = 0 if $min < 0;
660 :    
661 :     $max += $pad;
662 :    
663 :     $max = $len if $max > $len;
664 :    
665 :     # warn "pad=$pad min=$min max=$max\n";
666 :     }
667 :    
668 :     # Each segment is a list [$reference, $class, $start, $end]
669 :    
670 :     my $seg_id = $contig;
671 :     $seg_id =~ s/:/--/g;
672 :     push(@$segments, [$seg_id, "CDS", $min, $max, $pegs]);
673 : olson 1.12 }
674 :    
675 :     if (!$types or @$types == 0)
676 :     {
677 : parrello 1.19 $types = ['CDS:curated'];
678 : olson 1.12 }
679 :     else
680 :     {
681 : parrello 1.19 return $self->emit_segments_only_for_match($segments);
682 : olson 1.12 }
683 :     my %types;
684 :    
685 :     #
686 :     # Build a hash of the types we're looking for.
687 :     #
688 :    
689 :     map { $types{$_}++ } @$types;
690 :    
691 : parrello 1.19 $fig = $self->{fig};
692 : olson 1.12
693 :     # warn "features: segments=", Dumper($segments);
694 :     # warn " types=", Dumper($types);
695 :     # warn " feature_ids=", Dumper($feature_ids);
696 :    
697 :     #
698 :     # Determine what tracks we need to extract by matching against
699 :     # the list of types requested. Also determine, for the tracks
700 :     # we are keeping, if we need to map segments to features.
701 :     #
702 :    
703 :     my $need_features = 0;
704 :    
705 :     my @bindings = values(%{$self->{bindings}});
706 :     my @all_tracks = map { $_->[1] } @bindings;
707 :     my @tracks = ();
708 :    
709 :     # warn "Types: ", Dumper(\%types);
710 :    
711 :     #
712 :     # We keep track of all the tracks of method Related so that
713 :     # they can be processed first. This is a little ugly, but
714 :     # there is some cross-track processing that needs to be supported.
715 :     #
716 :    
717 :     my @related_tracks;
718 :     for my $track (@all_tracks)
719 :     {
720 : parrello 1.19 # warn "Testing track ", Dumper($track);
721 :     if ($types{$track->{type}})
722 :     {
723 :     # warn "Adding track $track->{type}\n";
724 :    
725 :     if ($track->{method} eq "Related")
726 :     {
727 :     push(@related_tracks, $track);
728 :     }
729 :     else
730 :     {
731 :     push(@tracks, $track);
732 :     }
733 :     $need_features++ if $track->{need_pegs};
734 :     }
735 : olson 1.12 }
736 :    
737 :     #
738 :     # Stuff the related tracks onto the front of the track list.
739 :     #
740 :    
741 :     unshift(@tracks, @related_tracks);
742 :    
743 :     #
744 :     # Walk the segments. For each segment, we emit the segment header,
745 :     # and determine if we have any active tracks. If we do,
746 :     # invoke that code to generate features. Then emit the segment footer.
747 :     #
748 :     # At the start of each segment we determine if we have any peg mapping to
749 :     # determine, and do it once.
750 :     #
751 :    
752 :     # warn "have segs ", Dumper($segments);
753 :     my $output = [];
754 :     for my $segment (@$segments)
755 :     {
756 : parrello 1.19 #
757 :     # Each segment is a list [$reference, $class, $start, $end, $peg_list]
758 :     #
759 :     # $reference is the contig.
760 :     # $class is the class of this segment.
761 :     # $start and $end are the starting and ending points on the contig.
762 :     #
763 :    
764 :     my($seg_id, $class, $start, $end, $peg_list) = @$segment;
765 :    
766 :     my $pegs = [map { $_->[0] } @$peg_list];
767 :    
768 :     # warn "segment $seg_id has pegs ", Dumper($pegs);
769 :    
770 :     my $contig = $seg_id;
771 :     $contig =~ s/--/:/g;
772 :    
773 :     push(@$output, qq(<SEGMENT id="$seg_id" start="$start" stop="$end" version="1.0">));
774 :    
775 :     my $idx = 1;
776 :     for my $peg (@$pegs)
777 :     {
778 :     $self->{peg_index}->{$peg} = $idx++;
779 :     }
780 :    
781 :     for my $track (@tracks)
782 :     {
783 :     my $func = $track->{perl_binding_ref};
784 :    
785 :     if ($func)
786 :     {
787 :     my $track_output = $self->$func($track, $segment, $pegs);
788 :     push(@$output, @$track_output);
789 :     }
790 :     else
791 :     {
792 :     die "No function found for ", Dumper($track);
793 :     }
794 :     }
795 :     push(@$output, qq(</SEGMENT>));
796 : olson 1.12 }
797 :     return $output;
798 :     }
799 :    
800 :     sub emit_segments_only_for_match
801 :     {
802 :     my($self, $segments) = @_;
803 :    
804 :     my $output = [];
805 :    
806 :     for my $segment (@$segments)
807 :     {
808 : parrello 1.19 #
809 :     # Each segment is a list [$reference, $class, $start, $end, $peg_list]
810 :     #
811 :     # $reference is the contig.
812 :     # $class is the class of this segment.
813 :     # $start and $end are the starting and ending points on the contig.
814 :     #
815 :    
816 :     my($seg_id, $class, $start, $end, $peg_list) = @$segment;
817 :    
818 :     for my $peg_ent (@$peg_list)
819 :     {
820 :     my($peg, $loc, $peg_start, $peg_end) = @$peg_ent;
821 :    
822 :     my $contig = $seg_id;
823 :     $contig =~ s/--/:/g;
824 :    
825 :     my $label = $self->get_peg_name($peg);
826 :    
827 :     if ($peg_start > $peg_end)
828 :     {
829 :     ($peg_start, $peg_end) = ($peg_end, $peg_start);
830 :     }
831 :    
832 :     push(@$output, qq(<SEGMENT id="$seg_id" start="$peg_start" stop="$peg_end" Note="$label" version="1.0">));
833 :     push(@$output, qq(</SEGMENT>));
834 :     }
835 : olson 1.12 }
836 :     return $output;
837 :     }
838 : olson 1.3
839 :    
840 :     #
841 :     # Track binding functions
842 :     #
843 :    
844 :     =pod
845 :    
846 :     =head1 track_cds
847 :    
848 :     Generate the CDS features for a segment
849 :    
850 :     =cut
851 :    
852 :     sub track_cds
853 :     {
854 :     my($self, $track, $segment, $pegs) = @_;
855 :    
856 :     #
857 :     # $track is the track configuration hash.
858 :     # $segments is a list of segments, including the SEED features on that segment.
859 :     # @$pegs is a list of pegs in that segment.
860 :     #
861 :    
862 :     my $fig = $self->{fig};
863 :     my $type = $track->{type};
864 :     my $method = $track->{method};
865 :     my $category = $track->{category};
866 :    
867 : olson 1.9 # warn "track_cds: track=", Dumper($track);
868 :     # warn "track_cds: segment=", Dumper($segment);
869 :     # warn "track_cds: pegs=", Dumper($pegs);
870 :    
871 : olson 1.4 my($feature_type) = @{$track->{perl_binding_args}};
872 :    
873 : olson 1.3 my $features = [];
874 :    
875 :     my($contig, $class, $start, $end) = @$segment;
876 :    
877 :     foreach my $peg (@$pegs)
878 : olson 1.4 {
879 : parrello 1.19 next unless $peg =~ /\.$feature_type\./;
880 :     my $loc = $self->get_feature_location($peg);
881 :     my $start = $fig->beg_of($loc);
882 :     my $stop = $fig->end_of($loc);
883 :     my $orientation;
884 :    
885 :     # warn "Checking peg $peg\n";
886 :    
887 :     if ($start < $stop)
888 :     {
889 :     $orientation = '+';
890 :     }
891 :     else
892 :     {
893 :     $orientation = '-';
894 :     my $t = $stop;
895 :     $stop = $start;
896 :     $start = $t;
897 :     }
898 :    
899 :     my $len = $stop - $start + 1;
900 :     my $phase = $len % 3;
901 :    
902 :     #
903 :     # Determine function and label.
904 :     #
905 :    
906 :     my($func, $name);
907 :    
908 :     $name = $self->get_peg_name($peg);
909 :     $func = $self->get_peg_func($peg);
910 :    
911 :     my $score;
912 :     if (defined($self->{peg_in_related_track}->{$peg}))
913 :     {
914 :     $score = $self->{peg_index}->{$peg};
915 :     }
916 :     else
917 :     {
918 :     $score = 0;
919 :     }
920 :    
921 :     #
922 :     # oh, don't look too closely here.
923 :     #
924 :     # we want to use relative urls here, so that proxying does the right thing.
925 :     #
926 :     # These links are exposed from gbrowse, where the urls look like
927 :     # http://www.nmpdr.org/nmpdr_v3/FIG/gbrowse.cgi/GB_195099.1?ref=195099.1--CP000025;start=925978;stop=966477;nav4=1;plugin=
928 :     #
929 :     # the brwoser doesn't know that /GB_195099.1 isn't a real part of the path, so
930 :     # the usual relative links don't work. we put in .. to make that work.
931 :     #
932 :    
933 :     my $link_url;
934 :     if ($feature_type eq "peg")
935 :     {
936 :     $link_url = "../protein.cgi?prot=$peg&SPROUT=1";
937 :     }
938 :     else
939 :     {
940 :     $link_url = "../feature.cgi?prot=$peg;feature=$peg&SPROUT=1";
941 :     }
942 :     my $link = qq(<LINK href="$link_url">$peg</LINK>);
943 :    
944 :     my $feature = <<END;
945 : olson 1.3 <FEATURE id="Gene:$peg" label="$name">
946 : parrello 1.19 <TYPE id="$type" category="$category" reference="no" subparts="no">$type</TYPE>
947 :     <METHOD id="$method">$method</METHOD>
948 :     <START>$start</START>
949 :     <END>$stop</END>
950 :     <SCORE>$score</SCORE>
951 :     <ORIENTATION>$orientation</ORIENTATION>
952 :     <NOTE>$func</NOTE>
953 :     <PHASE>$phase</PHASE>
954 :     $link
955 : olson 1.3 </FEATURE>
956 :     END
957 :     push(@$features, $feature);
958 :     }
959 :     return $features;
960 :     }
961 :    
962 : olson 1.4
963 : olson 1.6 #
964 :     # Return a track containing features that have a given attribute.
965 :     #
966 :     # If the feature has the attribute, the score contains the value of that attribute.
967 :     #
968 :     sub track_attribute
969 :     {
970 :     my($self, $track, $segment, $pegs) = @_;
971 :    
972 :     #
973 :     # $track is the track configuration hash.
974 :     # $segments is a list of segments, including the SEED features on that segment.
975 :     # @$pegs is a list of pegs in that segment.
976 :     #
977 :    
978 :     my $fig = $self->{fig};
979 :     my $type = $track->{type};
980 :     my $method = $track->{method};
981 :     my $category = $track->{category};
982 :    
983 :     my($attribute) = @{$track->{perl_binding_args}};
984 :    
985 :     my $features = [];
986 :    
987 :     my($contig, $class, $start, $end) = @$segment;
988 :    
989 :     foreach my $peg (@$pegs)
990 :     {
991 : parrello 1.19 # Old mechansim.
992 :     #
993 :     # my $val = $fig->get_attribute($peg, $attribute);
994 :     # next unless defined($val);
995 :     # next unless $val > 0;
996 :     #
997 :     # Now we use feature_attributes, which returns the attributes for a feature
998 :     # as a list of triples [att-name, value, associated-url].
999 :     # RAE: now a list of fourples: [peg, att-name, value, associated-url]
1000 :     #
1001 :    
1002 :     my @atts = $fig->feature_attributes($peg);
1003 :    
1004 :     #
1005 :     # Search @atts for the one we're looking for with a true value.
1006 :     #
1007 :    
1008 :     @atts = grep { $_->[0] eq $attribute and $_->[1] != 0 } @atts;
1009 :    
1010 :     next unless @atts > 0;
1011 :     # RAE: I changed $fig->feature_attributes to return peg, tag, value, url
1012 :     #my $val = $atts[0]->[1];
1013 :     #my $url = $atts[0]->[2];
1014 :    
1015 :     my $val = $atts[0]->[2];
1016 :     my $url = $atts[0]->[3];
1017 :    
1018 :     my $link;
1019 :    
1020 :     if ($url)
1021 :     {
1022 :     $link = qq(<LINK href="$url">$peg</LINK>);
1023 :     }
1024 :    
1025 :     #
1026 :     # Computed val, finish up the rest of the feature computation.
1027 :     #
1028 :    
1029 :     my $loc = $self->get_feature_location($peg);
1030 :     my $start = $fig->beg_of($loc);
1031 :     my $stop = $fig->end_of($loc);
1032 :     my $orientation;
1033 :    
1034 :     $orientation = ".";
1035 :     if ($start > $stop)
1036 :     {
1037 :     my $t = $stop;
1038 :     $stop = $start;
1039 :     $start = $t;
1040 :     }
1041 :    
1042 :     my $len = $stop - $start + 1;
1043 :     my $phase = $len % 3;
1044 :    
1045 :     my $name = $self->get_peg_name($peg);
1046 :    
1047 :     my $feature = <<END;
1048 : olson 1.10 <FEATURE id="Gene:$peg" label="$name">
1049 : parrello 1.19 <TYPE id="$type" category="$category" reference="no" subparts="no">$type</TYPE>
1050 :     <METHOD id="$method">$method</METHOD>
1051 :     <START>$start</START>
1052 :     <END>$stop</END>
1053 :     <SCORE>$val</SCORE>
1054 :     <ORIENTATION>$orientation</ORIENTATION>
1055 :     <PHASE>$phase</PHASE>
1056 :     $link
1057 : olson 1.6 </FEATURE>
1058 :     END
1059 :     push(@$features, $feature);
1060 :     }
1061 :     return $features;
1062 :     }
1063 :    
1064 : olson 1.7 #
1065 :     # Return a track containing pegs grouped into functional clusters.
1066 :     #
1067 :     # Two pegs with the same score are in the same cluster.
1068 :     #
1069 :     sub track_functional_cluster
1070 :     {
1071 :     my($self, $track, $segment, $pegs) = @_;
1072 :    
1073 :     #
1074 :     # $track is the track configuration hash.
1075 :     # $segments is a list of segments, including the SEED features on that segment.
1076 :     # @$pegs is a list of pegs in that segment.
1077 :     #
1078 :    
1079 :     my $fig = $self->{fig};
1080 :     my $type = $track->{type};
1081 :     my $method = $track->{method};
1082 :     my $category = $track->{category};
1083 :    
1084 :     my $features = [];
1085 :    
1086 :     my($contig, $class, $start, $end) = @$segment;
1087 :    
1088 :     my(%pool, %assigned);
1089 :    
1090 :     map { $pool{$_}++ } @$pegs;
1091 : olson 1.14
1092 :     #
1093 :     # Compute the clusters for these pegs.
1094 :     #
1095 :    
1096 :     my %clusters;
1097 :    
1098 :     for my $peg (@$pegs)
1099 :     {
1100 : parrello 1.19 my @cluster = grep { defined($pool{$_}) } $fig->in_cluster_with($peg);
1101 :     $clusters{$peg} = [$peg, int(@cluster), [@cluster]];
1102 : olson 1.14 }
1103 :    
1104 :    
1105 : olson 1.7 my $score = 1;
1106 :     while (%pool > 0)
1107 :     {
1108 : parrello 1.19 #
1109 :     # Scan the remaining pegs, and find the largest cluster.
1110 :     #
1111 :    
1112 :     my @peg_clusters;
1113 :    
1114 :     #
1115 :     # Need to remove pegs that have already been assigned
1116 :     # from the original clusters.
1117 :     #
1118 :    
1119 :     for my $c (@clusters{keys(%pool)})
1120 :     {
1121 :     my($peg, $len, $cluster) = @$c;
1122 :     my @new_cluster = grep { defined($pool{$_}) } @$cluster;
1123 :     push(@peg_clusters, [$peg, int(@new_cluster), [@new_cluster]]);
1124 :     }
1125 :    
1126 :     # print "peg cluster len=", int(@peg_clusters), "\n";
1127 :    
1128 :     my $winner = (sort { $b->[1] <=> $a->[1] } @peg_clusters)[0];
1129 :     # print "Winner=", Dumper($winner);
1130 :    
1131 :     my($peg, $count, $cluster) = @$winner;
1132 :    
1133 :     delete $pool{$peg};
1134 :    
1135 :     # print "Cluster is @$cluster\n";
1136 :    
1137 :     for my $c (@$cluster)
1138 :     {
1139 :     if ($pool{$c})
1140 :     {
1141 :     delete $pool{$c};
1142 :     $assigned{$c} = $score;
1143 :     }
1144 :     }
1145 :     $assigned{$peg} = $score;
1146 :     $score++;
1147 : olson 1.7 }
1148 :    
1149 :     foreach my $peg (@$pegs)
1150 :     {
1151 : parrello 1.19 my $val = $assigned{$peg};
1152 :     next unless defined($val);
1153 :     next unless $val > 0;
1154 :    
1155 :     my $loc = $self->get_feature_location($peg);
1156 :     my $start = $fig->beg_of($loc);
1157 :     my $stop = $fig->end_of($loc);
1158 :     my $orientation;
1159 :    
1160 :     if ($start < $stop)
1161 :     {
1162 :     $orientation = '+';
1163 :     }
1164 :     else
1165 :     {
1166 :     $orientation = '-';
1167 :     my $t = $stop;
1168 :     $stop = $start;
1169 :     $start = $t;
1170 :     }
1171 :    
1172 :     my $len = $stop - $start + 1;
1173 :     my $phase = $len % 3;
1174 :    
1175 :     my $name = $self->get_peg_name($peg);
1176 :    
1177 :     my $feature = <<END;
1178 : olson 1.10 <FEATURE id="Gene:$peg" label="$name">
1179 : parrello 1.19 <TYPE id="$type" category="$category" reference="no" subparts="no">$type</TYPE>
1180 :     <METHOD id="$method">$method</METHOD>
1181 :     <START>$start</START>
1182 :     <END>$stop</END>
1183 :     <SCORE>$val</SCORE>
1184 :     <ORIENTATION>$orientation</ORIENTATION>
1185 :     <PHASE>$phase</PHASE>
1186 : olson 1.7 </FEATURE>
1187 :     END
1188 :     push(@$features, $feature);
1189 :     }
1190 :     return $features;
1191 :     }
1192 :    
1193 : olson 1.6
1194 : olson 1.3 sub track_contig
1195 :     {
1196 :     my($self, $track, $segment, $pegs) = @_;
1197 :    
1198 :     #
1199 :     # $track is the track configuration hash.
1200 :     # $segments is a list of segments, including the SEED features on that segment.
1201 :     # @$pegs is a list of pegs in that segment.
1202 :     #
1203 :    
1204 :     my $fig = $self->{fig};
1205 :     my $type = $track->{type};
1206 :     my $method = $track->{method};
1207 :     my $category = $track->{category};
1208 :    
1209 :     my $features = [];
1210 :    
1211 :     my($contig, $class, $start, $end) = @$segment;
1212 :    
1213 :     my $feature = <<EOF;
1214 :     <FEATURE id="sequence:$contig/1" label="$contig">
1215 : parrello 1.19 <TYPE id="$type" category="$category" reference="no" subparts="no">$type</TYPE>
1216 :     <METHOD id="Component">Component</METHOD>
1217 :     <START>$start</START>
1218 :     <END>$end</END>
1219 :     <SCORE>0</SCORE>
1220 :     <ORIENTATION>+</ORIENTATION>
1221 :     <PHASE>0</PHASE>
1222 :     <LINK href="http://mcs.anl.gov">$contig</LINK>
1223 :     <NOTE>"Test Note"</NOTE>
1224 :     <GROUP id="Sequence:$contig" type="Sequence" />
1225 : olson 1.3 </FEATURE>
1226 :     EOF
1227 :    
1228 :     push(@$features, $feature);
1229 :    
1230 :     return $features;
1231 :     }
1232 :    
1233 : olson 1.8
1234 :     #
1235 :     # Map a set of aliases to a canonical name, based on a hierarchy
1236 :     # of naming preferences.
1237 :     #
1238 : olson 1.3 sub anti_alias() {
1239 : olson 1.8 my(@aliases) = @_;
1240 :    
1241 :     my(%names);
1242 :    
1243 :     for $_ (@aliases)
1244 :     {
1245 : parrello 1.19 #
1246 :     # Match the various forms of names, and mark in %matches.
1247 :     #
1248 :    
1249 :     if (/^gi\|/) {
1250 :     $names{"gi"} = $_;
1251 :     }
1252 :     elsif (/^uni\|/)
1253 :     {
1254 :     $names{"uni"} = $_;
1255 :     }
1256 :     elsif (/^[a-z]{2,3}[A-Z][0-9]?/)
1257 :     {
1258 :     $names{"gene_name"} = $_;
1259 :     }
1260 :     elsif (/^[A-Za-z]{1,2}\d+$/)
1261 :     {
1262 :     $names{seq} = $_;
1263 :     }
1264 :     elsif (/^[A-Z]{2}_/) {
1265 :     $names{"ref"} = $_;
1266 :     }
1267 :     elsif (/^fig\|/)
1268 :     {
1269 :     $names{fig} = $_;
1270 :     }
1271 : olson 1.8
1272 :     }
1273 :    
1274 :     # warn "generated map ", join(" ", %names), "\n";
1275 :    
1276 :     my(@prec) = qw(gene_name seq uni ref fig kegg);
1277 :    
1278 :     my @matches = grep { $_ ne "" } @names{@prec};
1279 :    
1280 :     # warn "mapped to @matches\n";
1281 :     return $matches[0];
1282 : olson 1.1 }
1283 :    
1284 : olson 1.4 =pod
1285 :    
1286 :     4. Bruce and Bob need to help implement a "track" for showing each
1287 :     genome related to the given one. The logic for this track is based
1288 :     on a subroutine that takes as input
1289 :    
1290 :     the given genome X
1291 :     the given contig
1292 :     the begin (leftmost coord) of the map B
1293 :     the end of the map E
1294 :     a second genome Y
1295 :    
1296 :     The routine must return a set of features with coordinates in the
1297 :     range begin-end (i.e., coordinates in the given genome's contig)
1298 :    
1299 :     This routine should have roughly this logic:
1300 :    
1301 :     a. Let S1 be the features in X in the region B-E
1302 :     b. Let O1 be the BBHs of those features in Y
1303 :     c. Let O1r be the largest subset of features in O1 that occur
1304 :     within a single contig and within a distance of E-B.
1305 :     d. pick the bounding pairs (of S1 entries with O1r entries).
1306 :     Align the midpoints of each bounding pair. Now find all
1307 :     genes in Y that fall within the appropriate range around
1308 :     the aligned midpoints (remember the genomes may have
1309 :     reversed directions, so pick the appropriate coords).
1310 :     e. adjust all the coordinates to fall in the given genomes
1311 :     range.
1312 :     f. return the features.
1313 :    
1314 :     =cut
1315 :    
1316 :     sub track_related
1317 :     {
1318 :     my($self, $track, $segment, $pegs) = @_;
1319 :    
1320 :     # warn "Track_related ", Dumper(\@_);
1321 : olson 1.5 # warn "Segment: ", Dumper($segment);
1322 : olson 1.4 # warn "pegs ", Dumper($pegs);
1323 :     #
1324 :     # $track is the track configuration hash.
1325 :     # $segments is a list of segments, including the SEED features on that segment.
1326 :     # @$pegs is a list of pegs in that segment.
1327 :     #
1328 : olson 1.5
1329 : olson 1.4 my $fig = $self->{fig};
1330 :     my $type = $track->{type};
1331 :     my $method = $track->{method};
1332 :     my $category = $track->{category};
1333 :    
1334 :     my($y) = @{$track->{perl_binding_args}};
1335 :    
1336 :     my $features = [];
1337 :    
1338 :     my($contig, $class, $start, $end) = @$segment;
1339 :    
1340 :     my $s1 = $pegs;
1341 :    
1342 :     my %contig;
1343 :    
1344 :     #
1345 :     # Remember the index of the related peg.
1346 :     #
1347 :    
1348 :     my %related_index;
1349 : olson 1.5 my %related_peg_to_original_peg;
1350 :    
1351 :     my %org_peg_info;
1352 : olson 1.4
1353 :     my($my_peg_min, $my_peg_max);
1354 :    
1355 :     for my $peg (@$s1)
1356 :     {
1357 : parrello 1.19 # print "Handling $peg\n";
1358 :    
1359 :     my $peg_bbhs = $fig->bbh_list($y, [$peg]);
1360 :    
1361 :     # warn "bbhs for $peg: ", Dumper ($peg_bbhs);
1362 : olson 1.5
1363 : parrello 1.19 #
1364 :     # Is it one peg or a list?
1365 :     #
1366 :     my @o1;
1367 :     my $tmp = $peg_bbhs->{$peg};
1368 :     if (ref($tmp) eq "ARRAY")
1369 :     {
1370 :     @o1 = @$tmp;
1371 :     }
1372 :     elsif ($tmp ne "")
1373 :     {
1374 :     @o1 = ($tmp);
1375 :     }
1376 : olson 1.5
1377 : parrello 1.19 for my $o (@o1)
1378 :     {
1379 :     $related_peg_to_original_peg{$o} = $peg;
1380 :    
1381 :     my $loc = $self->get_feature_location($o);
1382 :     # warn "o=$o loc=$loc\n";
1383 :    
1384 :     if ($loc eq '')
1385 :     {
1386 :     #
1387 :     # BUG.
1388 :     #
1389 :     # We hit a feature not in the db. it may have been deleted
1390 :     # (the bug as of 2005-0830 is that we renamed PEGs and
1391 :     # didnt recompute bbhs).
1392 :     #
1393 :    
1394 :     next;
1395 :    
1396 :     }
1397 :    
1398 :     my $b = $fig->beg_of($loc);
1399 :     my $e = $fig->end_of($loc);
1400 :     my $contig = $fig->contig_of($loc);
1401 :    
1402 :     #
1403 :     # Determine min and max and stash those for later use.
1404 :     #
1405 :    
1406 :     my($min, $max);
1407 :     if ($b < $e)
1408 :     {
1409 :     $min = $b;
1410 :     $max = $e;
1411 :     }
1412 :     else
1413 :     {
1414 :     $min = $e;
1415 :     $max = $b;
1416 :     }
1417 :    
1418 :     push(@{$contig{$contig}}, [$o, $b, $e, $min, $max]);
1419 :    
1420 :     $related_index{$o} = $self->{peg_index}->{$peg};
1421 :     }
1422 :    
1423 :     #
1424 :     # While we're scanning, compute min/max extent of our pegs.
1425 :     #
1426 :    
1427 :     my $ploc = $self->get_feature_location($peg);
1428 :    
1429 :     #
1430 :     # Work around the bug here too.
1431 :     #
1432 :     next if $ploc eq '';
1433 :    
1434 :     my $b = $fig->beg_of($ploc);
1435 :     my $e = $fig->end_of($ploc);
1436 : olson 1.4
1437 : parrello 1.19 my($max, $min);
1438 :     if ($b < $e)
1439 :     {
1440 :     $min = $b;
1441 :     $max = $e;
1442 :     }
1443 :     else
1444 :     {
1445 :     $min = $e;
1446 :     $max = $b;
1447 :     }
1448 :    
1449 :     $org_peg_info{$peg} = [$b, $e, $min, $max];
1450 :    
1451 :     if (!defined($my_peg_min))
1452 :     {
1453 :     $my_peg_min = $min;
1454 :     $my_peg_max = $max;
1455 :     }
1456 :     else
1457 :     {
1458 :     $my_peg_min = $min if $min < $my_peg_min;
1459 :     $my_peg_max = $max if $max > $my_peg_max;
1460 :     }
1461 :    
1462 : olson 1.4 }
1463 :    
1464 : olson 1.5 # warn "Contig info ", Dumper(\%contig);
1465 :    
1466 : olson 1.4 #
1467 :     # Pick the contig with the largest set of pegs.
1468 :     #
1469 :    
1470 :     #
1471 : olson 1.5 # Flatten contig map, generating an array @x
1472 :     # of pairs [contig-name, list of [peg, beg, end, min, max] tuples]
1473 : olson 1.4 #
1474 :     my @x = map { [$_, $contig{$_}] } keys(%contig);
1475 :    
1476 : olson 1.5 #
1477 :     # If we found no bbhs, just return.
1478 :     #
1479 : olson 1.4 return [] if @x == 0;
1480 :    
1481 :     #
1482 : olson 1.5 # Sort by number of hits per contig.
1483 : olson 1.4 #
1484 :    
1485 : olson 1.5 @x = sort { @{$b->[1]} <=> @{$a->[1]} } @x;
1486 :    
1487 :     #
1488 :     # And select the largest.
1489 :     #
1490 : olson 1.4
1491 :     my($o1r_contig, $o1r_pegs) = @{$x[0]};
1492 :    
1493 :     # print Dumper($o1r_contig, $o1r_pegs);
1494 :    
1495 :     #
1496 :     # We now have a candidate contig and a set of pegs on that contig.
1497 :     # The peg set is a list of tuples [peg, start, stop, min, max].
1498 :     # start > stop if on the - strand.
1499 :     #
1500 :    
1501 :     #
1502 :     # Sort the candidates by lowest-position.
1503 :     #
1504 :    
1505 :     my @csorted = sort { $a->[3] <=> $b->[3] } @$o1r_pegs;
1506 :    
1507 :     # print Dumper(\@csorted);
1508 :    
1509 :     #
1510 :     # Now scan them in order, creating an array that contains
1511 :     # for each element, the number of pegs in a span-wide section that
1512 :     # overlap.
1513 :     #
1514 :    
1515 :     my @overlaps;
1516 :    
1517 :     #
1518 :     # Remember, $start and $end bound the contig we're searching from.
1519 :     #
1520 :    
1521 :     my $search_size = $end - $start + 1;
1522 :    
1523 :     for my $base_idx (0..@csorted - 1)
1524 :     {
1525 : parrello 1.19 my $dist = 0;
1526 :     my $base = $csorted[$base_idx];
1527 :     my $reach = $base->[4] + $search_size;
1528 :     my $end_idx;
1529 :    
1530 :     # warn "Base: $base_idx reach=$reach base=@$base\n";
1531 :    
1532 :     for my $extent_idx ($base_idx + 1.. @csorted - 1)
1533 :     {
1534 :     my $extent = $csorted[$extent_idx];
1535 :    
1536 :     # warn "Look: $extent_idx extent=@$extent\n";
1537 :    
1538 :     #
1539 :     # If this peg is out of range, we're done.
1540 :     #
1541 :     if ($extent->[3] > $reach)
1542 :     {
1543 :     last;
1544 :     }
1545 :     $end_idx = $extent_idx;
1546 :     }
1547 :     if ($end_idx)
1548 :     {
1549 :     # print "Reached end $end_idx\n";
1550 :     push(@overlaps, [$base_idx, $end_idx, $end_idx - $base_idx]);
1551 :    
1552 :     #
1553 :     # Special case: If we cover all of the pegs, we don't have to search any more
1554 :     #
1555 :    
1556 :     if ($base_idx == 0 and $end_idx == (@csorted - 1))
1557 :     {
1558 :     last;
1559 :     }
1560 :     }
1561 : olson 1.4 }
1562 :    
1563 :     #
1564 :     # Overlaps is now a list of tuples (starting peg, ending peg, difference);
1565 :     # Sort on the number of pegs there. The maximal size is the set
1566 :     # we will map to.
1567 :     #
1568 :    
1569 :     my @sorted_overlaps = sort { $b->[2] <=> $a->[2] } @overlaps;
1570 :    
1571 : olson 1.5 return [] if @sorted_overlaps == 0;
1572 :    
1573 : olson 1.4 my ($start_peg, $end_peg) = @{$sorted_overlaps[0]};
1574 :    
1575 :     # print "got $start_peg $end_peg ", Dumper(\@sorted_overlaps);
1576 :    
1577 :     #
1578 :     # We now have the set of target pegs.
1579 : olson 1.5 # Compute the midpoint of this set of pegs, and the midpoint
1580 :     # of the related pegs. This defines the translation
1581 : olson 1.4 # between the contigs.
1582 :     #
1583 :    
1584 :     my @related_pegs = @csorted[$start_peg .. $end_peg];
1585 :    
1586 :     # print Dumper(\@related_pegs);
1587 :    
1588 :     #
1589 :     # Scan our pegs, searching for the min and max points.
1590 :     #
1591 :    
1592 : olson 1.5 my($rel_min, $rel_max, $org_min, $org_max);
1593 :     $rel_min = $org_min = 1e10;
1594 :     $org_max = $rel_max = -1;
1595 : olson 1.4
1596 : olson 1.5 my $flip;
1597 :    
1598 :     for my $rel_info (@related_pegs)
1599 : olson 1.4 {
1600 : parrello 1.19 my($rel_peg, $rel_b, $rel_e, $rel_pmin, $rel_pmax) = @$rel_info;
1601 :     my $org_peg = $related_peg_to_original_peg{$rel_peg};
1602 :    
1603 :     my($org_b, $org_e, $org_pmin, $org_pmax) = @{$org_peg_info{$org_peg}};
1604 :    
1605 :     print "Rel peg: @$rel_info\n";
1606 :     print "Org peg: $org_peg @{$org_peg_info{$org_peg}}\n";
1607 :    
1608 :     if ($rel_pmin < $rel_min)
1609 :     {
1610 :     $rel_min = $rel_pmin;
1611 :     $org_min = $org_pmin;
1612 :     }
1613 :    
1614 :     if ($rel_pmax > $rel_max)
1615 :     {
1616 :     $rel_max = $rel_pmax;
1617 :     $org_max = $org_pmax;
1618 :     }
1619 :    
1620 :     #$org_min = $org_pmin if $org_pmin < $org_min;
1621 :     #$org_max = $org_pmax if $org_pmax > $org_max;
1622 :    
1623 :     $self->{peg_in_related_track}->{$org_peg}++;
1624 :    
1625 :     #
1626 :     # Check one peg to see if we need to flip.
1627 :     #
1628 :    
1629 :     if (!defined($flip))
1630 :     {
1631 :     my $rel_strand = $rel_b > $rel_e;
1632 :     my $org_strand = $org_b > $org_e;
1633 :     $flip = ($rel_strand == $org_strand) ? 1 : -1;
1634 :     }
1635 : olson 1.4 }
1636 :    
1637 : olson 1.6 my $my_center = int(($org_min + $org_max) / 2);
1638 :     my $rel_center = int(($rel_min + $rel_max) / 2);
1639 : olson 1.4
1640 :     #
1641 :     # The centers induce the mapping of peg locations.
1642 :     #
1643 :     # Given a location in related-org coords L,
1644 :     # the corresponding location in my-org coords
1645 :     # is (L - $rel_center) + $my_center
1646 :     # or L + ($my_center - $rel_center)
1647 :     #
1648 : olson 1.5 # However, in order to support flip, we need this:
1649 :     #
1650 :     # $org = $my_center + flip * (L - $rel_center)
1651 :     #
1652 :     # To map backwards:
1653 :     #
1654 :     # L = $flip * ($org - $my_center) + $rel_center
1655 : olson 1.4 #
1656 :    
1657 :     my $offset = int($my_center - $rel_center);
1658 :    
1659 : olson 1.5
1660 :     #
1661 :     # Map from related coords to orig coords
1662 :     #
1663 :     my $map = sub {
1664 : parrello 1.19 my($val) = @_;
1665 :     return $my_center + $flip * ($val - $rel_center);
1666 : olson 1.5 };
1667 :    
1668 :     #
1669 :     # Map from org coords to related coords
1670 :     #
1671 :     my $revmap = sub {
1672 : parrello 1.19 my($val) = @_;
1673 :     return ($val - $my_center) + $rel_center;
1674 : olson 1.5 };
1675 :    
1676 : olson 1.4 #
1677 :     # Now walk the pegs in the area corresponding to the chosen contig,
1678 :     # returning translated features
1679 :     #
1680 :    
1681 : olson 1.5
1682 :     # warn "my_center=$my_center rel_center=$rel_center offset=$offset\n";
1683 : olson 1.4 my ($rel_pegs, $rel_start, $rel_end) = $self->{fig}->genes_in_region($y,
1684 : parrello 1.19 $o1r_contig,
1685 :     &$revmap($start),
1686 :     &$revmap($end));
1687 : olson 1.4
1688 : olson 1.5 # warn Dumper($rel_pegs);
1689 :     # warn Dumper(\%related_index);
1690 :     # warn Dumper($self);
1691 : olson 1.4 for my $peg (@$rel_pegs)
1692 :     {
1693 : parrello 1.19 next if $fig->ftype($peg) ne "peg";
1694 :     my $loc = $self->get_feature_location($peg);
1695 :    
1696 :     my $pstart = $fig->beg_of($loc);
1697 :     my $pstop = $fig->end_of($loc);
1698 :    
1699 :     my($orientation, $max, $min);
1700 : olson 1.4
1701 : parrello 1.19 if ($pstart < $pstop)
1702 :     {
1703 :     $max = $pstop;
1704 :     $min = $pstart;
1705 :     $orientation = "+";
1706 :     }
1707 :     else
1708 :     {
1709 :     $max = $pstart;
1710 :     $min = $pstop;
1711 :     $orientation = "-";
1712 :     }
1713 :    
1714 :     my $len = $max - $min + 1;
1715 :     my $name = $self->get_peg_name($peg);
1716 : olson 1.4
1717 :    
1718 : parrello 1.19 if ($flip == -1)
1719 :     {
1720 :     $orientation = $orientation eq "-" ? "+" : "-";
1721 :     }
1722 :     $pstart = &$map($pstart);
1723 :     $pstop = &$map($pstop);
1724 :    
1725 :     my $score = $related_index{$peg};
1726 :     $score = 0 unless defined($score);
1727 :     my $phase = $len % 3;
1728 :    
1729 :     #
1730 :     # This computes links that go to the SEED protein page.
1731 :     #
1732 :    
1733 :     my $protein_url = "../protein.cgi?prot=$peg&SPROUT=1";
1734 :    
1735 :     #
1736 :     # Compute a link back to gbrowse, but on *this* organism instead
1737 :     # of the reference. We keep the same zoom factor by mapping
1738 :     # the endpoints of the reference contig area to this org's contig.
1739 :     #
1740 :    
1741 :     {
1742 :     my $genome;
1743 :     $genome = $1 if $peg =~ /fig\|(\d+\.\d+)/;
1744 :    
1745 :     my $lstart = &$revmap($start);
1746 :     my $lend = &$revmap($end);
1747 :     my $seg = $o1r_contig;
1748 :     $seg =~ s/:/--/g;
1749 :     $protein_url = "../gbrowse.cgi/GB_$genome?ref=$seg&start=$lstart&stop=$lend";
1750 :     }
1751 :    
1752 :     my $link = qq(<LINK href="$protein_url">$name</LINK>);
1753 : olson 1.6
1754 : parrello 1.19 my $feature = <<END;
1755 : olson 1.4 <FEATURE id="Gene:$peg" label="$name">
1756 : parrello 1.19 <TYPE id="$type" category="$category" reference="no" subparts="no">$type</TYPE>
1757 :     <METHOD id="$method">$method</METHOD>
1758 :     <START>$pstart</START>
1759 :     <END>$pstop</END>
1760 :     <SCORE>$score</SCORE>
1761 :     <ORIENTATION>$orientation</ORIENTATION>
1762 :     <PHASE>$phase</PHASE>
1763 :     $link
1764 : olson 1.4 </FEATURE>
1765 :     END
1766 :     push(@$features, $feature);
1767 : parrello 1.19
1768 : olson 1.4 }
1769 :    
1770 :     return $features;
1771 :     }
1772 :    
1773 : olson 1.10 #
1774 :     # Some utility routines.
1775 :     #
1776 :    
1777 :     sub get_peg_func
1778 :     {
1779 :     my($self, $peg) = @_;
1780 :    
1781 :     if (my $cached = $self->{cache_peg_func}->{$peg})
1782 :     {
1783 : parrello 1.19 return $cached;
1784 : olson 1.10 }
1785 :    
1786 :     my $func = $self->{fig}->function_of($peg);
1787 :     $self->{cache_peg_func}->{$peg} = $func;
1788 :     return $func;
1789 :     }
1790 :    
1791 :    
1792 :     sub get_peg_name
1793 :     {
1794 :     my($self, $peg) = @_;
1795 :    
1796 :     if (my $cached = $self->{cache_peg_name}->{$peg})
1797 :     {
1798 : parrello 1.19 return $cached;
1799 : olson 1.10 }
1800 :    
1801 :     my @aliases = $self->{fig}->feature_aliases($peg);
1802 :     my $name = &anti_alias(@aliases);
1803 :     if (!$name) {
1804 : parrello 1.19 $name = $peg;
1805 :     $name =~ s,^fig\|\d+\.\d+\.peg\.,,;
1806 : olson 1.10 }
1807 :     $self->{cache_peg_name}->{$peg} = $name;
1808 :    
1809 :     return $name;
1810 :     }
1811 :    
1812 : olson 1.11 sub get_feature_location
1813 :     {
1814 :     my($self, $peg) = @_;
1815 : olson 1.10
1816 : olson 1.11 if (my $cached = $self->{cache_feature_location}->{$peg})
1817 :     {
1818 : parrello 1.19 return $cached;
1819 : olson 1.11 }
1820 :    
1821 :     my $loc = $self->{fig}->feature_location($peg);
1822 :     $self->{cache_feature_location}->{$peg} = $loc;
1823 :     return $loc;
1824 :     }
1825 : olson 1.10
1826 : olson 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3