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

Annotation of /FigWebServices/families_on_tree.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (view) (download)

1 : overbeek 1.1 ########################################################################
2 :     use CGI;
3 :    
4 :    
5 :     if (-f "$FIG_Config::data/Global/why_down")
6 :     {
7 :     local $/;
8 :     open my $fh, "<$FIG_Config::data/Global/why_down";
9 :     my $down_msg = <$fh>;
10 :    
11 :     print CGI::header();
12 :     print CGI::head(CGI::title("SEED Server down"));
13 :     print CGI::start_body();
14 :     print CGI::h1("SEED Server down");
15 :     print CGI::p("The seed server is not currently running:");
16 :     print CGI::pre($down_msg);
17 :     print CGI::end_body();
18 :     exit;
19 :     }
20 :    
21 :     if ($FIG_Config::readonly)
22 :     {
23 :     CGI::param("user", undef);
24 :     }
25 :     ########################################################################
26 :     use CGI;
27 :    
28 :    
29 :     if (-f "$FIG_Config::data/Global/why_down")
30 :     {
31 :     local $/;
32 :     open my $fh, "<$FIG_Config::data/Global/why_down";
33 :     my $down_msg = <$fh>;
34 :    
35 :     print CGI::header();
36 :     print CGI::head(CGI::title("SEED Server down"));
37 :     print CGI::start_body();
38 :     print CGI::h1("SEED Server down");
39 :     print CGI::p("The seed server is not currently running:");
40 :     print CGI::pre($down_msg);
41 :     print CGI::end_body();
42 :     exit;
43 :     }
44 :    
45 :     if ($FIG_Config::readonly)
46 :     {
47 :     CGI::param("user", undef);
48 :     }
49 :     ########################################################################
50 :     # -*- perl -*-
51 :     #
52 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
53 :     # for Interpretations of Genomes. All Rights Reserved.
54 :     #
55 :     # This file is part of the SEED Toolkit.
56 :     #
57 :     # The SEED Toolkit is free software. You can redistribute
58 :     # it and/or modify it under the terms of the SEED Toolkit
59 :     # Public License.
60 :     #
61 :     # You should have received a copy of the SEED Toolkit Public License
62 :     # along with this program; if not write to the University of Chicago
63 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
64 :     # Genomes at veronika@thefig.info or download a copy from
65 :     # http://www.theseed.org/LICENSE.TXT.
66 :     #
67 :    
68 :     use DB_File;
69 :     use URI::Escape; # uri_escape
70 :     use gjoseqlib;
71 :     use HTML;
72 :     use strict;
73 :     use CGI;
74 :     my $cgi = new CGI;
75 :     use SeedEnv;
76 :     use tree_utilities;
77 :    
78 :     if (0) {
79 :     print $cgi->header;
80 :     my @params = $cgi->param;
81 :     print "<pre>\n";
82 :     foreach $_ (@params) {
83 :     print "$_\t:",join(",",$cgi->param($_)),":\n";
84 :     }
85 :     exit;
86 :     }
87 :    
88 :     my $html = [];
89 :    
90 :     my $node = $cgi->param('node');
91 :     my $family = $cgi->param('family');
92 :     my $ali = $cgi->param('alignment');
93 :     my $tree = $cgi->param('tree');
94 :     my $request = $cgi->param('request');
95 :     my $keywords = $cgi->param('keywords');
96 :     my $function = $cgi->param('function');
97 :     my $dataD = $cgi->param('dataD');
98 :     my $dataDir = "/homes/overbeek/Ross/MakeCS/Data/CS";
99 :     my $dataDF = "$dataDir/$dataD";
100 :    
101 :     if ($request eq "show_otus")
102 :     {
103 :     &show_otus($cgi,$dataDir);
104 :     exit;
105 :     }
106 :     elsif (($request eq "show_options_for_otu") && $dataD)
107 :     {
108 :     &show_options_for_otu($cgi,$dataD);
109 :     exit;
110 :     }
111 :    
112 :     if (($request eq "show_virulence_functions") && (-s "$dataDF/virulence.functions"))
113 :     {
114 :     &show_virulence_functions($cgi,$dataDF,$html);
115 :     }
116 :     elsif (($request eq 'show_indexed_funcs') && $keywords)
117 :     {
118 :     &show_indexed_funcs($cgi,$dataDF,$html,$keywords);
119 :     }
120 :     elsif (($request eq "show_func") && $function)
121 :     {
122 :     $function =~ s/^\s+//;
123 :     $function =~ s/\s+$//;
124 :     &show_func($cgi,$dataDF,$html,$function);
125 :     }
126 :     elsif (($request eq "show_ali_or_occurs_tree") && $ali)
127 :     {
128 :     &show_ali($cgi);
129 :     exit;
130 :     }
131 :     elsif (($request eq "show_ali_or_occurs_tree") && $tree)
132 :     {
133 :     &show_occurs_tree($cgi,$dataDF);
134 :     exit;
135 :     }
136 :     elsif (($request eq "show_family_pegs") && $family)
137 :     {
138 :     &show_family1($cgi,$dataDF,$html,$family);
139 :     }
140 :     elsif (($request eq "show_family_tree") && $family)
141 :     {
142 :     &show_family2($cgi,$dataDF,$html,$family);
143 :     }
144 :     elsif (($request eq "show_node") && $node)
145 :     {
146 :     &show_changes($cgi,$dataDF,$html,$node);
147 :     }
148 :     elsif ($request eq "show_otu_tree")
149 :     {
150 :     &show_otu_tree($cgi,$dataDF,$html,'families');
151 :     }
152 :     elsif ($request eq "show_adjacency")
153 :     {
154 :     &show_otu_tree($cgi,$dataDF,$html,'adjacency');
155 :     }
156 :     elsif ($request eq "show_clusters")
157 :     {
158 :     &show_clusters($cgi,$dataDF,$html);
159 :     }
160 :     elsif ($request eq "show_kovbassa")
161 :     {
162 :     &show_kovbassa_signatures($cgi,$dataDF,$html);
163 :     }
164 :     elsif ($request eq "compute_kovbassa_sigs")
165 :     {
166 :     &compute_kovbassa_signatures($cgi,$dataDF,$html);
167 :     }
168 :     else
169 :     {
170 :     push(@$html,"<h1>Invalid request</h1>");
171 :     }
172 :     &HTML::show_page($cgi,$html);
173 :     exit;
174 :    
175 :     sub show_otu_tree {
176 :     my($cgi,$dataDF,$html,$type) = @_;
177 :    
178 :     $dataDF =~ /([^\/]+)$/;
179 :     my $dataD = $1;
180 :     my @tree;
181 :     my @tmp = `cat $dataDF/readable.tree`;
182 :     foreach $_ (@tmp)
183 :     {
184 :     if ($_ =~ /^(.*\+ )(n\d+)$/)
185 :     {
186 :     my $start = $1;
187 :     my $node = $2;
188 :     my $link = &show_node_link($dataD,$node,$type);
189 :     push(@tree,"$start$link\n");
190 :     }
191 :     elsif ($_ =~ /^(.*\- )(\d+\.\d+)(:.*)$/)
192 :     {
193 :     my $start = $1;
194 :     my $node = $2;
195 :     my $end = $3;
196 :     push(@tree,"$start",&show_node_link($dataD,$node,$type),"$end\n");
197 :     }
198 :     else
199 :     {
200 :     push(@tree,$_);
201 :     }
202 :     }
203 :    
204 :     my $tree = join("",@tree);
205 :     push(@$html,"<pre>\n$tree\n</pre>\n");
206 :     }
207 :    
208 : overbeek 1.4 sub fam_in_set {
209 :     my($fam,$set) = @_;
210 :    
211 :     my @fams = split(/,/,$set);
212 :     my $i;
213 :     for ($i=0; ($i < @fams) && ($fam ne $fams[$i]); $i++) {}
214 :     return ($i < @fams);
215 :     }
216 :    
217 : overbeek 1.1 sub show_family1 {
218 : overbeek 1.4 my($cgi,$dataDF,$html,$families) = @_;
219 : overbeek 1.1
220 :     $dataDF =~ /([^\/]+)$/;
221 :     my $dataD = $1;
222 :    
223 :     my $func;
224 :     foreach $_ (`cut -f1,2 $dataDF/families.all`)
225 :     {
226 : overbeek 1.4 # if (($_ =~ /^(\d+)\t(\S[^\t]*\S)/) && ($1 == $families)) { $func = $2 }
227 :     if (($_ =~ /^(\d+)\t(\S[^\t]*\S)/) && &fam_in_set($1,$families)) { $func = $2 }
228 : overbeek 1.1 }
229 :     my %genome_names = map { ($_ =~ /^(\S+)\t(\S.*\S)/) ? ($1 => $2) : () } `cat $dataDF/genome.names`;
230 :     my $col_hdrs = ['','Genome','Genome Name','Peg'];
231 :     my @tuples = map { my $tuple = $_; $tuple->[3] = &peg_link($tuple->[3]); $tuple }
232 :     sort { $a->[1] cmp $b->[1] }
233 : overbeek 1.4 map { (($_ =~ /^(\d+)\t([^\t]*)\t(fig\|(\d+\.\d+)\.peg\.\d+)/) && &fam_in_set($1,$families)) ?
234 : overbeek 1.1 [$cgi->checkbox(-name => "check.$3",-checked => 1,-label => ''),
235 :     $4,
236 :     $genome_names{$4},
237 :     $3] : () }
238 :     `cut -f1,2,4 $dataDF/families.all`;
239 : overbeek 1.4 push(@$html,$cgi->start_form(-action => "./families5_on_tree.cgi"));
240 :     push(@$html,"<br><br>\n",&HTML::make_table($col_hdrs,\@tuples,"Distribution of Family $families: $func"),$cgi->hr,"\n");
241 : overbeek 1.1 push(@$html,"<input type=hidden name=dataD value=$dataD>\n");
242 :     push(@$html,"<input type=hidden name=request value=show_ali_or_occurs_tree>\n");
243 :     push(@$html,$cgi->submit('alignment'),
244 :     $cgi->checkbox(-name => "dna",-checked => 0,-label => 'DNA'),
245 :     "<br><br>",
246 :     $cgi->submit('tree'),
247 :     $cgi->end_form());
248 :     }
249 :    
250 :     sub show_family2 {
251 :     my($cgi,$dataDF,$html,$family) = @_;
252 :     # if union is specified, you want to treat all families with the same function
253 :     # as "present"
254 :     my $union = $cgi->param('union');
255 :     my $func; # used only if union is requested
256 :     my @fam_func_peg_tuples = map { chomp; [split(/\t/,$_)] } `cut -f1,2,4 $dataDF/families.all`;
257 :     if ($union) {
258 :     my @tmp = grep { $_->[0] == $family } @fam_func_peg_tuples;
259 :     $func = $tmp[0]->[1];
260 :     }
261 :    
262 :     $dataDF =~ /([^\/]+)$/;
263 :     my $dataD = $1;
264 :    
265 :     push(@$html,"<br><hr>",&show_fam_table_link($dataDF,$family),"<br><br>");
266 :     my @tuples = grep { $union ? ($_->[1] eq $func) : ($_->[0] == $family) }
267 :     @fam_func_peg_tuples;
268 :     my $func = (@tuples > 0) ? $tuples[0]->[1] : 'hypothetical protein';
269 :     my %has = map { ($_->[2] =~ /fig\|(\d+\.\d+)/) ? ($1 => 1) : () } @tuples;
270 :     my @tree;
271 :     my %node_vals;
272 :     if (! -s "$dataDF/fams.on.tree.db")
273 :     {
274 :     %node_vals = map { (($_ =~ /^(\S+)\t(n\d+)\t(\S+)/) && ($family eq $1)) ? ($2 => $3) : () } `cat $dataDF/families.on.tree`;
275 :     }
276 :     else
277 :     {
278 :     if ($dataD ne "bacteria") { die "only bacteria have a db"; }
279 :     my $bacteriaD = "/homes/overbeek/Ross/MakeCS/Data/CS/bacteria";
280 :     my $fams_on_tree_db = "$bacteriaD/fams.on.tree.db";
281 :     my %fams_on_tree;
282 :     tie %fams_on_tree,"DB_File", $fams_on_tree_db, O_RDWR, 0666, $DB_HASH
283 :     or die "Cannot open file '$fams_on_tree_db': $!\n";
284 :     my $node_vals = $fams_on_tree{$family};
285 :     %node_vals = map { $_ =~ s/://g; ($_ =~ /^(\S+),(\S+)/) ? ($1 => $2) : () } split("::",$node_vals);
286 :     }
287 :    
288 :     my @tmp = `cat $dataDF/readable.tree`;
289 :     foreach $_ (@tmp)
290 :     {
291 :     if ($_ =~ /^(.*\- )(\d+\.\d+)(:.*)$/)
292 :     {
293 :     my $start = $1;
294 :     my $genome = $2;
295 :     my $end = $3;
296 :     if ($has{$genome})
297 :     {
298 :     push(@tree,"$start<b>$genome$end</b>\n");
299 :     }
300 :     else
301 :     {
302 :     push(@tree,"$start$genome$end\n");
303 :     }
304 :     }
305 :     elsif (($_ =~ /^(.*\+ )(n\d+)$/) && (! $union))
306 :     {
307 :     my $start = $1;
308 :     my $node = $2;
309 :     my $status = $node_vals{$node};
310 :     my $to_show = $node . ":$status";
311 :     push(@tree,"$start$to_show\n");
312 :     }
313 :     else
314 :     {
315 :     push(@tree,$_);
316 :     }
317 :     }
318 :     my $tree = join("",@tree);
319 :     push(@$html,"<br><h2>$func</h2><br><pre>\n$tree\n</pre>\n");
320 :     }
321 :    
322 :     sub show_changes {
323 :     my($cgi,$dataDF,$html,$node) = @_;
324 :    
325 :     my $type = $cgi->param('type');
326 :     if ($type eq 'families')
327 :     {
328 :     &show_changes_families($cgi,$dataDF,$html,$node);
329 :     }
330 :     else
331 :     {
332 :     &show_changes_adjacency($cgi,$dataDF,$html,$node);
333 :     }
334 :     }
335 :    
336 :     sub show_changes_adjacency {
337 :     my($cgi,$dataDF,$html,$node) = @_;
338 :     $dataDF =~ /([^\/]+)$/;
339 :     my $dataD = $1;
340 :     my $col_hdrs = ['Family','Function','Ancestral','New','Compare'];
341 :     my @events = grep { $node eq $_->[1] }
342 :     map { ($_ =~ /(\S+)\t(\S+)\t(\d+):\S+\t(\d+):\S+\t(\d+)/) ? [$1,$2,$3,$4,$5] : () }
343 :     `cat $dataDF/placed.events`;
344 :     my %families = map { ($_->[2] => [$_->[3],$_->[4]])} @events;
345 :     my %pegs_needed = map { (($_->[2] => 1), ($_->[3] => 1),($_->[4] => 1)) } @events;
346 :     my %fam_peg = map { my $x;
347 :     (($_ =~ /^(\d+)\t\S+\t(\d+)\t\S+\t\S+\t(\S+)\t(\S+)/) &&
348 :     ($x = $families{$1}) && (($x->[0] eq $2) || ($x->[1] eq $2))) ? ("$1,$2" => $3) : () }
349 :     `cat $dataDF/adjacency.of.unique`;
350 :     my %peg_to_func = map { (($_ =~ /^([^t]+)\t([^\t]*)\t(\S+)/) && $pegs_needed{$1}) ? ($3 => $2) : () } `cut -f1,2,4 $dataDF/families.all`;
351 :     my @rows;
352 :     my $ancestor;
353 :     foreach my $event (@events)
354 :     {
355 :     my($anc,$node,$fam,$fam1,$fam2) = @$event;
356 :     $ancestor = $anc;
357 :     my $peg1 = $fam_peg{"$fam,$fam1"};
358 :     my $peg2 = $fam_peg{"$fam,$fam2"};
359 :     my $func = $peg_to_func{$peg1};
360 :     if ($peg1 && $peg2 && $func)
361 :     {
362 :     push(@rows,[&show_fam_table_link($dataDF,$fam),
363 :     $func,
364 :     &peg_link($peg1),
365 :     &peg_link($peg2),
366 :     &compare_link([$peg1,$peg2])]);
367 :     }
368 :     }
369 :     push(@$html,&HTML::make_table($col_hdrs,\@rows,"Changes in Adjacency from $ancestor"));
370 :     }
371 :    
372 :    
373 :     sub show_changes_families {
374 :     my($cgi,$dataDF,$html,$node) = @_;
375 :    
376 :     $dataDF =~ /([^\/]+)$/;
377 :     my $dataD = $1;
378 :    
379 :     my %func = map { ($_ =~ /^(\d+)\t(\S[^\t]*\S)/) ? ($1 => $2) : () } `cut -f1,2 $dataDF/families.all`;
380 :     my $col_hdrs = ['Show Where','Show PEGs','Family','Function','Clusters','Coupling'];
381 :     my @tmp = grep { (($_ =~ /^\S+\t(\S+)/) && ($1 eq $node)) }
382 :     `cat /$dataDF/where.shifts.occurred`;
383 :     my @tabG = sort { ($a->[4] cmp $b->[4]) or ($a->[3] <=> $b->[3]) }
384 :     map { ($_ =~ /^(\S+)\t\S+\t(\S+)\t0\t1/) ? [&show_fam_links($dataDF,$2),$1,$2,$func{$2}] : () }
385 :     @tmp;
386 :     # tabG entries are [linkT,linkP,ancestor,fam,func]
387 :    
388 :     # try to pick up the ancestor node from the first entry in @tabG
389 :     # If you cannot get it, try to take it from @tabL
390 :     my $anc = (@tabG > 0) ? $tabG[0]->[-3] : undef;
391 :     foreach $_ (@tabG) { splice(@$_,2,1) } ### get rid of ancestor
392 :     ## tabG entries are [linkT,linkP,fam,func]
393 :    
394 :     my @tabL = sort { ($a->[4] cmp $b->[4]) or ($a->[3] <=> $b->[3]) }
395 :     map { ($_ =~ /^(\S+)\t\S+\t(\S+)\t1\t0/) ? [&show_fam_links($dataDF,$2),$1,$2,$func{$2}] : () }
396 :     @tmp;
397 :     if (! $anc)
398 :     {
399 :     $anc = (@tabL > 0) ? $tabL[0]->[-3] : '';
400 :     }
401 :     foreach $_ (@tabL) { splice(@$_,2,1) } ### get rid of ancestor
402 :    
403 :     ## @tabG and @tabL are of the form [link-to-tree,link-to-peg-display,family,function]]
404 :     ## we now add coupling data.
405 :    
406 :     my $with_couplingL = &build_table(\@tabL,$dataDF);
407 :     my $with_couplingG = &build_table(\@tabG,$dataDF);
408 :    
409 :     push(@$html,&HTML::make_table($col_hdrs,$with_couplingG,"Families Gained from Ancestor $anc"),$cgi->hr,"\n");
410 :     push(@$html,&HTML::make_table($col_hdrs,$with_couplingL,"Families Lost from Ancestor $anc"),$cgi->hr,"\n");
411 :     }
412 :    
413 :     sub coupling_data {
414 :     my($dataDF,$famsH) = @_;
415 :    
416 :     my $coupled = {};
417 :     foreach $_ (`cat $dataDF/coupled.families`)
418 :     {
419 :     if (($_ =~ /(\S+)\t(\S+)\t(\d+)/) && $famsH->{$1} && $famsH->{$2} && ($1 ne $2))
420 :     {
421 :     push(@{$coupled->{$1}},$2);
422 :     }
423 :     }
424 :     return $coupled;
425 :     }
426 :    
427 :     sub cluster_link_and_cluster_html {
428 :     my($family,$coupledH,$fam_to_func,$dataD) = @_;
429 :    
430 :     my $cluster_link = '';
431 :     my $coupled = $coupledH->{$family};
432 :     my $coupled_html = "";
433 :     if ($coupled && (@$coupled > 0))
434 :     {
435 :     $cluster_link = &show_clusters_link($dataD,$family,$coupled);
436 :     $coupled_html = join("<br>",map { &show_fam_table_link($dataD,$_) . '(' .
437 :     &show_func_link($dataD,$fam_to_func->{$_}) . ')'
438 :     } @$coupled);
439 :     }
440 :     return ($cluster_link,$coupled_html);
441 :     }
442 :    
443 :     sub build_table {
444 :     my($tab,$dataDF) = @_;
445 :     $dataDF =~ /([^\/]+)$/;
446 :     my $dataD = $1;
447 :    
448 :     my %famH = map { ($_->[-2] => 1) } @$tab;
449 :     my %fam_to_func = map { ($_->[2] => $_->[3]) } @$tab;
450 :     my $coupledH = &coupling_data($dataDF,\%famH);
451 :     my @with_coupling;
452 :     foreach my $tuple (@$tab)
453 :     {
454 :     my($link1,$link2,$family,$function) = @$tuple;
455 :     $tuple->[3] = &show_func_link($dataD,$function);
456 :     my($cluster_link,$coupled_html) = &cluster_link_and_cluster_html($family,$coupledH,\%fam_to_func,$dataD);
457 :     $tuple->[4] = $cluster_link;
458 :     $tuple->[5] = $coupled_html;
459 :     push(@with_coupling,$tuple);
460 :     }
461 :     return \@with_coupling;
462 :     }
463 :    
464 :     sub show_fam_links {
465 :     my($dataDF,$fam) = @_;
466 :    
467 :     $dataDF =~ /([^\/]+)$/;
468 :     my $dataD = $1;
469 :    
470 : overbeek 1.2 return ("<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_family_tree&family=$fam&dataD=$dataD>tree</a>",
471 :     "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_family_pegs&family=$fam&pegs=1&dataD=$dataD>PEGs:$fam</a>");
472 : overbeek 1.1 }
473 :    
474 :     sub show_fam_table_link {
475 :     my($dataDF,$fam) = @_;
476 :    
477 :     $dataDF =~ /([^\/]+)$/;
478 :     my $dataD = $1;
479 :    
480 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?family=$fam&pegs=1&request=show_family_pegs&dataD=$dataD>$fam: PEGs</a>";
481 : overbeek 1.1 }
482 :    
483 :     sub show_func {
484 :     my($cgi,$dataDF,$html,$function) = @_;
485 :     $dataDF =~ /([^\/]+)$/;
486 :     my $dataD = $1;
487 :    
488 :     my %counts;
489 :     my %fams = map { my $x;
490 :     if (($_ =~ /^(\d+)\t(\S[^\t]*\S)/) && ($2 eq $function))
491 :     {
492 :     $counts{$1}++;
493 :     ($1 => 1);
494 :     }
495 :     else
496 :     {
497 :     ()
498 :     }
499 :     } `cat $dataDF/families.all`;
500 :     my $col_hdrs = ['Family','PEGs in Family','Distribution on Tree','Size'];
501 :     my @tab = map { [$_,
502 :     &show_fam_table_link($dataDF,$_),
503 :     &show_fam_on_tree_link($dataDF,$_),
504 :     $counts{$_}
505 :     ] } sort { $a <=> $b } keys(%fams);
506 :     push(@$html,"<h1>$function</h1>", &HTML::make_table($col_hdrs,\@tab,"Families with Function"),"<hr>");
507 :    
508 :     my @fams = keys(%fams);
509 :     if (@fams > 1)
510 :     {
511 :     push(@$html,"<br>",&show_fam_on_tree_link($dataDF,$fams[0],"union")," with union of families<br><br>");
512 :     }
513 :     push(@$html,"<hr>");
514 :    
515 :     my @gained = map { (($_ =~ /^(\S+)\t(\S+)\t(\d+)\t0\t1/) && $fams{$3}) ?
516 :     [$1,&show_node_link($dataD,$2,'families'),&show_fam_table_link($dataDF,$3)] : () }
517 :     `cat $dataDF/where.shifts.occurred`;
518 :     push(@$html,&HTML::make_table(['Ancestor','Descendant','Family'],\@gained,"Where $function was GAINED"));
519 :     push(@$html,"<hr><br>");
520 :     my @lost = map { (($_ =~ /^(\S+)\t(\S+)\t(\d+)\t1\t0/) && $fams{$3}) ?
521 :     [$1,&show_node_link($dataD,$2,'families'),&show_fam_table_link($dataDF,$3)] : () }
522 :     `cat $dataDF/where.shifts.occurred`;
523 :     push(@$html,&HTML::make_table(['Ancestor','Descendant','Family'],\@lost,"Where $function was LOST"));
524 :     }
525 :    
526 :     sub show_indexed_funcs {
527 :     my($cgi,$dataDF,$html,$keywords) = @_;
528 :    
529 :     $dataDF =~ /([^\/]+)$/;
530 :     my $dataD = $1;
531 :    
532 :     my $functions_in_fams = &functions_in_at_least_one_family($dataDF);
533 :     # $keywords = "$dataD " . $keywords; ### tell the user to add it,if necessary
534 :    
535 :     my %funcs_to_show;
536 :     foreach my $func (`svr_sphinx_indexing -k \'$keywords\' | cut -f1 | svr_function_of | cut -f2`)
537 :     {
538 :     chomp $func;
539 : overbeek 1.3 $func =~ s/\s*\#.*$//;
540 : overbeek 1.1 if ($functions_in_fams->{$func})
541 :     {
542 :     $funcs_to_show{$func}++;
543 :     }
544 :     }
545 :     my @funcs = sort { $funcs_to_show{$b} <=> $funcs_to_show{$a} } keys(%funcs_to_show);
546 :     if (@funcs == 0)
547 :     {
548 :     push(@$html,"<h1>Sorry, no functions matched</h1>\n");
549 :     }
550 :     else
551 :     {
552 :     my @links = map { [&show_func_link($dataD,$_)] } @funcs;
553 :     push(@$html,&HTML::make_table(['Possible Functions'],\@links,"Possible functions - Select to find nodes where shifts occurred"));
554 :     }
555 :     }
556 :    
557 :     sub show_virulence_functions {
558 :     my($cgi,$dataDF,$html) = @_;
559 :    
560 :     $dataDF =~ /([^\/]+)$/;
561 :     my $dataD = $1;
562 :    
563 :     my $functions_in_fams = &functions_in_at_least_one_family($dataDF);
564 :     my @virulence_functions = map { chomp; $functions_in_fams->{$_} ? $_ : () } `cat $dataDF/virulence.functions`;
565 :     my @links = map { [&show_func_link($dataD,$_)] } sort @virulence_functions;
566 :     push(@$html,&HTML::make_table(['Function Sometimes Associated with Virulence'],
567 :     \@links,
568 :     'Functions Known to Be Associated with Virulence in Some Organisms'));
569 :     }
570 :     sub functions_in_at_least_one_family {
571 :     my($dataDF) = @_;
572 :    
573 :     my %functions_in_fams = map { chomp; ($_ => 1) } `cut -f2 $dataDF/families.all`;
574 :     return \%functions_in_fams;
575 :     }
576 :    
577 :     sub subsystems_of {
578 :     my($fig,$reaction_to_roles,$reaction) = @_;
579 :    
580 :     my $roles_and_pegs = $reaction_to_roles->{$reaction};
581 :     my @roles = map { my $x = $_; $x =~ s/^[^:]+://; $x } @{$reaction_to_roles->{$reaction}};
582 :     my $printable_roles = "";
583 :     if (@roles > 0)
584 :     {
585 :     my %tmp = map { $_ =~ /^([^:]+):(\S.*\S)/; ($2 => $1) } @$roles_and_pegs;
586 :     my @tmp = map { &peg_link($tmp{$_}) . "<br>" . $_ } sort keys(%tmp);
587 :     $printable_roles = join(",",@tmp);
588 :     }
589 :     my %subsys;
590 :     foreach my $role (@roles)
591 :     {
592 :     foreach my $s ($fig->role_to_subsystems($role))
593 :     {
594 :     $subsys{$s} = 1;
595 :     }
596 :     }
597 :     return join("\t",sort keys(%subsys)) . "<br>" . $printable_roles;
598 :     }
599 :    
600 :     sub show_ali {
601 :     my($cgi) = @_;
602 :     my @checked = map { ($_ =~ s/^check.//) ? $_ : () } $cgi->param();
603 :     if (@checked > 1)
604 :     {
605 :     my %checked = map { $_ => 1 } @checked;
606 :     my $dataD = $cgi->param('dataD');
607 :     my $dataDir = "/homes/overbeek/Ross/MakeCS/Data/CS";
608 :     my $dataDF = "$dataDir/$dataD";
609 :    
610 :     my $tmp_seqs_in = "$FIG_Config::temp/tmp$$.seqs.in";
611 :     my $tmp_seqs_ali = "$FIG_Config::temp/tmp$$.ali";
612 :     open(SEQS,">$tmp_seqs_in") || die "could not open $tmp_seqs_in";
613 :     my $dna = $cgi->param('dna');
614 :     opendir(GENOMES,"$dataDF/Seqs") || die "could not open $dataDF/Seqs";
615 :     my @genomes = grep { $_ !~ /^\./ } readdir(GENOMES);
616 :     closedir(GENOMES);
617 :     foreach my $g (@genomes)
618 :     {
619 :     my @seqs;
620 :     if($dna)
621 :     {
622 :     @seqs = grep { $checked{$_->[0]} } &gjoseqlib::read_fasta("$dataDF/PegDNA/$g");
623 :     }
624 :     else
625 :     {
626 :     @seqs = grep { $checked{$_->[0]} } &gjoseqlib::read_fasta("$dataDF/Seqs/$g");
627 :     }
628 :     foreach my $tuple(@seqs)
629 :     {
630 :     my($peg,undef,$seq) = @$tuple;
631 :     print SEQS ">$peg\n$seq\n";
632 :     }
633 :     }
634 :     close(SEQS);
635 :    
636 :     &SeedUtils::run("svr_align_seqs < $tmp_seqs_in | svr_ali_to_html > $tmp_seqs_ali");
637 :     open(SEQS,"<$tmp_seqs_ali") || die "could not open $tmp_seqs_ali";
638 :     print $cgi->header;
639 :     while (defined($_ = <SEQS>))
640 :     {
641 :     $_ =~ s/\b(fig\|\d+\.\d+\.peg\.\d+)/<a target=_blank href=http:\/\/pubseed.theseed.org\/seedviewer.cgi?page=Annotation&feature=$1>$1<\/a>/g;
642 :     print $_;
643 :     }
644 :     close(SEQS);
645 :     unlink($tmp_seqs_in,$tmp_seqs_ali);
646 :     }
647 :     }
648 :    
649 :     sub show_occurs_tree {
650 :     my($cgi,$dataDF) = @_;
651 :    
652 :     my @checked = map { ($_ =~ s/^check.//) ? $_ : () } $cgi->param();
653 :     if (@checked > 3)
654 :     {
655 :     my %genome_name =
656 :     map { ($_ =~ /^(\d+\.\d+)\t(\S.*\S)/) ? ($1 => $2) : () }
657 :     `cat $dataDF/genome.names`;
658 :     my $tmp_labels = "$FIG_Config::temp/tmp$$.labels";
659 :     open(LABELS,">$tmp_labels") || die "could not open $tmp_labels";
660 :     foreach my $peg (@checked)
661 :     {
662 :     if ($peg =~ /^fig\|(\d+\.\d+)/)
663 :     {
664 :     my $g = $genome_name{$1};
665 :     my $link = &peg_link($peg);
666 :     print LABELS "$peg\t$link: $g\n";
667 :     }
668 :     }
669 :     close(LABELS);
670 :     my $tmp_seqs = "$FIG_Config::temp/tmp$$.seqs";
671 :     open(SEQS,"| svr_fasta -fasta -protein | svr_align_seqs | svr_tree | sketch_tree -m -l $tmp_labels > $tmp_seqs")
672 :     || die "could not open $tmp_seqs";
673 :     foreach $_ (@checked)
674 :     {
675 :     print SEQS "$_\n";
676 :     }
677 :     close(SEQS);
678 :     open(SEQS,"<$tmp_seqs") || die "could not open $tmp_seqs";
679 :     print $cgi->header;
680 :     print "<pre>\n";
681 :     while (defined($_ = <SEQS>))
682 :     {
683 :     print $_;
684 :     }
685 :     print "</pre>\n";
686 :     close(SEQS);
687 :     unlink($tmp_seqs,$tmp_labels);
688 :     }
689 :     }
690 :     sub show_options_for_otu {
691 :     my($cgi,$g) = @_;
692 :    
693 :     print $cgi->header;
694 :     print "<h3>$g</h3>\n";
695 :     print "<ol>\n";
696 : overbeek 1.2 print "<li> <a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_otu_tree&dataD=$g>Gains/losses of Gene Families</a>\n";
697 :     print "<li><a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_kovbassa&dataD=$g>Genes that Distinguishing Families (Genes as Signatures)</a>\n";
698 :     print "<li><a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_adjacency&dataD=$g>Changes in Adjacency</a>\n";
699 : overbeek 1.1 print "</ol>\n";
700 :    
701 :     push(@$html,"<br>",
702 : overbeek 1.2 $cgi->start_form(-action => "./families_on_tree.cgi"),
703 : overbeek 1.1 $cgi->textfield(-name => 'keywords',-size => 100),
704 :     $cgi->hidden(-override => 1,-name => 'request',-value => 'show_indexed_funcs'),
705 :     $cgi->hidden(-override => 1,-name => 'dataD',-value => $dataD),
706 :     $cgi->submit('show functions matching keywords'),$cgi->end_form(),
707 :     "<hr>",
708 :     &virulence_functions_link($cgi,$dataDF));
709 :     foreach $_ (@$html) { print $_ }
710 :    
711 :     my $txt = <<END;
712 :     <pre><br><hr>
713 :     We think of changes as belonging to three broad classes:
714 :    
715 :     1. Genes that are gained or lost,
716 :    
717 :     2. Rearrangements, which we often think of as "changes in adjacency", and
718 :    
719 :     3. SNPs
720 :    
721 :     Users concerned about issues relating to virulence will probably focus on the
722 :     first two classes, unless that have high-quality sequences for a number of
723 :     very close genomes (and they can separate these into virulent and nonvirulent).
724 :    
725 :     On the other hand users interested in what makes a production strain display
726 :     desirable properties may well have a number of very, very close strains, and for them
727 :     SNPs might well be the most interesting form of change.
728 :    
729 :     Lurking behind efforts to accurately display these three types of changes all
730 :     depend heavily on our ability to match up genese between genomes. That is,
731 :     there are (almost inevitably) a set of protein families that are the
732 :     basis of comparison, and these are often of dubious quality. Be aware that
733 :     we have attempted to build fairly accurate families of isofunctional homologs
734 :     (i.e., genes that have the same function and a common ancestor), but they are
735 :     far from perfect.
736 :    
737 :     END
738 :     print $txt;
739 :     }
740 :    
741 :     sub show_otus {
742 :     my($cgi,$datadir) = @_;
743 :    
744 :     print $cgi->header;
745 :     if (opendir(GENERA,$dataDir))
746 :     {
747 :     my @genera = grep { $_ !~ /^\./ } readdir(GENERA);
748 :     closedir(GENERA);
749 :     print "<h1>What Changed?</h1>\n";
750 :     print "<h2><a target=_blank href=\"http://bioseed.mcs.anl.gov/~overbeek/what_changed.html\">Getting Started: a short Tutorial</a></h2>\n";
751 :     print "<h2>Genera Available</h2>\n";
752 :     foreach my $g (sort @genera)
753 :     {
754 : overbeek 1.2 print "<h3><a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_options_for_otu&dataD=$g>$g</a>\n";
755 : overbeek 1.1 }
756 :     }
757 :     else
758 :     {
759 :     print "<h1>The dataD parameter is invalid\n";
760 :     }
761 :     }
762 :    
763 :     sub peg_link {
764 :     my($peg) = @_;
765 :    
766 :     return "<a target=_blank href=http://pubseed.theseed.org/seedviewer.cgi?page=Annotation&feature=$peg>$peg</a>";
767 :     }
768 :    
769 :     sub compare_link {
770 :     my($pegs) = @_;
771 :     my @genomes = map { ($_ =~ /^fig\|(\d+\.\d+)/) ? $1 : () } @$pegs;
772 :     my $args = join("&",map { "show_genome=$_" } @genomes);
773 :     return "<a target=_blank href=http://pubseed.theseed.org/seedviewer.cgi?page=Annotation&feature=" .
774 :     $pegs->[0] . "&$args>Compare Regions</a>";
775 :     }
776 :    
777 :     sub show_func_link {
778 :     my($dataD,$function) = @_;
779 :    
780 :     if (! $function) { return '' }
781 :     my $functionQ = uri_escape($function);
782 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_func&function=$functionQ&dataD=$dataD>$function</a>"
783 : overbeek 1.1 }
784 :    
785 :     sub virulence_functions_link {
786 :     my($cgi,$dataDF) = @_;
787 :    
788 :     if ((-s "$dataDF/virulence.functions") && ($dataDF =~ /([^\/]+)$/))
789 :     {
790 :     my $dataDQ = uri_escape($1);
791 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_virulence_functions&dataD=$dataDQ>Some Posssible Functions Associated with Virulence</a>";
792 : overbeek 1.1 }
793 :     return '';
794 :     }
795 :    
796 :     sub show_fam_on_tree_link {
797 :     my($dataDF,$fam,$union) = @_;
798 :    
799 :     $union = defined($union) ? $union : 0;
800 :     if ($dataDF =~ /([^\/]+)$/)
801 :     {
802 :     my $dataDQ = uri_escape($1);
803 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?family=$fam&union=$union&dataD=$dataDQ&request=show_family_tree>Family on Tree</a>";
804 : overbeek 1.1 }
805 :     return '';
806 :     }
807 :    
808 :     sub show_node_link {
809 :     my($dataD,$node,$type) = @_;
810 :     if ($type eq "families")
811 :     {
812 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_node&dataD=$dataD&node=$node&type=families>$node</a>";
813 : overbeek 1.1 }
814 :     else
815 :     {
816 : overbeek 1.2 return "<a target=_blank href=http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_node&dataD=$dataD&node=$node&type=adjacency>$node</a>";
817 : overbeek 1.1 }
818 :     }
819 :    
820 :     sub show_clusters_link {
821 :     my($dataD,$family,$coupled) = @_;
822 :     my %tmp = map { $_ => 1 } ($family,@$coupled);
823 :     my $families = join(",",sort { $a <=> $b} keys(%tmp));
824 : overbeek 1.2 return "<a target=_blank href=\"http://bioseed.mcs.anl.gov/ross/FIG/families_on_tree.cgi?request=show_clusters&dataD=$dataD&families=$families\">$family</a>";
825 : overbeek 1.1 }
826 :    
827 :     sub show_clusters {
828 :     my($cgi,$dataDF,$html) = @_;
829 :    
830 :     my $families = $cgi->param('families');
831 :     my @families = split(/,/,$families);
832 :     my %families = map { $_ => 1 } @families;
833 :     my %genome_names = map { ($_ =~ /^(\S+)\t(\S.*\S)/) ? ($1 => $2) : () } `cat $dataDF/genome.names`;
834 :     my @genome_pegN_fam_func = sort { ($a->[0] <=> $b->[0]) or ($a->[1] <=> $b->[1]) }
835 :     map { (($_ =~ /^(\S+)\t([^\t]*)\t[^\t]*\tfig\|(\d+\.\d+)\.peg\.(\d+)/) && $families{$1}) ?
836 :     [$3,$4,$1,$2] : ()
837 :     } `cat $dataDF/families.all`;
838 :     push(@$html,$cgi->h1('Relevant Clusters'));
839 :     my $col_hdrs = ['Family','Function','PEG'];
840 :     my $last = shift @genome_pegN_fam_func;
841 :     while ($last)
842 :     {
843 :     my $last_g = $last->[0];
844 :     my $last_pegN = $last->[1];
845 :     my @set;
846 :     while ($last && ($last_g == $last->[0]) && &close($last_pegN,$last->[1]))
847 :     {
848 :     $last_pegN = $last->[1];
849 :     push(@set,[$last->[2],$last->[3],&peg_link("fig|" . $last_g . ".peg." . $last_pegN)]);
850 :     $last = shift @genome_pegN_fam_func;
851 :     }
852 :     if (@set > 1)
853 :     {
854 :     push(@$html,&HTML::make_table($col_hdrs,\@set,"Cluster for $last_g: $genome_names{$last_g}"));
855 :     push(@$html,"<hr><br><br>\n");
856 :     }
857 :     }
858 :     }
859 :    
860 :     sub close {
861 :     my($pegN1,$pegN2) = @_;
862 :    
863 :     return abs($pegN2 - $pegN1) <= 7;
864 :     }
865 :    
866 :     sub show_kovbassa_signatures {
867 :     my($cgi,$dataDF,$html) = @_;
868 :    
869 :     $dataDF =~ /([^\/]+)$/;
870 :     my $dataD = $1;
871 :     my @tree;
872 :     my @tmp = `cat $dataDF/readable.tree`;
873 :     foreach $_ (@tmp)
874 :     {
875 :     if ($_ =~ /^(.*\+ )(n\d+)$/)
876 :     {
877 :     my $start = $1;
878 :     my $node = $2;
879 :     my $link = $cgi->radio_group(-nolabels => 1,
880 :     -name => "radio_$node",
881 :     -override => 1,
882 :     -values => [1,0,2],
883 :     -default => 0);
884 :     push(@tree,"$start$node$link\n");
885 :     }
886 :     elsif ($_ =~ /^(.*\- )(\d+\.\d+)(:.*)$/)
887 :     {
888 :     my $start = $1;
889 :     my $node = $2;
890 :     my $end = $3;
891 :     my $link = $cgi->radio_group(-nolabels => 1,
892 :     -name => "radio_$node",
893 :     -override => 1,
894 :     -values => [1,0,2],
895 :     -default => 0);
896 :     push(@tree,"$start",$link,"$end\n");
897 :     }
898 :     else
899 :     {
900 :     push(@tree,$_);
901 :     }
902 :     }
903 :     my $tree = join("",@tree);
904 :     push(@$html,"<br>",
905 : overbeek 1.2 $cgi->start_form(-action => "./families_on_tree.cgi"),
906 : overbeek 1.1 $cgi->hidden(-override => 1,-name => 'request',-value => 'compute_kovbassa_sigs'),
907 :     $cgi->hidden(-override => 1,-name => 'dataD',-value => $dataD),
908 :     "<pre>\n$tree\n</pre>\n",
909 :     $cgi->submit('Compute Gene Signatures'));
910 :     }
911 :    
912 :     sub compute_kovbassa_signatures {
913 :     my($cgi,$dataDF,$html) = @_;
914 :    
915 :     $dataDF =~ /([^\/]+)$/;
916 :     my $dataD = $1;
917 :     my @out;
918 :     my @in;
919 :     my @params = grep { $_ =~ /^radio/ } $cgi->param();
920 :     foreach my $param (@params)
921 :     {
922 :     if ($param =~ /^radio_(\S+)/)
923 :     {
924 :     my $node = $1;
925 :     if ($cgi->param($param) == 1)
926 :     {
927 :     push(@out,$node);
928 :     }
929 :     elsif ($cgi->param($param) == 2)
930 :     {
931 :     push(@in,$node);
932 :     }
933 :     }
934 :     }
935 :     my $families = "$dataDF/families.all";
936 :     my $tree = &parse_newick_tree(join("",`cat $dataDF/labeled.tree`));
937 :     my $indexP = &tree_utilities::tree_index_tables($tree);
938 :     my $genomes = {};
939 :     &in_set(\@out,$tree,$indexP,1,$genomes);
940 :     &in_set(\@in,$tree,$indexP,2,$genomes);
941 :     my $s1N = grep { $genomes->{$_} == 1 } keys(%$genomes);
942 :     my $s2N = grep { $genomes->{$_} == 2 } keys(%$genomes);
943 :    
944 :     my %by_family;
945 :     my %fam2func;
946 :     foreach $_ (map { (($_ =~ /^(\S+)\t([^\t]*)\tfig\|(\d+\.\d+)/) && $genomes->{$3}) ? [$1,$2,$3] : () }
947 :     `cut -f1,2,4 $families`)
948 :     {
949 :     my($fam,$func,$g) = @$_;
950 :     $by_family{$fam}->{$g} = 1;
951 :     $fam2func{$fam} = $func;
952 :     }
953 :    
954 :     my @scored_families;
955 :     my %famH;
956 :     foreach my $f (keys(%by_family))
957 :     {
958 :     my $fH = $by_family{$f};
959 :     my @hits = keys(%$fH);
960 :     my $s1Hits = grep { $genomes->{$_} == 1 } @hits;
961 :     my $s2Hits = grep { $genomes->{$_} == 2 } @hits;
962 :     my $counts1 = [$s1Hits,$s1N-$s1Hits];
963 :     my $counts2 = [$s2Hits,$s2N-$s2Hits];
964 :     my($din,$dout) = &discriminates($counts2,$counts1);
965 :     my $disc = $din+$dout;
966 :     if ($disc > 1.5)
967 :     {
968 :     $famH{$f} = 1;
969 :     push(@scored_families,[sprintf("%0.3f",$disc),$f]);
970 :     }
971 :     }
972 :     my $coupledH = &coupling_data($dataDF,\%famH);
973 :     my $col_hdrs = ['Score','Tree','Family','Mean Length','Function','Clusters','Coupling'];
974 :     my %mean_len = map { ($_ =~ /^(\S+)\t(\S+)/) ? ($1 => $2) : () } `cut -f1,6 $dataDF/families.all`;
975 :     my $rows = [];
976 :     foreach my $tuple (sort { $b->[0] <=> $a->[0] } @scored_families)
977 :     {
978 :     my($sc,$fam) = @$tuple;
979 :     my $func = $fam2func{$fam};
980 :     my $func_link = &show_func_link($dataD,$func);
981 :     my($cluster_link,$coupled_html) = &cluster_link_and_cluster_html($fam,$coupledH,\%fam2func,$dataD);
982 :     push(@$rows,[$sc,
983 :     &show_fam_links($dataDF,$fam),
984 :     $mean_len{$fam},
985 :     $func_link,
986 :     $cluster_link,
987 :     $coupled_html
988 :     ]);
989 :     }
990 :     push(@$html,&HTML::make_table($col_hdrs,$rows,"Families that Distinguish"));
991 :     }
992 :    
993 :     sub in_set {
994 :     my($nodes,$tree,$indexP,$which,$genomes) = @_;
995 :    
996 :     foreach my $x (@$nodes)
997 :     {
998 :     if ($x =~ /^(\d+\.\d+)/)
999 :     {
1000 :     $genomes->{$1} = $which;
1001 :     }
1002 :     elsif (($x =~ /^(n\d+)/) && (my $node = &label_to_node($indexP,$1)))
1003 :     {
1004 :     my $tips = &tips_of_tree($node);
1005 :     foreach my $tip (@$tips)
1006 :     {
1007 :     $genomes->{$tip} = $which;
1008 :     }
1009 :     }
1010 :     }
1011 :     }
1012 :    
1013 :     sub discriminates {
1014 :     my($xin,$xout) = @_;
1015 :    
1016 :     my $sx = &vector_sum($xin);
1017 :     my $sy = &vector_sum($xout);
1018 :     my $xy = &scalar_product($xin,$xout);
1019 :     my $xx = &scalar_product($xin,$xin);
1020 :     my $yy = &scalar_product($xout,$xout);
1021 :    
1022 :     my $din = (($sx != 0) && ($yy != 0)) ? (1 - (($sy * $xy) / ($sx * $yy))) : 0;
1023 :     my $dout = (($sy != 0) && ($xx != 0)) ? (1 - (($sx * $xy) / ($sy * $xx))) : 0;
1024 :     return ($din,$dout);
1025 :     }
1026 :    
1027 :     sub vector_sum {
1028 :     my($v) = @_;
1029 :    
1030 :     my $sum = 0;
1031 :     foreach $_ (@$v)
1032 :     {
1033 :     $sum += $_;
1034 :     }
1035 :     return $sum;
1036 :     }
1037 :    
1038 :     sub scalar_product {
1039 :     my($x,$y) = @_;
1040 :    
1041 :     if (@$x != @$y)
1042 :     {
1043 :     print STDERR &Dumper($x,$y);
1044 :     die "incompatable";
1045 :     }
1046 :     my $i;
1047 :     my $sum = 0;
1048 :     for ($i=0; ($i < @$x); $i++)
1049 :     {
1050 :     $sum += $x->[$i] * $y->[$i];
1051 :     }
1052 :     return $sum;
1053 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3