[Bio] / FigKernelScripts / compute_atomic_regulons.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/compute_atomic_regulons.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #########################################################################
2 :     # -*- perl -*-
3 :     #
4 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
5 :     # for Interpretations of Genomes. All Rights Reserved.
6 :     #
7 :     # This file is part of the SEED Toolkit.
8 :     #
9 :     # The SEED Toolkit is free software. You can redistribute
10 :     # it and/or modify it under the terms of the SEED Toolkit
11 :     # Public License.
12 :     #
13 :     # You should have received a copy of the SEED Toolkit Public License
14 :     # along with this program; if not write to the University of Chicago
15 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
16 :     # Genomes at veronika@thefig.info or download a copy from
17 :     # http://www.theseed.org/LICENSE.TXT.
18 :     #
19 :    
20 :     use DB_File;
21 :     use FileHandle;
22 :     use File::stat;
23 :     use File::Basename;
24 :    
25 :     use SSserver;
26 :     use SeedEnv;
27 :     use Data::Dumper;
28 :     use Statistics::Descriptive;
29 :     use SAPserver;
30 :     my $sapO = SAPserver->new;
31 :     use strict;
32 :    
33 :     opendir(ORGS,"$FIG_Config::organisms") || die "could not open $FIG_Config::organisms";
34 :     my @genomes = grep { ($_ =~ /^\d+\.\d+$/) &&
35 : overbeek 1.2 (-s "$FIG_Config::organisms/$_/UserSpace/ExpressionData/final_on_off_calls.txt" ) &&
36 :     (! -s "$FIG_Config::organisms/$_/UserSpace/ExpressionData/mix.values" )
37 :     }
38 : overbeek 1.1 readdir(ORGS);
39 :     closedir(ORGS);
40 :    
41 :     my $mix = {};
42 :     foreach my $genome (@genomes)
43 :     {
44 :     &process_genome($sapO,$genome,$mix);
45 :     }
46 :     &build_scored_corr($sapO,$mix);
47 :    
48 :     sub process_genome {
49 :     my($sapO,$genome,$mix) = @_;
50 :     &build_atomic_regulons($sapO,$genome);
51 :     &build_experiment_vectors($sapO,$genome);
52 :     &compute_exp_distances($sapO,$genome);
53 :     &build_ar_vectors($sapO,$genome);
54 :     &compute_ar_distances($sapO,$genome);
55 :     &build_ar_trees($sapO,$genome);
56 :     &build_exp_trees($sapO,$genome);
57 :     &build_mix_values($sapO,$genome,$mix);
58 :     }
59 :    
60 :     sub build_scored_corr {
61 :     my($sapO,$mixH) = @_;
62 :    
63 :     my $dir = "/home/overbeek/Proj/FIGdisk/FIG/Data/ExpressionData";
64 :     my %pearson_hash;
65 : overbeek 1.2 my $pearson_hash_tie = tie %pearson_hash,'DB_File',"$dir/pearson.DB",O_RDONLY,0666,$DB_BTREE;
66 : overbeek 1.1 my @sets = map { chomp; [split(/,/,$_)] } `cat $dir/gene.sets`;
67 :     my $sets = \@sets;
68 :     my $setI;
69 :     my %peg_to_set;
70 :     for ($setI=0; ($setI < @$sets); $setI++)
71 :     {
72 :     my $set = $sets->[$setI];
73 :     foreach my $peg (@$set)
74 :     {
75 :     $peg_to_set{$peg} = $setI;
76 :     }
77 :     }
78 :     my $corrH;
79 :     foreach my $peg (map { my($genome,$pegN) = split(/,/,$_); "fig|$genome.$pegN" } keys(%pearson_hash))
80 :     {
81 :     my($set1,$set2);
82 :     if (($mix->{$peg}) && defined($set1 = $peg_to_set{$peg}))
83 :     {
84 :     my $connH = &conn_by_pc($peg,\%pearson_hash);
85 :     foreach my $peg2 (keys(%$connH))
86 :     {
87 :     if ((my $x = $connH->{$peg2}) && ($set2 = $peg_to_set{$peg2}))
88 :     {
89 :     $x = $x * $mix->{$peg};
90 :     if (! $corrH->{$set1}->{$set2})
91 :     {
92 :     $corrH->{$set1}->{$set2} = $x;
93 :     }
94 :     else
95 :     {
96 :     $corrH->{$set1}->{$set2} += $x;
97 :     }
98 :     }
99 :     }
100 :     }
101 :     }
102 :    
103 :     open(CORR,">$dir/sets.corr") || die "could not open $dir/sets.corr";
104 :     foreach my $set1 (sort {$a <=> $b } keys(%$corrH))
105 :     {
106 :     my $hash = $corrH->{$set1};
107 :     foreach my $set2 (keys(%$hash))
108 :     {
109 :     my $sc = sprintf("%0.3f",$hash->{$set2});
110 :     if ($sc > 0.001)
111 :     {
112 :     print CORR join("\t",$set1,$set2,$sc),"\n";
113 :     }
114 :     }
115 :     }
116 :     close(CORR);
117 :     }
118 :    
119 :     sub build_mix_values {
120 :     my($sapO,$genome,$mixH) = @_;
121 :    
122 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
123 :     open(MIX,">$expD/mix.values") || die "could not open mix.values";
124 :     my($vecs,$sz1) = &read_on_off($genome);
125 :     foreach my $peg (keys(%$vecs))
126 :     {
127 :     my $v = $vecs->{$peg};
128 :     my $on = &v_occurs($v,1);
129 :     my $off = &v_occurs($v,-1);
130 :     my $min = ($on > $off) ? $off : $on;
131 :     my $max = ($on < $off) ? $off : $on;
132 :     my $mix = sprintf("%0.3f", (($min+1)/($max+1)));
133 :     print MIX "$peg\t$mix\n";
134 :     $mixH->{$peg} = $mix;
135 :     }
136 :     }
137 :    
138 :     sub v_occurs {
139 :     my($v,$x) = @_;
140 :    
141 :     my $tot = 0;
142 :     foreach $_ (@$v)
143 :     {
144 :     if ($_ eq $x) { $tot++ }
145 :     }
146 :     return $tot;
147 :     }
148 :    
149 :     use Clustering;
150 :     use tree_utilities;
151 :     sub build_ar_trees {
152 :     my($sapO,$genome) = @_;
153 :    
154 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
155 :     my $linkage = "avg_dist";
156 :    
157 :     my %conn;
158 :     my %reg;
159 :     foreach $_ (`cat $expD/ar.vectors.distances`)
160 :     {
161 :     if ($_ =~ /^(\d+)\t(\d+)\t(\S+)/)
162 :     {
163 :     $conn{$1}->{$2} = $conn{$2}->{$1} = $3;
164 :     }
165 :     }
166 :    
167 :     my $maxD;
168 :     for ($maxD=0.05; ($maxD <= 0.35); $maxD += 0.05)
169 :     {
170 :     my($clusters,$trees) = &Clustering::cluster(\%conn,$maxD,$linkage);
171 :     open(TREES,">$expD/ar.trees-$maxD.nwk") || die "could not open $expD/ar.trees-$maxD.nwk";
172 :     foreach my $tree (@$trees)
173 :     {
174 :     print TREES &tree_utilities::to_newick($tree),"\n//\n";
175 :     }
176 :     close(TREES);
177 :     open(CLUST,">$expD/ar.clusters-$maxD") || die "could not open $expD/ar.clusters-$maxD";
178 :     foreach my $clust (@$clusters)
179 :     {
180 :     print CLUST join(",",@$clust),"\n";
181 :     }
182 :     close(CLUST);
183 :     }
184 :     }
185 :    
186 :     sub build_exp_trees {
187 :     my($sapO,$genome) = @_;
188 :    
189 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
190 :     my $linkage = "avg_dist";
191 :    
192 :     my %conn;
193 :     my %reg;
194 :     foreach $_ (`cat $expD/experiment.vectors.distances`)
195 :     {
196 :     if ($_ =~ /^(\d+)\t(\d+)\t(\S+)/)
197 :     {
198 :     $conn{$1}->{$2} = $conn{$2}->{$1} = $3;
199 :     }
200 :     }
201 :    
202 :     my $maxD;
203 :     for ($maxD=0.05; ($maxD <= 0.35); $maxD += 0.05)
204 :     {
205 :     my($clusters,$trees) = &Clustering::cluster(\%conn,$maxD,$linkage);
206 :     open(TREES,">$expD/experiment.trees-$maxD.nwk") || die "could not open $expD/experiment.trees-$maxD.nwk";
207 :     foreach my $tree (@$trees)
208 :     {
209 :     print TREES &tree_utilities::to_newick($tree),"\n//\n";
210 :     }
211 :     close(TREES);
212 :     open(CLUST,">$expD/experiment.clusters-$maxD") || die "could not open $expD/experiment.clusters-$maxD";
213 :     foreach my $clust (@$clusters)
214 :     {
215 :     print CLUST join(",",@$clust),"\n";
216 :     }
217 :     close(CLUST);
218 :     }
219 :     }
220 :    
221 :     sub read_atomic_regulons {
222 :     my($genome) = @_;
223 :    
224 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
225 :     my $ar = [];
226 :     foreach $_ (`cat $expD/atomic.regulons`)
227 :     {
228 :     if ($_ =~ /^(\d+)\t(\S+)/)
229 :     {
230 :     push(@{$ar->[$1-1]},$2);
231 :     }
232 :     }
233 :     return $ar;
234 :     }
235 :    
236 :     sub compute_exp_distances {
237 :     my($sapO,$genome) = @_;
238 :    
239 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
240 :     open(EXPDIST,">$expD/experiment.vectors.distances") || die "could not open $expD/experiment.vectors.distances";
241 :     my $expV = &read_exp_vecs($genome);
242 :     my @vecIDs = sort { $a <=> $b } keys(%$expV);
243 :     my($i,$j);
244 :     for ($i=0; ($i < $#vecIDs); $i++)
245 :     {
246 :     for ($j=$i+1; ($j < @vecIDs); $j++)
247 :     {
248 :     my $d = &distance($expV,$vecIDs[$i],$vecIDs[$j]);
249 :     print EXPDIST join("\t",($vecIDs[$i],$vecIDs[$j],$d)),"\n";
250 :     print EXPDIST join("\t",($vecIDs[$j],$vecIDs[$i],$d)),"\n";
251 :     }
252 :     }
253 :     close(EXPDIST);
254 :     }
255 :    
256 :     sub compute_ar_distances {
257 :     my($sapO,$genome) = @_;
258 :    
259 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
260 :     open(ARDIST,">$expD/ar.vectors.distances") || die "could not open $expD/ar.vectors.distances";
261 :     my $arV = &read_ar_vecs($genome);
262 :     my @vecIDs = sort { $a <=> $b } keys(%$arV);
263 :     my($i,$j);
264 :     for ($i=0; ($i < $#vecIDs); $i++)
265 :     {
266 :     for ($j=$i+1; ($j < @vecIDs); $j++)
267 :     {
268 :     my $d = &distance($arV,$vecIDs[$i],$vecIDs[$j]);
269 :     print ARDIST join("\t",($vecIDs[$i],$vecIDs[$j],$d)),"\n";
270 :     print ARDIST join("\t",($vecIDs[$j],$vecIDs[$i],$d)),"\n";
271 :     }
272 :     }
273 :     close(EXPDIST);
274 :     }
275 :    
276 :     sub read_ar_vecs {
277 :     my($genome) = @_;
278 :    
279 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
280 :     return {map { $_ =~ /^(\d+)\t(\S+)/; $1 => [split(/,/,$2)] } `cat $expD/ar.vectors`}
281 :     }
282 :    
283 :     sub read_exp_vecs {
284 :     my($genome) = @_;
285 :    
286 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
287 :     return {map { $_ =~ /^(\d+)\t(\S+)/; $1 => [split(/,/,$2)] } `cat $expD/experiment.vectors`}
288 :     }
289 :    
290 :     sub distance {
291 :     my($idH,$v1,$v2) = @_;
292 :    
293 :     # $sim can range from -1 to 1.
294 :     # I am defining distance as 1 - (($sim + 1)/2)
295 :     my $sim = &dot_product(&unit_vector($idH->{$v1}),&unit_vector($idH->{$v2}));
296 :     return sprintf("%0.4f",1 - (($sim + 1)/2));
297 :     }
298 :    
299 :     sub build_ar_vectors {
300 :     my($sapO,$genome) = @_;
301 :    
302 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
303 :     open(ARVEC,">$expD/ar.vectors") || die "could not open $expD/ar.vectors";
304 :    
305 :     my $ar = &read_atomic_regulons($genome);
306 :     my($vecs,$vec_sz) = &read_vecs($genome);
307 :     my $on_off_vecs = $vecs->[0];
308 :     my $i;
309 :     my @vals;
310 :     for ($i=0; ($i < @$ar); $i++)
311 :     {
312 :     my $j;
313 :     for ($j=0; ($j < $vec_sz); $j++)
314 :     {
315 :     $vals[$j] = &ar_status($ar->[$i],$j,$on_off_vecs);
316 :     }
317 :     print ARVEC $i+1,"\t",join(",",@vals),"\n";
318 :     }
319 :     close(ARVEC);
320 :     }
321 :    
322 :    
323 :     sub build_experiment_vectors {
324 :     my($sapO,$genome) = @_;
325 :    
326 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
327 :     open(EXPVEC,">$expD/experiment.vectors") || die "could not open $expD/experiment.vectors";
328 :    
329 :     my $ar = &read_atomic_regulons($genome);
330 :     my($vecs,$vec_sz) = &read_vecs($genome);
331 :     my $on_off_vecs = $vecs->[0];
332 :     my $i;
333 :     my @vals;
334 :     for ($i=0; ($i < $vec_sz); $i++)
335 :     {
336 :     my $j;
337 :     for ($j=0; ($j < @$ar); $j++)
338 :     {
339 :     $vals[$j] = &ar_status($ar->[$j],$i,$on_off_vecs);
340 :     }
341 :     print EXPVEC $i+1,"\t",join(",",@vals),"\n";
342 :     }
343 :     close(EXPVEC);
344 :     }
345 :    
346 :     sub ar_status {
347 :     my($ar1,$expI,$on_off_vecs) = @_;
348 :    
349 :     my @pegs = @$ar1;
350 :     my $on = 0;
351 :     my $off = 0;
352 :     foreach my $peg (@pegs)
353 :     {
354 :     my $v = $on_off_vecs->{$peg}->[$expI];
355 :     if ($v == 1) { $on++ }
356 :     elsif ($v == -1) { $off++ }
357 :     }
358 :    
359 :     if ($on > $off) { return 1 }
360 :     elsif ($on < $off) { return -1 }
361 :     else { return 0 }
362 :     }
363 :    
364 :     sub dot_product {
365 :     my($v1,$v2) = @_;
366 :    
367 :     $v1 = &unit_vector($v1);
368 :     $v2 = &unit_vector($v2);
369 :    
370 :     my $tot = 0;
371 :     my $i;
372 :     for ($i=0; ($i < @$v1); $i++)
373 :     {
374 :     if ($v1->[$i] && $v2->[$i])
375 :     {
376 :     $tot += $v1->[$i] * $v2->[$i];
377 :     }
378 :     }
379 :     return $tot;
380 :     }
381 :    
382 :     sub unit_vector {
383 :     my($v) = @_;
384 :    
385 :     my $tot = 0;
386 :     my $uv = [];
387 :     my $i;
388 :     for ($i = 0; ($i < @$v); $i++)
389 :     {
390 :     my $x = $v->[$i];
391 :     if (defined($x))
392 :     {
393 :     $tot += $x * $x;
394 :     }
395 :     }
396 :    
397 :     my $nf = sqrt($tot);
398 :     for ($i = 0; ($i < @$v); $i++)
399 :     {
400 :     my $x = $v->[$i];
401 :     $x = $x ? $x : 0;
402 :     push(@$uv,$nf ? sprintf("%0.4f",($x / $nf)) : 0);
403 :     }
404 :     return $uv;
405 :     }
406 :    
407 :     sub build_atomic_regulons {
408 :     my($sapO,$genome) = @_;
409 :    
410 :     my $genH = $sapO->all_features( -ids => [$genome], -type => ['peg'] );
411 :     my $all_pegs = $genH->{$genome};
412 :     my $pegH = $sapO->ids_to_functions( -ids => $all_pegs );
413 :    
414 :     my $expD = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
415 :     open(AR,">$expD/atomic.regulons") || die "could not open $expD/atomic.regulons";
416 :    
417 :     my($vecs,$vec_sz) = &read_vecs($genome);
418 :     my $on_off_vecs = $vecs->[0];
419 :     my $clusters = &atomic_regulons($genome);
420 :     ##
421 :     ## Here we keep only conjectures that have few contradictions. I set the value to 1 % of the arrays
422 :     ##
423 :     my @perfect = sort { @{$b->[1]} <=> @{$a->[1]} }
424 :     map { (($_->[0] / $vec_sz) <= 0.01) ? $_->[1] : () } &scored_clusters($clusters,$vecs);
425 :    
426 :     my @ar = &condense_sets(\@perfect);
427 :    
428 :     my $nxt = 1; ## I flipped and decided to renumber the sets ##
429 :     @ar = sort { @{$b->[1]} <=> @{$a->[1]} } @ar;
430 :     @ar = map { $_->[1] } @ar;
431 :    
432 :     #
433 :     # Now I am going to put all PEGs with identical on/off values into
434 :     # the same set by 1) merging sets containing them and 2) adding pegs
435 :     # that were not put into atomic regulons
436 :     #
437 :     @ar = &extend_ar($genome,$on_off_vecs,\@ar);
438 :     foreach my $pegs (sort { @$b <=> @$a } @ar)
439 :     {
440 :     foreach my $peg (@$pegs)
441 :     {
442 :     print AR join("\t",($nxt,$peg,&func_of($pegH,$peg))),"\n";
443 :     }
444 :     $nxt++;
445 :     }
446 :     close(AR);
447 :     }
448 :    
449 :     # $sets points to a list of 2-tuples [$set,$pegs]
450 :     sub condense_sets {
451 :     my($sets) = @_;
452 :    
453 :     my $set_to_pegs = {};
454 :     my $peg_to_sets = {};
455 :     my %sz;
456 :    
457 :     foreach my $tuple (@$sets)
458 :     {
459 :     my($set,$pegs) = @$tuple;
460 :     $sz{$set} = @$pegs;
461 :     foreach my $peg (@$pegs)
462 :     {
463 :     $set_to_pegs->{$set}->{$peg} = 1;
464 :     $peg_to_sets->{$peg}->{$set} = 1;
465 :     }
466 :     }
467 :    
468 :     foreach my $peg (keys(%$peg_to_sets))
469 :     {
470 :     my @in = sort { $sz{$b} <=> $sz{$a} } keys(%{$peg_to_sets->{$peg}});
471 :     if (@in > 1)
472 :     {
473 :     my $i;
474 :     for ($i=1; ($i < @in); $i++)
475 :     {
476 :     my @pegs1 = keys(%{$set_to_pegs->{$in[0]}});
477 :     my @pegs2 = keys(%{$set_to_pegs->{$in[$i]}});
478 :     my @pegs = (@pegs1,@pegs2);
479 :    
480 :     if (1) # (&compatible(\@pegs))
481 :     {
482 :     foreach my $peg1 (@pegs2)
483 :     {
484 :     $set_to_pegs->{$in[0]}->{$peg1} = 1;
485 :     delete $set_to_pegs->{$in[$i]}->{$peg1};
486 :     $peg_to_sets->{$peg1}->{$in[0]} = 1;
487 :     delete $peg_to_sets->{$peg1}->{$in[$i]};
488 :     }
489 :     }
490 :     }
491 :     }
492 :     }
493 :    
494 :     my @ar;
495 :     foreach my $set (keys(%$set_to_pegs))
496 :     {
497 :     my @pegs = sort { &SeedUtils::by_fig_id($a,$b) } keys(%{$set_to_pegs->{$set}});
498 :     if (@pegs > 1)
499 :     {
500 :     push(@ar,[$set,\@pegs]);
501 :     }
502 :     }
503 :     return @ar;
504 :     }
505 :    
506 :     sub compatible {
507 :     my($pegs) = @_;
508 :    
509 :     my $genome = &SeedUtils::genome_of($pegs->[0]);
510 :    
511 :     my $dir = "$FIG_Config::data/ExpressionData";
512 :     my %pearson_hash;
513 :     my $pearson_hash_tie = tie %pearson_hash,'DB_File',"$dir/pearson.DB",O_RDONLY,0666,$DB_BTREE;
514 :    
515 :     my ($i,$j);
516 :     my $ok = 1;
517 :     for ($i=0; $ok && ($i < @$pegs); $i++)
518 :     {
519 :     my $hash = &conn_by_pc($pegs->[$i],\%pearson_hash);
520 :     for ($j=0; $ok && ($j < @$pegs); $j++)
521 :     {
522 :     if ($pegs->[$i] ne $pegs->[$j])
523 :     {
524 :     my $pc = $hash->{$pegs->[$j]};
525 :     if (! (defined($pc = $hash->{$pegs->[$j]}) && ($pc >= 0.5)))
526 :     {
527 :     $ok = 0;
528 :     }
529 :     }
530 :     }
531 :     }
532 :     undef $pearson_hash_tie;
533 :     untie %pearson_hash;
534 :     return $ok;
535 :     }
536 :    
537 :     sub extend_ar {
538 :     my($genome,$on_off_vecs,$ars) = @_;
539 :    
540 :     my %by_on_off;
541 :     while (my($peg,$vals) = each(%$on_off_vecs))
542 :     {
543 :     push(@{$by_on_off{join("",@$vals)}},$peg);
544 :     }
545 :     my @new_sets = map { my $x = $by_on_off{$_}; (@$x > 1) ? $x : () } keys(%by_on_off);
546 :     # print STDERR &Dumper([map { &member('fig|300852.3.peg.1',$_) ? $_ : () } @new_sets]); die "aborted";
547 :     my $n = 1;
548 :     my @all = map { [$n++,$_] } (@$ars,@new_sets);
549 :     my @condensed = map { $_->[1] } &condense_sets(\@all);
550 :     return @condensed;
551 :     }
552 :    
553 :     sub func_of {
554 :     my($pegH,$peg) = @_;
555 :    
556 :     my $func = $pegH->{$peg};
557 :     $func = $func ? $func : 'hypothetical protein';
558 :     return $func;
559 :     }
560 :    
561 :     sub scored_clusters {
562 :     my($clusters,$vecs) = @_;
563 :     my @scored = map { [&scored_cluster($vecs,$_->[1]),$_] } @$clusters;
564 :     return @scored;
565 :     }
566 :    
567 :     sub scored_cluster {
568 :     my($vecs,$pegs) = @_;
569 :     my($i,$j);
570 :     my $high_sc = 0;
571 :    
572 :     for ($i=0; ($i < (@$pegs-1)); $i++)
573 :     {
574 :     for ($j=$i+1; ($j < @$pegs); $j++)
575 :     {
576 :     my($v1,$v2);
577 :     if (($v1 = $vecs->[0]->{$pegs->[$i]}) &&
578 :     ($v2 = $vecs->[0]->{$pegs->[$j]}))
579 :     {
580 :     my $sc = &scored_pair($v1,$v2);
581 :    
582 :     if ($sc > $high_sc)
583 :     {
584 :     $high_sc = $sc;
585 :     }
586 :     }
587 :     }
588 :     }
589 :     return $high_sc;
590 :     }
591 :    
592 :     sub scored_pair {
593 :     my($v1,$v2) = @_;
594 :    
595 :     if ((! $v1) || (! $v2)) { return 0 }
596 :    
597 :     my $i;
598 :     my $match = 0;
599 :     my $mismatch = 0;
600 :     my $sc = 0;
601 :     for ($i=0; ($i < @$v1); $i++)
602 :     {
603 :     if (($v1->[$i] != 0) && ($v2->[$i] != 0))
604 :     {
605 :     if ($v1->[$i] != $v2->[$i])
606 :     {
607 :     $mismatch++;
608 :     }
609 :     else
610 :     {
611 :     $match++;
612 :     }
613 :     }
614 :     }
615 :     return ($match+$mismatch) ? int(($mismatch * 100) / ($match+$mismatch)) : 0;
616 :     }
617 :    
618 :    
619 :     sub member {
620 :     my($x,$xL) = @_;
621 :    
622 :     my $i;
623 :     for ($i=0; ($i < @$xL) && ($x ne $xL->[$i]); $i++) {}
624 :     return ($i < @$xL);
625 :     }
626 :    
627 :    
628 :     sub read_vecs {
629 :     my($genome) = @_;
630 :    
631 :     my($vecs1,$sz1) = &read_on_off($genome);
632 :     my($vecs2,$sz2) = &read_normalized_values($genome);
633 :     if ($sz1 != $sz2) { die "sz1=$sz1 sz2=$sz2" }
634 :     my $cutoffs = &cutoffs($genome);
635 :     return ([$vecs1,$vecs2,$cutoffs],$sz1);
636 :     }
637 :    
638 :     ###################
639 :    
640 :     sub atomic_regulons {
641 :     my($genome) = @_;
642 :     my $dir = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
643 :     my $n = 1;
644 :     my @clusters = ();
645 :     foreach $_ (`cat $dir/coregulated.*`)
646 :     {
647 :     if ($_ =~ /^(fig\S+)(\t(\S.*\S))?/)
648 :     {
649 :     my $desc = $3 ? $3 : '';
650 :     my $pegs = [split(/,/,$1)];
651 :     push(@clusters,[$n++,$pegs,$desc]);
652 :     }
653 :     }
654 :     return \@clusters;
655 :     }
656 :    
657 :     sub get_pc_hash {
658 :     my($pegs) = @_;
659 :    
660 :     my $corrH = &get_corr(&SeedUtils::genome_of($pegs->[0]));
661 :     my $hash = &compute_pc($pegs,$corrH);
662 :     return $hash;
663 :     }
664 :    
665 :     sub read_on_off {
666 :     my($genome) = @_;
667 :     return &read_vecs1($genome,"final_on_off_calls.txt");
668 :     }
669 :    
670 :     sub read_normalized_values {
671 :     my($genome) = @_;
672 :    
673 :     &read_vecs1($genome,"raw_data.tab");
674 :     }
675 :    
676 :     sub cutoffs {
677 :     my($genome) = @_;
678 :    
679 :     my $dir = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
680 :     return [map { $_ =~ /^(\S+)\s+(\S+)/; [sprintf("%0.3f",$1),sprintf("%0.3f",$2)] } `cat $dir/cutoffs.txt`];
681 :     }
682 :    
683 :     ###################
684 :     sub read_vecs1 {
685 :     my($genome,$file) = @_;
686 :    
687 :     my $dir = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
688 :     my %seen;
689 :     my $vecs = {};
690 :     my $sz;
691 :     foreach $_ (`cat $dir/$file`)
692 :     {
693 :     if ($_ =~ /^(fig\|\d+\.\d+\.(peg|rna)\.\d+)\t(\S.*\S)/)
694 :     {
695 :     my $peg = $1;
696 :     my $vec = [split(/\t/,$3)];
697 :     if (! $sz)
698 :     {
699 :     $sz = @$vec;
700 :     }
701 :     if (! $seen{$peg})
702 :     {
703 :     $seen{$peg} = 1;
704 :     if (@$vec != $sz) { die 'malformed vectors' }
705 :     $vecs->{$peg} = $vec;
706 :     }
707 :     }
708 :     }
709 :     return ($vecs,$sz);
710 :     }
711 :    
712 :     sub conn_by_pc {
713 :     my($peg,$pearson_hash) = @_;
714 :    
715 :     my($genome,$pegN) = ($peg =~ /fig\|(\d+\.\d+)\.((peg|rna)\.\d+)$/);
716 :     my $to = {};
717 :     if ($_ = $pearson_hash->{"$genome,$pegN"})
718 :     {
719 :     $to = {map { $_ =~ /(\S+):(\S+)$/; ("fig\|$genome\.$2" => $1) } split(/,/,$_) };
720 :     }
721 :     return $to;
722 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3