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

Annotation of /FigKernelScripts/compute_atomic_regulons_for_dir.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3