[Bio] / FigWebServices / sgv.cgi Repository:
ViewVC logotype

Annotation of /FigWebServices/sgv.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (view) (download)

1 : overbeek 1.1 # -*- perl -*-
2 :     #
3 :     # 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 :     #
18 :    
19 :    
20 :     use SeedHTML;
21 :     use strict;
22 :     use SeedEnv;
23 :     use SeedV;
24 :     use ProtSims;
25 :     use GenoGraphics;
26 :    
27 :     use CGI;
28 :     my $cgi = new CGI;
29 :    
30 :     if (0)
31 :     {
32 :     my $VAR1;
33 :     eval(join("",`cat /tmp/sgv_parms`));
34 :     $cgi = $VAR1;
35 :     # print STDERR &Dumper($cgi);
36 :     }
37 :    
38 :     if (0)
39 :     {
40 :     print $cgi->header;
41 :     my @params = $cgi->param;
42 :     print "<pre>\n";
43 :     foreach $_ (@params)
44 :     {
45 :     print "$_\t:",join(",",$cgi->param($_)),":\n";
46 :     }
47 :    
48 :     if (0)
49 :     {
50 :     if (open(TMP,">/tmp/sgv_parms"))
51 :     {
52 :     print TMP &Dumper($cgi);
53 :     close(TMP);
54 :     }
55 :     }
56 :     exit;
57 :     }
58 :    
59 :     my $html = [];
60 :     unshift @$html, "<TITLE>Simple Genome Viewer</TITLE>\n";
61 :    
62 :     my $dir = $cgi->param('dir');
63 :     if ((! $dir) || ($dir !~ /\d+\.\d+$/))
64 :     {
65 :     push(@$html,$cgi->h2('You must set dir= to a SEED-format genome directory with a path ending in the genome ID'));
66 :     &SeedHTML::show_page($cgi,$html);
67 :     exit;
68 :     }
69 :    
70 :     my $request = $cgi->param('request');
71 :     if (! $request)
72 :     {
73 :     &start_computing_reference_genome_set($cgi,$dir,$html);
74 :     &make_subsystem_index($cgi,$html,$dir);
75 :     &basic_query($cgi,$html);
76 :     }
77 :     elsif ($request eq 'basic')
78 :     {
79 :     &basic_query($cgi,$html);
80 :     }
81 :     elsif ($request eq 'id')
82 :     {
83 :     &process_id($cgi,$html);
84 :     }
85 :     elsif ($request eq 'features')
86 :     {
87 :     &process_feature_search($cgi,$html);
88 :     }
89 :     elsif ($request eq 'feature')
90 :     {
91 :     &process_feature_display($cgi,$html);
92 :     }
93 :     elsif ($request eq 'subsystems')
94 :     {
95 :     &process_subsystems_search($cgi,$html);
96 :     }
97 :     elsif ($request eq 'peg2subsystems')
98 :     {
99 :     &process_peg2subsystems_search($cgi,$html);
100 :     }
101 :     elsif ($request eq 'compare')
102 :     {
103 :     &comparison($cgi,$html);
104 :     }
105 :     elsif ($request eq 'corresponding')
106 :     {
107 :     &process_corr_search($cgi,$html);
108 :     }
109 :     elsif ($request =~ 'query_only')
110 :     {
111 :     &process_query_only_search($cgi,$html);
112 :     }
113 :     elsif ($request eq 'reference_only')
114 :     {
115 :     &process_ref_only_search($cgi,$html);
116 :     }
117 :    
118 :     &SeedHTML::show_page($cgi,$html);
119 :    
120 :    
121 :     sub make_subsystem_index {
122 :     my($cgi,$html,$dir) = @_;
123 :    
124 :     my %ss = map { chomp; my($subsys,$var) = split(/\t/,$_); (($var ne "-1") && ($var ne 0)) ? ($subsys => $var) : () }
125 :     `cat $dir/Subsystems/subsystems`;
126 :    
127 :     my $sapObject = SAPserver->new;
128 :     my $ssH = $sapObject->classification_of( -ids => [keys(%ss)]);
129 :     open(INDEX,"| sort > $dir/Subsystems/subsystems.index") || die "could not open $dir/Subsystems/subsystems.index";
130 :     foreach $_ (`cat $dir/Subsystems/bindings`)
131 :     {
132 :     chomp;
133 :     my($subsys,$role,$peg) = split(/\t/,$_);
134 :     my $class = ($_ = $ssH->{$subsys}) ? join("; ",@$_) : "";
135 :     print INDEX join("\t",($class,$subsys,$role,$ss{$subsys},$peg)),"\n";
136 :     }
137 :     close(INDEX);
138 :     }
139 :    
140 :     sub id_search_form {
141 :     my($cgi,$dir,$html) = @_;
142 :    
143 :     my $queryG = $cgi->param('genome');
144 :     push(@$html,$cgi->start_form(-action => "sgv.cgi"),
145 :     '<br><b>Get Protein Page for FIG ID, or ACH Page for non-FIG IDs</b><br><br> ',
146 :     $cgi->textfield(-name => 'id', -size => 20),
147 :     $cgi->submit('go'),
148 :     $cgi->hidden(-name => 'request', -value => 'id', -override => 1),
149 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
150 :     $cgi->end_form,
151 :     $cgi->hr,$cgi->hr);
152 :     }
153 :    
154 :     sub start_computing_reference_genome_set {
155 :     my($cgi,$dir,$html) = @_;
156 :    
157 :     if (! -d "$dir/CorrToReferenceGenomes")
158 :     {
159 :     system "$FIG_Config::bin/get_neighbors_and_corr_to_ref $dir &";
160 :     }
161 :     }
162 :    
163 :     sub basic_query {
164 :     my($cgi,$html) = @_;
165 :    
166 :     my $queryG = $cgi->param('genome');
167 :     my $dir = $cgi->param('dir');
168 :    
169 :     &id_search_form($cgi,$dir,$html);
170 :    
171 :     push(@$html,$cgi->start_form(-action => "sgv.cgi"),
172 :     '<b>Query Features in Genome</b>: ',
173 :     $cgi->textfield(-name => 'pattern', -size => 30),
174 :     $cgi->submit('go'),
175 :     $cgi->hidden(-name => 'request', -value => 'features', -override => 1),
176 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
177 :     $cgi->end_form,
178 :     $cgi->hr,
179 :    
180 :     $cgi->start_form(-action => "sgv.cgi"),
181 :     '<b>Query Subsystems in Genome</b>: ',
182 :     $cgi->textfield(-name => 'pattern', -size => 30),
183 :     $cgi->submit('go'),
184 :     $cgi->hidden(-name => 'request', -value => 'subsystems', -override => 1),
185 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
186 :     $cgi->end_form,
187 :     $cgi->hr
188 :     );
189 :    
190 :     my $cache = "$dir/CorrToReferenceGenomes";
191 :     opendir(CACHE,$cache) || die "could not open $cache";
192 :     my @refG = map { ((-s "$cache/$_") && ($_ =~ /^(\d+\.\d+$)/)) ? $1 : () } readdir(CACHE);
193 :     closedir(CACHE);
194 :    
195 :     if (@refG > 0)
196 :     {
197 :     my $sapObject = SAPserver->new();
198 :     my($refG,$labels) = &build_labels(\@refG,$sapObject);
199 :    
200 :     push(@$html,$cgi->start_form(-action => "sgv.cgi"),
201 :     '<b>Compare Genome Against Reference Genome</b>: ',
202 :     $cgi->scrolling_list( -name => 'reference',
203 :     -values => $refG,
204 :     -labels => $labels,
205 :     -size => 4),
206 :     $cgi->submit('go'),
207 :     $cgi->hidden(-name => 'request', -value => "compare", -override => 1),
208 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
209 :     $cgi->end_form
210 :     );
211 :     }
212 :     }
213 :    
214 :     sub build_labels {
215 :     my($genomes,$sapObject) = @_;
216 :    
217 :     my $genomesH = $sapObject->all_genomes(-complete => 1);
218 :     my $metricsH = $sapObject->genome_metrics(-ids => $genomes);
219 :     my %labels = map { my($contigs,$sz) = @{$metricsH->{$_}};
220 :     my $lab = $genomesH->{$_} . " ($_): $sz bp, $contigs contigs";
221 :     $_ => $lab
222 :     }
223 :     @$genomes;
224 :    
225 :     my @genomes = sort { lc($labels{$a}) cmp lc($labels{$b}) } @$genomes;
226 :    
227 :     return (\@genomes,\%labels);
228 :     }
229 :    
230 :     sub process_feature_search {
231 :     my($cgi,$html) = @_;
232 :    
233 :     my $pattern = $cgi->param('pattern');
234 :     my $dir = $cgi->param('dir');
235 :     my $file = "$dir/assigned_functions";
236 :     my @hits = &process_index($file,$pattern);
237 :     &format_function_table($cgi,$html,\@hits);
238 :     }
239 :    
240 :     sub process_subsystems_search {
241 :     my($cgi,$html) = @_;
242 :    
243 :     my $pattern = $cgi->param('pattern');
244 :     my $dir = $cgi->param('dir');
245 :     my $file = "$dir/Subsystems/subsystems.index";
246 :     my @hits = &process_index($file,$pattern);
247 :     &format_subsystems_table($cgi,$html,\@hits);
248 :     }
249 :    
250 :     sub process_peg2subsystems_search {
251 :     my($cgi,$html) = @_;
252 :    
253 :     my $peg = $cgi->param('peg');
254 :     my $dir = $cgi->param('dir');
255 :     my $file = "$dir/Subsystems/subsystems.index";
256 :     my %subs = map { $_->[1] => 1 } &process_index($file,$peg);
257 :     my @hits = sort { ($a->[0] cmp $b->[0]) or ($a->[1] cmp $b->[1]) or ($a->[2] cmp $b->[2])
258 :     or &SeedUtils::by_fig_id($a->[4],$b->[4]) }
259 :     map { chop; [split(/\t/,$_)] }
260 :     grep { ($_ =~ /^[^\t]*\t([^\t]+)/) && $subs{$1} }
261 :     `cat $file`;
262 :     &format_subsystems_table($cgi,$html,\@hits);
263 :     }
264 :    
265 :     sub comparison {
266 :     my($cgi,$html) = @_;
267 :    
268 :     my $refG = $cgi->param('reference');
269 :     my $queryG = $cgi->param('genome');
270 :     my $dir = $cgi->param('dir');
271 :     push(@$html,$cgi->start_form(-action => "sgv.cgi"),
272 :     '<b>Corresponding Features</b>: ',
273 :     $cgi->textfield(-name => 'pattern', -size => 30),
274 :     $cgi->submit('go'),
275 :     $cgi->hidden(-name => 'request', -value => 'corresponding', -override => 1),
276 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
277 :     $cgi->hidden(-name => 'reference', -value => $refG, -override => 1),
278 :     $cgi->end_form,
279 :     $cgi->hr,
280 :    
281 :     $cgi->start_form(-action => "sgv.cgi"),
282 :     '<b>Features in Query, but Not Reference</b>: ',
283 :     $cgi->textfield(-name => 'pattern', -size => 30),
284 :     $cgi->submit('go'),
285 :     $cgi->hidden(-name => 'request', -value => 'query_only', -override => 1),
286 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
287 :     $cgi->hidden(-name => 'reference', -value => $refG, -override => 1),
288 :     $cgi->end_form,
289 :     $cgi->hr,
290 :    
291 :     $cgi->start_form(-action => "sgv.cgi"),
292 :     '<b>Features in Reference, but Not Query</b>: ',
293 :     $cgi->textfield(-name => 'pattern', -size => 30),
294 :     $cgi->submit('go'),
295 :     $cgi->hidden(-name => 'request', -value => 'reference_only', -override => 1),
296 :     $cgi->hidden(-name => 'dir', -value => $dir, -override => 1),
297 :     $cgi->hidden(-name => 'reference', -value => $refG, -override => 1),
298 :     $cgi->end_form,
299 :     $cgi->hr
300 :     );
301 :     }
302 :    
303 :     sub process_corr_search {
304 :     my($cgi,$html) = @_;
305 :    
306 :     my $pattern = $cgi->param('pattern');
307 :     my $dir = $cgi->param('dir');
308 :     my $refG = $cgi->param('reference');
309 :     my $genome = $cgi->param('genome');
310 :     my $file = "$dir/CorrToReferenceGenomes/$refG";
311 :     my @corr = sort { &SeedUtils::by_fig_id($a->[0],$b->[0]) }
312 :     map { chomp; [split(/\t/,$_)] } `cat $file`;
313 :    
314 :     my @genes = map { ($_->[0],$_->[1]) } @corr;
315 :     my $col_headings = ['PEG','Function','reference','Reference Function','P-sc',
316 :     'BBH','Context-Matches','AliasesR'];
317 :     my $tab = [];
318 :     foreach my $entry (@corr)
319 :     {
320 :     my($pegG,$pegR,$n_context,undef,$funcG,$funcR,$aliasesG,$aliasesR,$bbh,undef,$psc) = @$entry;
321 :    
322 :     push(@$tab,[&peg_link($cgi,$pegG),$funcG,&ref_peg_link($cgi,$pegR),$funcR,
323 :     $psc,$bbh,$n_context,$aliasesR]);
324 :     }
325 :     my $filtered = &filter_tab_entries($tab,$pattern);
326 :     push(@$html,&SeedHTML::make_table($col_headings,$filtered,'Correspondences'));
327 :     }
328 :    
329 :     sub process_query_only_search {
330 :     my($cgi,$html) = @_;
331 :    
332 :     my $pattern = $cgi->param('pattern');
333 :     my $dir = $cgi->param('dir');
334 :     my $refG = $cgi->param('reference');
335 :     my $queryG = $cgi->param('genome');
336 :     my $file = "$dir/CorrToReferenceGenomes/$refG";
337 :     my %in_corr = map { $_ =~ /^(\S+)/; $1 => 1 } `cat $file`;
338 :    
339 :     my @to_show = grep { ! $in_corr{$_} }
340 :     map { $_ =~ /^(\S+)/; $1 }
341 :     `cut -f1 $dir/Features/peg/tbl`;
342 :     my %functionsH = map { $_ =~ /^(\S+)\t(\S.*\S)/; $1 => $2 } `cat $dir/assigned_functions`;
343 :     my $col_hdrs = ["PEG","Function"];
344 :     my @tab = map { [&peg_link($cgi,$_),$functionsH{$_} ? $functionsH{$_} : ""] }
345 :     sort { &SeedUtils::by_fig_id($a,$b) } @to_show;
346 :     my $filtered = &filter_tab_entries(\@tab,$pattern);
347 :     push(@$html,&SeedHTML::make_table($col_hdrs,$filtered,"Genes Missing in Reference Genome"));
348 :     }
349 :    
350 :     sub process_ref_only_search {
351 :     my($cgi,$html) = @_;
352 :    
353 :     my $pattern = $cgi->param('pattern');
354 :     my $dir = $cgi->param('dir');
355 :     my $refG = $cgi->param('reference');
356 :     my $queryG = $cgi->param('genome');
357 :     my $file = "$dir/CorrToReferenceGenomes/$refG";
358 :    
359 :     my %in_ref = map { $_ =~ /^\S+\t(\S+)/; $1 => 1 } `cat $file`;
360 :    
361 :     my $sapObject = SAPserver->new();
362 :     my $genomeH = $sapObject->all_features(-ids => [$refG], -type => 'peg');
363 :     my @to_show = grep { ! $in_ref{$_} } @{$genomeH->{$refG}};
364 :    
365 :     my $functionsH = $sapObject->ids_to_functions(-ids => \@to_show);
366 :     my $col_hdrs = ["PEG","Function"];
367 :     my @tab = map { [&ref_peg_link($cgi,$_),$functionsH->{$_} ? $functionsH->{$_} : ""] }
368 :     sort { &SeedUtils::by_fig_id($a,$b) } @to_show;
369 :     my $filtered = &filter_tab_entries(\@tab,$pattern);
370 :     push(@$html,&SeedHTML::make_table($col_hdrs,$filtered,"Genes Present Only in Reference Genome"));
371 :     }
372 :    
373 :     sub process_index {
374 :     my($file,$pattern) = @_;
375 :    
376 :     my @lines = `cat $file`;
377 :     if ( ! $pattern)
378 :     {
379 :     return map { chop; [split(/\t/,$_)] } @lines;
380 :     }
381 :     elsif ($pattern =~ /^\s*\/(.*)\/\s*$/)
382 :     {
383 :     return &perl_patmatch(\@lines,$1);
384 :     }
385 :     else
386 :     {
387 :     return &substr_match(\@lines,$pattern);
388 :     }
389 :     }
390 :    
391 :     sub perl_patmatch {
392 :     my($lines,$pat) = @_;
393 :    
394 :     my @lines = grep { $_ =~ /$pat/i } @$lines;
395 :     return map { chop; [split(/\t/,$_)] } @lines;
396 :     }
397 :    
398 :     sub substr_match {
399 :     my($lines,$pat) = @_;
400 :    
401 :     $pat =~ s/^\s+//;
402 :     $pat =~ s/\s+$//;
403 :     my @words = split(/\s+/,$pat);
404 :     my @lines = @$lines;
405 :     foreach my $word (@words)
406 :     {
407 :     @lines = grep { &matchword($word,$_) } @lines;
408 :     }
409 :     return map { chop; [split(/\t/,$_)] } @lines;
410 :     }
411 :    
412 :     sub matchword {
413 :     my($word,$str) = @_;
414 :    
415 :     my $wordL = lc $word;
416 :     my $strL = lc $str;
417 :     if (index($strL,$wordL) >= 0)
418 :     {
419 :     if ($wordL =~ /^fig\|\d+\.\d+\.peg\.\d+$/i)
420 :     {
421 :     my $wordQ = quotemeta $wordL;
422 :     return ($strL =~ /$wordQ\b/);
423 :     }
424 :     return 1;
425 :     }
426 :     return 0;
427 :     }
428 :    
429 :     sub format_function_table {
430 :     my($cgi,$html,$entries) = @_;
431 :    
432 :     my $col_hdrs = ['ID','Type','Function','Psi-blast','Subsystems'];
433 :     my $tab = [];
434 :    
435 :     foreach my $entry (@$entries)
436 :     {
437 :     my($fid,$function) = @$entry;
438 :     $fid =~ /fig\|\d+\.\d+\.([^\.]+)\.\d+$/;
439 :     my $type = $1;
440 :     if ($type eq "peg")
441 :     {
442 :     push(@$tab,[
443 :     &comp_reg_link($fid,$cgi),
444 :     'peg',
445 :     $function,
446 :     &psi_blast_link($fid),
447 :     &subsys_link($cgi,$fid)
448 :     ]);
449 :     }
450 :     else
451 :     {
452 :     push(@$tab,[$fid,$type,$function,"",""]);
453 :     }
454 :     }
455 :    
456 :     if (@$tab > 0)
457 :     {
458 :     push(@$html,&SeedHTML::make_table($col_hdrs,$tab,"Features"));
459 :     }
460 :     else
461 :     {
462 :     push(@$html,$cgi->h3('no matches'));
463 :     }
464 :     push(@$html,$cgi->hr,&query_link($cgi));
465 :     }
466 :    
467 :     sub format_subsystems_table {
468 :     my($cgi,$html,$entries) = @_;
469 :    
470 :     my $col_hdrs = ['Classification','Subsystem','Role','Variant','PEG'];
471 :     my $tab = [];
472 :    
473 :     foreach my $entry (@$entries)
474 :     {
475 :     my($class,$subsys,$role,$variant,$peg) = @$entry;
476 :     push(@$tab,[
477 :     $class,
478 :     $subsys,
479 :     $role,
480 :     $variant,
481 :     &peg_link($cgi,$peg)
482 :     ]);
483 :     }
484 :     if (@$tab > 0)
485 :     {
486 :     push(@$html,&SeedHTML::make_table($col_hdrs,$tab,"Subsystems"));
487 :     }
488 :     else
489 :     {
490 :     push(@$html,$cgi->h3('no matches'));
491 :     }
492 :     push(@$html,$cgi->hr,&query_link($cgi));
493 :     }
494 :    
495 :    
496 :     sub comp_reg_link {
497 :     my($peg,$cgi) = @_;
498 :    
499 :     my $target = "target.$$";
500 :     my $dir = $cgi->param('dir');
501 :     my $url = $cgi->url() . "?request=feature&fid=$peg&dir=$dir";
502 :     return "<a target=$target href=$url>$peg</a>";
503 :     }
504 :    
505 :     sub psi_blast_link {
506 :     my($peg) = @_;
507 :    
508 :     my $url = "http://seed-viewer.theseed.org/protein.cgi?prot=$peg&request=use_protein_tool&tool=Psi-Blast";
509 :     my $target = "target.$$";
510 :     return "<a target=$target href=$url>Psi</a>";
511 :     }
512 :    
513 :     sub ach_link {
514 :     my($cgi,$peg) = @_;
515 :     my $url = "http://seed-viewer.theseed.org/seedviewer.cgi?page=ACHresults&query=$peg";
516 :     my $target = "target.$$";
517 :     return "<a target=$target href=$url>ACH</a>";
518 :     }
519 :    
520 :     sub subsys_link {
521 :     my($cgi,$peg) = @_;
522 :    
523 :     my $dir = $cgi->param('dir');
524 :     my $url = $cgi->url() . "?request=peg2subsystems&dir=$dir&peg=$peg";
525 :     my $target = "target.$$";
526 :     return "<a target=$target href=$url>sub</a>";
527 :     }
528 :    
529 :     sub peg_link {
530 :     my($cgi,$peg) = @_;
531 :    
532 :     my $dir = $cgi->param('dir');
533 :     my $url = $cgi->url() . "?request=features&dir=$dir&pattern=$peg";
534 :     my $target = "target.$$";
535 :     return "<a target=$target href=$url>$peg</a>";
536 :     }
537 :    
538 :     sub ref_peg_link {
539 :     my($cgi,$peg) = @_;
540 :    
541 :    
542 :     my $target = "target.$$";
543 :     my $url = "http://seed-viewer.theseed.org/seedviewer.cgi?page=Annotation&feature=$peg";
544 :     return "<a target=$target href=$url>$peg</a>";
545 :     }
546 :    
547 :     sub query_link {
548 :     my($cgi) = @_;
549 :    
550 :     my $dir = $cgi->param('dir');
551 :     my $url = $cgi->url() . "?request=basic&dir=$dir";
552 :     return "<a href=$url>Basic Query Form</a>";
553 :     }
554 :    
555 :     sub filter_tab_entries {
556 :     my($tab,$pattern) = @_;
557 :    
558 :     if (! $pattern) { return $tab }
559 :    
560 :     my $filtered = [];
561 :     foreach my $entry (@$tab)
562 :     {
563 :     my @tmp = &substr_match([join("\t",@$entry)],$pattern);
564 :     if (@tmp > 0)
565 :     {
566 :     push(@$filtered,$entry);
567 :     }
568 :     }
569 :     return $filtered;
570 :     }
571 :    
572 :     sub process_id {
573 :     my($cgi,$html) = @_;
574 :    
575 :     my $id = $cgi->param('id');
576 :     if ($id =~ /^\s*(fig\|\d+\.\d+\.peg\.\d+)\s*$/)
577 :     {
578 :     my $peg = $1;
579 :     print $cgi->redirect("http://seed-viewer.theseed.org/seedviewer.cgi?page=Annotation&feature=$peg");
580 :     exit;
581 :     }
582 :     elsif ($id =~ /^\s*(\S+)\s*$/)
583 :     {
584 :     print $cgi->redirect("http://seed-viewer.theseed.org/seedviewer.cgi?page=ACHresults&query=$1");
585 :     exit;
586 :     }
587 :     else
588 :     {
589 :     push(@$html,$cgi->h2('Invalid request'));
590 :     }
591 :     }
592 :    
593 :     sub process_feature_display {
594 :     my($cgi,$html) = @_;
595 :    
596 :     my $dir = $cgi->param('dir');
597 :     my $seedV = SeedV->new($dir);
598 :     my $sapObject = SAPserver->new();
599 :    
600 :     my $fid = $cgi->param('fid');
601 :     my $func = $seedV->function_of($fid);
602 :     my $loc = $seedV->feature_location($fid);
603 :     my $dna_seq = $seedV->dna_seq($loc);
604 :    
605 :     my ($contig,$beg,$end);
606 :     if ($loc =~ /^(\S+)_(\d+)_(\d+)$/)
607 :     {
608 :     ($contig,$beg,$end) = ($1,$2,$3);
609 :     $loc = "contig: $1, from $2 to $3"
610 :     }
611 :     push(@$html,$cgi->h1("Feature: $fid"),$cgi->h2("Function: $func"),$cgi->h3("Location: $loc"));
612 :     &push_seq($cgi,$html,$fid,$dna_seq);
613 :    
614 :     if (($fid =~ /\.peg\.\d+$/) && $contig)
615 :     {
616 :     my $pseq = $seedV->get_translation($fid);
617 :     &push_seq($cgi,$html,$fid,$pseq);
618 :     &push_compare_regions($cgi,$html,$fid,$sapObject,$seedV,$contig,$beg,$end);
619 :     }
620 :     }
621 :    
622 :     sub push_seq {
623 :     my($cgi,$html,$id,$pseq) = @_;
624 :    
625 :     push(@$html,"<pre>\n>$id\n");
626 :     my $i;
627 :     for ($i=0; ($i < length($pseq)); $i += 60)
628 :     {
629 :     my $piece = ($i < (length($pseq) - 60)) ? substr($pseq,$i,60) : substr($pseq,$i);
630 :     push(@$html,"$piece\n");
631 :     }
632 :     push(@$html,"</pre>\n");
633 :     }
634 :    
635 :     sub push_compare_regions {
636 :     my($cgi,$html,$fid,$sapObject,$seedV,$contig,$beg,$end) = @_;
637 :    
638 :     my $min = ($beg < $end) ? ($beg - 4000) : $end - 4000;
639 :     my $max = ($beg < $end) ? ($end + 4000) : $beg + 4000;
640 :     my ($genes,$minV,$maxV) = $seedV->genes_in_region($contig,$min,$max);
641 :     my %genesG = map { ($_ => 1 ) } @$genes;
642 :     my %locsG = map { $_ => $seedV->feature_location($_) } @$genes;
643 :    
644 :     my $dir = $cgi->param('dir');
645 :     my $cache = "$dir/CorrToReferenceGenomes";
646 :     my %connected;
647 :     my %color; # note that for simplicity, I am using peg ids in the given genome as symbolic colors
648 :    
649 :     foreach $_ (`cat $cache/* | cut -f1,2`)
650 :     {
651 :     if (($_ =~ /^(\S+)\t(\S+)/) && $genesG{$1})
652 :     {
653 :     push(@{$connected{$1}},$2);
654 :     $color{$2} = $1;
655 :     }
656 :     }
657 :    
658 :     my $pinned = $connected{$fid};
659 :     my $locH = {};
660 :     if ($pinned)
661 :     {
662 :     my $pinLocH = $sapObject->fid_locations( -boundaries => 1, -ids => $pinned);
663 :     my @locations = map { &format_location($pinLocH->{$_}) } keys(%$pinLocH);
664 :     $locH = $sapObject->genes_in_region( -locations => \@locations, -includeLocation => 1);
665 :     }
666 :     # print &Dumper(\%genesG,\%locsG,$pinned,$locH,\%color);
667 :     push(@$html,@{&build_maps($seedV,$sapObject,$fid,\%genesG,\%locsG,$pinned,$locH,\%color)});
668 :     }
669 :    
670 :     sub format_location {
671 :     my($loc) = @_;
672 :    
673 :     if ($loc =~ /^(\d+\.\d+):(\S+)_(\d+)([+-])(\d+)$/)
674 :     {
675 :     my($genome,$contig,$beg,$strand,$ln) = ($1,$2,$3,$4,$5);
676 :    
677 :     my($min,$max);
678 :     if ($strand eq "+")
679 :     {
680 :     $min = &SeedUtils::max(1,$beg-4000);
681 :     $max = $beg+$ln+4000;
682 :     }
683 :     else
684 :     {
685 :     $min = &SeedUtils::max(1,$beg-($ln +4000));
686 :     $max = $beg + 4000;
687 :     }
688 :     return "$genome:$contig\_$min\_$max";
689 :     }
690 :     else
691 :     {
692 :     die "bad location: $loc";
693 :     }
694 :     }
695 :    
696 :     sub build_maps {
697 :     my($seedV,$sapObject,$pegG,$genesG,$locsG,$pinned,$locH,$color) = @_;
698 :    
699 :    
700 :     my @genome_ids = map { &SeedUtils::genome_of($_) } @$pinned;
701 :     my $genomeH = $sapObject->genome_names( -ids => \@genome_ids);
702 :     #
703 :     # first, compute a list of what we use to build each map. This will include
704 :     #
705 :     # [pinned_gene, Contig,Beg,End,GenusSpecies]]
706 :     # [[gene,Contig,Beg,End,color],...] sorted in order
707 :     #
708 :     my @map_data = ();
709 :     push(@map_data,&data_for_given_genome($pegG,$genesG,$locsG,$seedV));
710 :    
711 :     foreach my $pegR (@$pinned)
712 :     {
713 :     push(@map_data,&data_for_pinned($pegR,$locH,$color,$genomeH));
714 :     }
715 :     &set_colors($pegG,\@map_data);
716 :     my $functionH = &function_hash($sapObject,$seedV,\@map_data);
717 :    
718 :     my $gg = [];
719 : overbeek 1.2 my $sz_region = 12000;
720 : overbeek 1.1
721 :     foreach my $map_set (@map_data)
722 :     {
723 :     my($pin_data,$gene_data) = @$map_set;
724 :     my($peg,$contig,$beg,$end,$genus_species) = @$pin_data;
725 :    
726 :     if ($contig && $beg && $end) {
727 :     my $mid = int(($beg + $end) / 2);
728 :     my $min = int($mid - ($sz_region / 2));
729 :     my $max = int($mid + ($sz_region / 2));
730 :     my $genes = [];
731 :     foreach my $entry (@$gene_data)
732 :     {
733 :     my($fid1,$contig1,$beg1,$end1,$color) = @$entry;
734 :     my $beg1 = &in_bounds($min,$max,$beg1);
735 : overbeek 1.2 my $end1 = &in_bounds($min,$max,$end1);
736 : overbeek 1.1 my $function = $functionH->{$fid1};
737 :    
738 :     my $info = join('<br/>', "<b>PEG:</b> $fid1",
739 :     "<b>Contig:</b> $contig1",
740 :     "<b>Begin:</b> $beg1",
741 :     "<b>End:</b> $end1",
742 :     $function ? "<b>Function:</b> $function" : ()
743 :     );
744 :    
745 :     my $shape = "Rectangle";
746 :     if (($fid1 !~ /\.bs\./) && ($beg1 < $end1)) { $shape = "rightArrow" }
747 :     elsif (($fid1 !~ /\.bs\./) && ($beg1 > $end1)) { $shape = "leftArrow" }
748 :    
749 :     push(@$genes,[&min($beg1,$end1),
750 :     &max($beg1,$end1),
751 :     $shape,
752 :     ($fid1 !~ /\.bs\./) ? $color : 'black',
753 :     "",
754 :     $fid1,
755 :     $info
756 :     ]
757 :     );
758 :     }
759 :    
760 :     # Sequence title can be replaced by [ title, url, popup_text, menu, popup_title ]
761 :    
762 :     my $desc = "Genome: $genus_species<br />Contig: $contig";
763 :     my $map = [ [ SeedUtils::abbrev( $genus_species ), undef, $desc, undef, 'Contig' ],
764 :     0,
765 :     $max+1 - $min,
766 :     ($beg < $end) ? &decr_coords($genes,$min) : &flip_map($genes,$min,$max)
767 :     ];
768 :    
769 :     push(@$gg,$map);
770 :     }
771 :     }
772 :    
773 :     &GenoGraphics::disambiguate_maps($gg);
774 :     return &GenoGraphics::render($gg,700,4,0,2);
775 :     }
776 :    
777 :     sub data_for_given_genome {
778 :     my($peg,$pegs,$locs,$seedV) = @_;
779 :    
780 :     my @gene_data = ();
781 :     foreach my $peg1 (keys(%$pegs))
782 :     {
783 :     push(@gene_data,[$peg1,&split_loc($locs->{$peg1}),$peg1]);
784 :     }
785 :     my @gene_data = sort { ($a->[2]+$a->[3]) <=> ($b->[2]+$b->[3]) } @gene_data;
786 :     return [[$peg,&split_loc($locs->{$peg}),$seedV->genus_species],[@gene_data]];
787 :     }
788 :    
789 :     sub data_for_pinned {
790 :     my($peg,$locH,$color,$genomeH) = @_;
791 :    
792 :     my $genome = &SeedUtils::genome_of($peg);
793 :     my @tmp = grep { $locH->{$_}->{$peg} } keys(%$locH);
794 :     my $locH1 = $locH->{$tmp[0]};
795 :     my $pinned_data = [$peg,&split_new_loc($locH1->{$peg}->[0]),$genomeH->{&SeedUtils::genome_of($peg)}];
796 :     my @gene_data = ();
797 :     foreach my $peg1 (keys(%$locH1))
798 :     {
799 :     # print STDERR &Dumper($peg1,$color->{$peg1}); die "aborted";
800 :     push(@gene_data,[$peg1,&split_new_loc($locH1->{$peg1}->[0]),$color->{$peg1}]);
801 :     }
802 :     my @gene_data = sort { ($a->[2]+$a->[3]) <=> ($b->[2]+$b->[3]) } @gene_data;
803 :     return [$pinned_data,[@gene_data]];
804 :     }
805 :    
806 :     sub split_loc {
807 :     my($loc) = @_;
808 :    
809 :     if ($loc && ($loc =~ /^(.*:)?(\S+)_(\d+)_(\d+)$/))
810 :     {
811 :     return ($2,$3,$4);
812 :     }
813 :     die "bad_loc: $loc";
814 :     }
815 :    
816 :     sub split_new_loc {
817 :     my($loc) = @_;
818 :    
819 :     if ($loc =~ /^(\d+\.\d+):(\S+)_(\d+)([+-])(\d+)$/)
820 :     {
821 :     my($genome,$contig,$beg,$strand,$ln) = ($1,$2,$3,$4,$5);
822 :     if ($strand eq "+")
823 :     {
824 :     return ($contig,$beg,$beg+$ln-1);
825 :     }
826 :     else
827 :     {
828 :     return ($contig,$beg,$beg-($ln-1));
829 :     }
830 :     }
831 :     die "bad_loc: $loc";
832 :     }
833 :    
834 :     sub in_bounds {
835 :     my($min,$max,$x) = @_;
836 :    
837 :     if ($x < $min) { return $min }
838 :     elsif ($x > $max) { return $max }
839 :     else { return $x }
840 :     }
841 :    
842 :     sub decr_coords {
843 :     my($genes,$min) = @_;
844 :     my($gene);
845 :    
846 :     foreach $gene (@$genes) {
847 :     $gene->[0] -= $min;
848 :     $gene->[1] -= $min;
849 :     }
850 :     return $genes;
851 :     }
852 :    
853 :     sub function_hash {
854 :     my($sapObject,$seedV,$map_data) = @_;
855 :    
856 :     my $functionH = {};
857 :     my $gene_data = $map_data->[0]->[1];
858 :     foreach my $tuple (@$gene_data)
859 :     {
860 :     my $fid = $tuple->[0];
861 :     my $func = $seedV->function_of($fid);
862 :     $functionH->{$fid} = $func;
863 :     }
864 :    
865 :     my $i;
866 :     my @ids = ();
867 :     for ($i=1; ($i < @$map_data); $i++)
868 :     {
869 :     $gene_data = $map_data->[$i]->[1];
870 :     push(@ids,map { $_->[0] } @$gene_data);
871 :     }
872 :    
873 :     my $fH = $sapObject->ids_to_functions( -ids => \@ids);
874 :     while (my($id,$func) = each(%$fH))
875 :     {
876 :     $functionH->{$id} = $func;
877 :     }
878 :     return $functionH;
879 :     }
880 :    
881 :     sub flip_map {
882 :     my($genes,$min,$max) = @_;
883 :     my($gene);
884 :    
885 :     foreach $gene (@$genes) {
886 :     ($gene->[0],$gene->[1]) = ($max - $gene->[1],$max - $gene->[0]);
887 :     if ($gene->[2] eq "rightArrow") { $gene->[2] = "leftArrow" }
888 :     elsif ($gene->[2] eq "leftArrow") { $gene->[2] = "rightArrow" }
889 :     }
890 :     return $genes;
891 :     }
892 :    
893 :     sub set_colors {
894 :     my($red_peg,$map_data) = @_;
895 :    
896 :     my %colors;
897 :     foreach my $map (@$map_data)
898 :     {
899 :     my $genes = $map->[1];
900 :     foreach $_ (@$genes)
901 :     {
902 :     if ($_->[4])
903 :     {
904 :     $colors{$_->[4]}++;
905 :     }
906 :     }
907 :     }
908 :     my @by_occ = sort { $colors{$b} <=> $colors{$a} } keys(%colors);
909 :     my $i;
910 :     my %to_color;
911 :     for ($i=1; ($i <= @by_occ); $i++)
912 :     {
913 :     $to_color{$by_occ[$i-1]} = "color$i";
914 :     }
915 :     $to_color{$red_peg} = "red";
916 :     foreach my $map (@$map_data)
917 :     {
918 :     my $genes = $map->[1];
919 :     foreach $_ (@$genes)
920 :     {
921 :     if ($_->[4])
922 :     {
923 :     $_->[4] = $to_color{$_->[4]};
924 :     }
925 :     else
926 :     {
927 :     $_->[4] = 'grey';
928 :     }
929 :     }
930 :     }
931 :     }
932 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3