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

Annotation of /FigKernelPackages/FC.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1
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 :     package FC;
19 :    
20 :     use strict;
21 :     use ERDB;
22 :     use Data::Dumper;
23 : overbeek 1.8 use Carp;
24 : overbeek 1.3 use FIG;
25 : overbeek 1.7 use PartitionSeqs;
26 : overbeek 1.1
27 : overbeek 1.12 sub co_occurring_FIGfams {
28 :     my($db,$ff) = @_;
29 :    
30 :     my $query = $db->Get("Family HasMember Feature IsInPair Pairing Determines PairSet","Family(id) = ?",[$ff]);
31 :    
32 :     my(%hits,%pairsets,$row,$pairset,$sc);
33 :     while ($row = $query->Fetch())
34 :     {
35 :     ($pairset,$sc) = $row->Values(['PairSet(id)','PairSet(score)']);
36 :     $pairsets{$pairset} = $sc;
37 :     }
38 :    
39 :     while (($pairset,$sc) = each(%pairsets))
40 :     {
41 : overbeek 1.14 $query = $db->Get("IsDeterminedBy Pairing IsPairOf Feature IsMemberOf Family","IsDeterminedBy(from-link) = ?",[$pairset]);
42 : overbeek 1.12 while ($row = $query->Fetch())
43 :     {
44 : overbeek 1.14 my($ff1,$fam_func) = $row->Values(['IsMemberOf(to-link)'
45 :     # ,'Family(family-function)'
46 :     ]);
47 : overbeek 1.12 if ($ff1 ne $ff)
48 :     {
49 : overbeek 1.14 if ((! $hits{$ff1}) || ($hits{$ff1}->[0] < $sc))
50 : overbeek 1.12 {
51 : overbeek 1.14 $hits{$ff1} = [$sc,$fam_func];
52 : overbeek 1.12 }
53 :     }
54 :     }
55 :     }
56 :     return map { [$_,$hits{$_}] } sort { $hits{$b} <=> $hits{$a} } keys(%hits);
57 :     }
58 :    
59 : overbeek 1.1 sub co_occurs {
60 :     my($db,$peg) = @_;
61 :    
62 :     my @tuples;
63 :     my $query = $db->Get("IsInPair Pairing Determines PairSet",
64 :     "IsInPair(from-link) = ?", [$peg]);
65 :     while (my $row = $query->Fetch())
66 :     {
67 : overbeek 1.13 my ($pair, $pairset, $score) = $row->Values(['Determines(from-link)','PairSet(id)','PairSet(score)']);
68 : overbeek 1.1 my ($fid) = grep { $_ ne $peg } split(/:/, $pair);
69 : overbeek 1.13 push @tuples, [$score,$fid,$pairset];
70 : overbeek 1.1 }
71 :     return @tuples;
72 :     }
73 :    
74 : overbeek 1.15 sub co_occurs_batch {
75 :     my($db,$pegs) = @_;
76 :    
77 :     my @tuples;
78 :     my $query = $db->Get("IsInPair Pairing Determines PairSet",
79 :     "IsInPair(from-link) IN (" . join(", ", map { "'$_'" } @$pegs) . ")", []);
80 :     while (my $row = $query->Fetch())
81 :     {
82 :     my ($peg1,$pair,$score) = $row->Values(['IsInPair(from-link)','Determines(from-link)','PairSet(score)']);
83 :     my($peg2) = grep { $_ ne $peg1 } split(/:/,$pair);
84 :     push(@tuples, [$peg1,$peg2,$score]);
85 :     }
86 :     return wantarray ? @tuples : \@tuples;
87 :     }
88 :    
89 :     sub coupled_to {
90 :     my($db,$peg) = @_;
91 :    
92 :     return map { [$_->[1],$_->[0]] } &FC::co_occurs($db,$peg);
93 :     }
94 :    
95 :     sub coupled_to_batch {
96 :     my($db,@pegs) = @_;
97 :    
98 :     return &FC::co_occurs_batch($db,\@pegs);
99 :     }
100 :    
101 : overbeek 1.9 sub in_co_occurrence_cluster {
102 :     my($db,$peg) = @_;
103 :    
104 :     my $query = $db->Get("OccursIn Cluster IsOccurrenceOf",
105 :     "OccursIn(from-link) = ?", [$peg]);
106 :     my @pegs = ();
107 :     while (my $row = $query->Fetch())
108 :     {
109 :     my($peg1) = $row->Values(['IsOccurrenceOf(to-link)']);
110 :     push(@pegs,$peg1);
111 :     }
112 :     return (@pegs > 0) ? \@pegs : undef;
113 :     }
114 :    
115 : overbeek 1.15 sub co_occurrence_evidence {
116 :     my($db,$peg1,$peg2) = @_;
117 :    
118 :     my $flipped = 0;
119 :     if ($peg1 gt $peg2)
120 :     {
121 :     $flipped = 1;
122 :     ($peg1,$peg2) = ($peg2,$peg1);
123 :     }
124 :     my $pairing_id = join(":",($peg1,$peg2));
125 :     my $query = $db->Get("Pairing Determines PairSet IsDeterminedBy",
126 :     "Pairing (id) = ?", [$pairing_id]);
127 :     my $tuples;
128 :     while (my $row = $query->Fetch())
129 :     {
130 :     my($pair_id,$inverted) = $row->Values(['IsDeterminedBy(to-link)','IsDeterminedBy(inverted)']);
131 :     my($peg3,$peg4) = split(":",$pair_id);
132 :     if ($flipped != $inverted) { ($peg3,$peg4) = ($peg4,$peg3) }
133 :     if ($peg3 ne $peg1)
134 :     {
135 :     push(@$tuples,[$peg3,$peg4]);
136 :     }
137 :     }
138 :     return $tuples;
139 :     }
140 :    
141 :     sub coupling_evidence {
142 :     my($db,$peg1,$peg2) = @_;
143 :    
144 :     my $tuples = &co_occurrence_evidence($db,$peg1,$peg2);
145 :     if (! $tuples) { return () }
146 :    
147 :     my $i;
148 :     my %seen;
149 :     for ($i=0; ($i < @$tuples); $i++)
150 :     {
151 :     my($peg3,$peg4) = @{$tuples->[$i]};
152 :     $peg3 =~ /^fig\|(\d+\.\d+)/;
153 :     my $genome = $1;
154 :     my $otu = $db->OTU($genome);
155 :     if (! $seen{$otu})
156 :     {
157 :     $seen{$otu} = 1;
158 :     push(@{$tuples->[$i]},1);
159 :     }
160 :     else
161 :     {
162 :     push(@{$tuples->[$i]},0);
163 :     }
164 :     }
165 :     return wantarray ? @{$tuples} : $tuples;
166 :     }
167 :    
168 : overbeek 1.1 sub co_occurrence_and_evidence {
169 :     my($db,$peg) = @_;
170 :    
171 :     my @tuples;
172 :     my $query = $db->Get("IsInPair Pairing Determines PairSet IsDeterminedBy",
173 :     "IsInPair(from-link) = ?", [$peg]);
174 :     my @tuoples = ();
175 :     while (my $row = $query->Fetch())
176 :     {
177 :     push(@tuples,[$row->Values(['Determines(from-link)',
178 :     'PairSet(score)',
179 : overbeek 1.3 'IsDeterminedBy(inverted)',
180 : overbeek 1.1 'IsDeterminedBy(to-link)'])]);
181 :     }
182 : overbeek 1.2 @tuples = sort { $a->[0] <=> $b->[0] } @tuples;
183 :     my @co_occurs = ();
184 :    
185 :     my $x = shift @tuples;
186 :     while ($x)
187 :     {
188 :     my $curr = $x->[0];
189 :     my $sc = $x->[1];
190 :     my @set = ();
191 :     while ($x && ($x->[0] eq $curr))
192 :     {
193 :     my($peg3,$peg4) = split(/:/,$x->[3]);
194 :     push(@set,($x->[2]) ? [$peg3,$peg4] : [$peg4,$peg3]);
195 :     $x = shift @tuples;
196 :     }
197 :     my $i;
198 :     for ($i=0; ($i < @set) && ($peg ne $set[$i]->[0]); $i++) {}
199 :     if ($i == @set)
200 :     {
201 :     @set = map { [$_->[1],$_->[0]] } @set;
202 :     }
203 :     my($peg2) = grep { $_ ne $peg } split(/:/,$curr);
204 : overbeek 1.3 push(@co_occurs,[$sc,$peg2,[grep { $_->[0] ne $peg } @set]]);
205 : overbeek 1.2 }
206 :     return sort { $b->[0] <=> $a->[0] } @co_occurs;
207 : overbeek 1.1 }
208 :    
209 :     sub co_occurrence_set {
210 :     my($db,$set) = @_;
211 :    
212 :     my $pairs;
213 :     my $sc;
214 : overbeek 1.8
215 :     my $query = $db->Get("PairSet",
216 :     "PairSet(id) = ?", [$set]);
217 :     my $row = $query->Fetch();
218 :     if ($row)
219 : overbeek 1.1 {
220 : overbeek 1.8 ($sc) = $row->Values(['PairSet(score)']);
221 :     $pairs = [];
222 :     my $query = $db->Get("PairSet IsDeterminedBy",
223 :     "PairSet(id) = ?", [$set]);
224 :     while (my $row = $query->Fetch())
225 :     {
226 :     my($pairing,$inverted) = $row->Values(['IsDeterminedBy(to-link)','IsDeterminedBy(inverted)']);
227 :     my($peg1,$peg2) = split(/:/,$pairing);
228 :     push(@$pairs,$inverted ? [$peg2,$peg1] : [$peg1,$peg2]);
229 :     }
230 : overbeek 1.1 }
231 :     return ($sc,$pairs);
232 :     }
233 :    
234 : overbeek 1.3 sub all_co_occurrence_pair_sets {
235 :     my($db) = @_;
236 :    
237 :     my $query = $db->Get("PairSet","",[]);
238 :     my @PairSets = ();
239 :     while (my $row = $query->Fetch())
240 :     {
241 :     my($pair_set) = $row->Values(['PairSet(id)']);
242 :     push(@PairSets,$pair_set);
243 :     }
244 :     return @PairSets;
245 :     }
246 :    
247 :     sub all_co_occurrence_clusters {
248 :     my($db) = @_;
249 :    
250 :     my $query = $db->Get("Cluster","",[]);
251 :     my @Clusters = ();
252 :     while (my $row = $query->Fetch())
253 :     {
254 :     my($id) = $row->Values(['Cluster(id)']);
255 :     push(@Clusters,$id);
256 :     }
257 :     return @Clusters;
258 :     }
259 :    
260 : overbeek 1.9 sub largest_co_occurrence_clusters {
261 :     my($db,$peg) = @_;
262 :    
263 :     my $fig = new FIG;
264 : overbeek 1.11 my %pegs = map { $_->id2 => $_->psc } $fig->sims($peg,1000,1.0e-30,'fig');
265 : overbeek 1.9
266 : overbeek 1.11 my @pegs = sort { $pegs{$a} <=> $pegs{$b} } grep { $_ ne $peg } keys(%pegs);
267 :     unshift @pegs,$peg;
268 :     $pegs{$peg} = 0;
269 : overbeek 1.9 my @clusters = ();
270 :     foreach my $peg1 (@pegs)
271 :     {
272 :     my $cluster = &in_co_occurrence_cluster($db,$peg1);
273 :     if (defined($cluster))
274 :     {
275 : overbeek 1.11 my $tuple = [$peg1,$pegs{$peg1}, [grep { $_ ne $peg1 } @$cluster]];
276 : overbeek 1.9 push(@clusters,$tuple);
277 :     }
278 :     }
279 : overbeek 1.11 return sort { @{$b->[2]} <=> @{$a->[2]} } @clusters;
280 : overbeek 1.9 }
281 :    
282 : overbeek 1.1 sub co_occurrence_cluster {
283 :     my($db,$cluster) = @_;
284 :    
285 :     my $query = $db->Get("Cluster IsOccurrenceOf",
286 :     "Cluster(id) = ?", [$cluster]);
287 :     my @pegs = ();
288 :     while (my $row = $query->Fetch())
289 :     {
290 :     my($peg) = $row->Values(['IsOccurrenceOf(to-link)']);
291 :     push(@pegs,$peg);
292 :     }
293 :     return @pegs;
294 :     }
295 :    
296 :     sub co_occurrence_evidence {
297 :     my($db,$peg1,$peg2) = @_;
298 :    
299 :     my $key = ($peg1 lt $peg2) ? "$peg1:$peg2" : "$peg2:$peg1";
300 :     my $query = $db->Get("Determines PairSet",
301 :     "Determines(from-link) = ?", [$key]);
302 :     my($sc,$pairs);
303 :    
304 :     if ($query)
305 :     {
306 :     my $row = $query->Fetch();
307 :     if ($row)
308 :     {
309 :     my($set) = $row->Values(['PairSet(id)']);
310 :     ($sc,$pairs) = &co_occurrence_set($db,$set);
311 :     my $i;
312 :     for ($i=0; ($i < @$pairs) && ($pairs->[$i]->[0] ne $peg1); $i++) {}
313 :     if ($i == @$pairs)
314 :     {
315 :     $pairs = [map { [$_->[1],$_->[0]] } @$pairs];
316 :     for ($i=0; ($i < @$pairs) && ($pairs->[$i]->[0] ne $peg1); $i++) {}
317 :     if ($i == @$pairs)
318 :     {
319 :     print STDERR &Dumper($pairs,$i,$peg1,$peg2);
320 :     die "Something is screwed up";
321 :     }
322 :     }
323 :     splice(@$pairs,$i,1);
324 :     }
325 :     return $pairs;
326 :     }
327 :     return [];
328 :     }
329 : overbeek 1.3
330 : overbeek 1.8 sub is_co_occurrence_pair {
331 :     my($db,$peg1,$peg2) = @_;
332 :    
333 :     my $key = ($peg1 lt $peg2) ? "$peg1:$peg2" : "$peg2:$peg1";
334 :     my $query = $db->Get("Pairing",
335 :     "Pairing(id) = ?",[$key]);
336 :     return ($query && $query->Fetch()) ? 1 : 0;
337 :     }
338 :    
339 : overbeek 1.3 sub in_pair_set {
340 :     my($db,$peg1,$peg2) = @_;
341 :    
342 :     my($key,$inv,$set);
343 : parrello 1.5 $key = ($peg1 lt $peg2) ? "$peg1:$peg2" : "$peg2:$peg1";
344 : overbeek 1.3
345 :     my $query = $db->Get("Pairing Determines",
346 :     "Pairing(id) = ?",[$key]);
347 :     if ($query)
348 :     {
349 :     my $row = $query->Fetch();
350 :     if ($row)
351 :     {
352 :     my($pair,$inverted) = $row->Values(['Determines(to-link)','Determines(inverted)']);
353 :     if ($peg1 gt $peg2)
354 :     {
355 :     $inverted = ! $inverted;
356 :     }
357 :     return ($pair,$inverted);
358 :     }
359 :     }
360 :     return undef;
361 :     }
362 :    
363 : overbeek 1.4 # This routine inserts a Pairing for $peg1 and $peg2
364 : overbeek 1.3 sub insert_pair {
365 :     my($db,$peg1,$peg2,$set) = @_;
366 : parrello 1.5 # Compute the key and the inversion flag.
367 :     my ($invertFlag, $key);
368 :     if ($peg1 lt $peg2) {
369 :     $key = "$peg1:$peg2";
370 :     $invertFlag = 0;
371 :     } else {
372 :     $key = "$peg2:$peg1";
373 :     $invertFlag = 1;
374 :     }
375 :     # Create the pairing.
376 : overbeek 1.8
377 :     my $query = $db->Get("Pairing",
378 :     "Pairing(id) = ?",[$key]);
379 :     if ((! $query) || (! $query->Fetch()))
380 :     {
381 :     $db->InsertObject('Pairing', id => $key);
382 :     # Connect it to its constituent features.
383 :     $db->InsertObject('IsInPair', from_link => $peg1, to_link => $key);
384 :     $db->InsertObject('IsInPair', from_link => $peg2, to_link => $key);
385 :     }
386 :    
387 :     $query = $db->Get("Determines",
388 :     "Determines(from-link) = ?",[$key]);
389 :     if ($query && $query->Fetch())
390 :     {
391 :     return 0; # it is already connected to another PairSet
392 :     }
393 :    
394 : parrello 1.5 # Insert it into the pair set.
395 :     $db->InsertObject('IsDeterminedBy', from_link => $set, to_link => $key,
396 :     inverted => $invertFlag);
397 : overbeek 1.3 print STDERR "inserted $peg1:$peg2 into $set\n";
398 : overbeek 1.8 return 1;
399 : overbeek 1.3 }
400 :    
401 : overbeek 1.4 # This routine creates a new PairSet, returning the ID
402 : overbeek 1.3 sub new_set {
403 :     my($db) = @_;
404 : parrello 1.5 # Create the new pair set.
405 :     my $retVal = $db->InsertNew('PairSet', score => 0);
406 :     # Return its ID.
407 :     return $retVal;
408 : overbeek 1.3 }
409 : overbeek 1.1
410 : overbeek 1.4
411 :     # This routine deletes a PairSet and relationships to Pairings
412 :     # (but not the Pairings).
413 :    
414 : overbeek 1.3 sub delete_PairSet {
415 :     my($db,$pair_set) = @_;
416 : parrello 1.5 # Disconnect the pair set from all of its pairings.
417 :     $db->Disconnect('IsDeterminedBy', PairSet => $pair_set);
418 :     # Delete the pair set. Because we've already disconnected, the pairings
419 :     # are safe.
420 :     $db->Delete(PairSet => $pair_set);
421 : overbeek 1.3 print STDERR "deleted PairSet $pair_set\n";
422 :     }
423 :    
424 : overbeek 1.4 # This routine deletes the relationship between a Pairing and a PairSet
425 :     # (but does not delete either entity).
426 :    
427 : overbeek 1.3 sub delete_pair_from_PairSet {
428 :     my($db,$pair,$PairSet) = @_;
429 : parrello 1.5 # Delete the relationship connecting this set to this pair.
430 :     $db->DeleteRow('IsDeterminedBy', $PairSet, $pair);
431 : overbeek 1.3 }
432 :    
433 : overbeek 1.4 # This routine updates the "score" field in a PairSet
434 : overbeek 1.3 sub set_pair_set_sc {
435 :     my($db,$pair_set,$sc) = @_;
436 : parrello 1.5 # Update the score in place.
437 :     $db->UpdateEntity(PairSet => $pair_set, score => $sc);
438 : overbeek 1.3 print STDERR "reset score of PairSet to $sc\n";
439 :     }
440 :    
441 : overbeek 1.8 sub check_and_rescore {
442 :     my($db,$pair_set) = @_;
443 :    
444 :     my $sc = &rescore($db,$pair_set);
445 :     if ($sc < 5)
446 :     {
447 :     print STDERR "Rescored PairSet $pair_set to $sc, so we are deleting it\n";
448 :     &delete_PairSet($db,$pair_set);
449 :     return 0;
450 :     }
451 :     return 1;
452 :     }
453 :    
454 : overbeek 1.3 sub rescore {
455 :     my($db,$pair_set) = @_;
456 :    
457 :     my $fig = new FIG;
458 :     my(undef,$pairs) = &co_occurrence_set($db,$pair_set);
459 : overbeek 1.8 defined($pairs) || confess $pair_set;
460 :    
461 :     my $sz = @$pairs;
462 : overbeek 1.3 my %genome_sets;
463 :     foreach my $pair (@$pairs)
464 :     {
465 : overbeek 1.6 $genome_sets{$fig->get_representative_genome(&FIG::genome_of($pair->[0]))} = 1;
466 : overbeek 1.3 }
467 : overbeek 1.8 my $sz_reps = keys(%genome_sets);
468 :     print STDERR "rescoring: $sz pairs in $pair_set, score=$sz_reps\n";
469 : overbeek 1.3 my $sc = keys(%genome_sets);
470 :     &set_pair_set_sc($db,$pair_set,$sc);
471 : overbeek 1.8 return $sc;
472 : overbeek 1.3 }
473 :    
474 :     sub extend_pairs_and_pair_sets {
475 :     my($db,$pchF) = @_;
476 :    
477 :     open(PCHS,"sort $pchF |") || die "could not open $pchF";
478 :    
479 :     my $line = <PCHS>;
480 :     while ($line && ($line =~ /^((\S+)\t(\S+))\t(\S+)\t(\S+)/))
481 :     {
482 :     my $curr = $1;
483 : overbeek 1.7 print STDERR "PAIR: $curr\n";
484 : overbeek 1.3 my $set = [];
485 :     while ($line && ($line =~ /^((\S+)\t(\S+))\t(\S+)\t(\S+)/) && ($1 eq $curr))
486 :     {
487 :     push(@$set,[$2,$3,$4,$5]);
488 :     $line = <PCHS>;
489 :     }
490 : overbeek 1.8 print STDERR "processing set ",&Dumper($set);
491 : overbeek 1.3 &add_pch_set($db,$set);
492 :     }
493 :     close(PCHS);
494 :     }
495 :    
496 :     sub add_pch_set {
497 :     my($db,$set) = @_;
498 :    
499 :     my $fig = new FIG;
500 :     my @pairs = ([$set->[0]->[0],$set->[0]->[1]],map { [$_->[2],$_->[3]] } @$set);
501 :    
502 : overbeek 1.8 my($pair,@unplaced,%in,%members);
503 :     my $ok = 1;
504 : overbeek 1.3 foreach $pair (@pairs)
505 :     {
506 :     my($pair_set,$inverted) = &in_pair_set($db,@$pair);
507 : overbeek 1.8 $inverted = $inverted ? 1 : 0; # force to {0,1} for boolean
508 : overbeek 1.3 if (! defined($pair_set))
509 :     {
510 :     push(@unplaced,$pair);
511 :     }
512 :     else
513 :     {
514 : overbeek 1.7 print STDERR "IN: ",join(":",@$pair)," SET: $pair_set INVERTED: $inverted\n";
515 : overbeek 1.8 my($keyPI) = "$pair_set:$inverted";
516 :     $in{"$keyPI"}++;
517 :     push(@{$members{$keyPI}},join(",",@$pair));
518 :    
519 :     my $opposite = $inverted ? 0 : 1;
520 :     if ($in{"$pair_set:$opposite"})
521 :     {
522 :     # We have a case in which a set contains both peg1:peg2 and peg2:peg1
523 :     print STDERR "We have a set containing peg1:peg2 and peg2:peg1 ", &Dumper(\%members,$keyPI,$pair);
524 :     $ok = 0;
525 :     }
526 : overbeek 1.3 }
527 :     }
528 : overbeek 1.8
529 :     if (! $ok)
530 :     {
531 :     print STDERR "BAILING OUT\n";
532 :     return;
533 :     }
534 :    
535 :     print STDERR "characterized pairs ",&Dumper(\@unplaced,\%in);
536 : overbeek 1.7 my @in_sets = map { [split(/:/,$_)] } sort { $in{$b} <=> $in{$a} } keys(%in);
537 :     print STDERR "in_sets: ", &Dumper(\@in_sets);
538 :    
539 : overbeek 1.3 if (@in_sets > 1)
540 :     {
541 :     my $i;
542 :     for ($i=1; ($i < @in_sets); $i++)
543 :     {
544 : overbeek 1.7 my(undef,$in_old) = &co_occurrence_set($db,$in_sets[$i]->[0]);
545 : overbeek 1.3 if ($in_sets[$i]->[1] ne $in_sets[0]->[1])
546 :     {
547 :     $in_old = [ map { [$_->[1],$_->[0]] } @$in_old];
548 :     }
549 : overbeek 1.7 print STDERR "copying pairs from $in_sets[$i]->[0] to $in_sets[0]->[0]\n";
550 : overbeek 1.3 foreach $pair (@$in_old)
551 :     {
552 : overbeek 1.7 print STDERR "\tInserting pair ",join(":",@$pair)," into $in_sets[0]->[0]\n";
553 : overbeek 1.8 &insert_pair($db,@$pair,$in_sets[0]->[0]);
554 : overbeek 1.3 }
555 : overbeek 1.7 print STDERR "Deleting PairSet $in_sets[$i]->[0]\n";
556 : overbeek 1.8 &delete_PairSet($db,$in_sets[$i]->[0]);
557 : overbeek 1.3 }
558 : overbeek 1.7 print STDERR "rescore PairSet $in_sets[0]->[0]\n";
559 : overbeek 1.8 &check_and_rescore($db,$in_sets[0]->[0]);
560 : overbeek 1.3 }
561 :    
562 :     if (@in_sets == 1)
563 :     {
564 : overbeek 1.8 print STDERR "just inserting unplaced\n";
565 :     if (@unplaced > 0)
566 : overbeek 1.3 {
567 : overbeek 1.8 foreach $pair (@unplaced)
568 :     {
569 :     print STDERR &Dumper(['inserting pair',$pair,$in_sets[0]]);
570 :     my($peg3,$peg4) = $in_sets[0]->[1] ? ($pair->[1],$pair->[0]) : @$pair;
571 :     print STDERR "Inserting Pairing $peg3:$peg4\n";
572 :     &insert_pair($db,$peg3,$peg4,$in_sets[0]->[0]);
573 :     }
574 :     print STDERR "rescore PairSet $in_sets[0]->[0]\n";
575 :     &check_and_rescore($db,$in_sets[0]->[0]);
576 : overbeek 1.3 }
577 :     }
578 :     elsif (@unplaced >= 5)
579 :     {
580 : overbeek 1.8 print STDERR "No existing sets, maybe we will create one", &Dumper(['unplaced',\@unplaced]);
581 : overbeek 1.7
582 : overbeek 1.3 my %genome_sets;
583 :     foreach my $pair (@unplaced)
584 :     {
585 : overbeek 1.6 $genome_sets{$fig->get_representative_genome(&FIG::genome_of($pair->[0]))} = 1;
586 : overbeek 1.3 }
587 :     my $sc = keys(%genome_sets);
588 :     if ($sc >= 5)
589 :     {
590 : overbeek 1.8 print STDERR "sc=$sc, so we will create one\n";
591 : overbeek 1.3 my $new_set = &new_set($db);
592 : overbeek 1.8 print STDERR "new PairSet=$new_set\n";
593 :    
594 :     if (! $new_set)
595 :     {
596 :     confess "Failed to acquire a new PairSet\n";
597 :     }
598 :     else
599 :     {
600 :     print STDERR "got new PairSet=$new_set\n",&Dumper(&co_occurrence_set($db,$new_set));
601 :     }
602 : overbeek 1.3 foreach my $pair (@unplaced)
603 :     {
604 : overbeek 1.7 print STDERR "Inserting pair ",join(":",@$pair)," into $new_set\n";
605 : overbeek 1.8 my $rc = &insert_pair($db,@$pair,$new_set);
606 :     if (! $rc) { print STDERR "\tinsertion failed: Pairing already exists\n"; }
607 : overbeek 1.3 }
608 : overbeek 1.8 print STDERR "rescoring new set $new_set\n";
609 :     &check_and_rescore($db,$new_set);
610 : overbeek 1.3 }
611 :     }
612 :     }
613 :    
614 :     sub cleanup_pair_sets {
615 :     my($db) = @_;
616 :    
617 : overbeek 1.8 my $fig = new FIG;
618 : overbeek 1.3 foreach my $pair_set (&all_co_occurrence_pair_sets)
619 :     {
620 : overbeek 1.7 print STDERR "cleaning up $pair_set\n";
621 : overbeek 1.3 my(undef,$set) = &co_occurrence_set($db,$pair_set);
622 : overbeek 1.8 next if (! $set);
623 :     my @partitioned = sort { @$b <=> @$a } &partition_pair_set($fig,$db,$set);
624 : overbeek 1.7
625 :     my $N = @partitioned;
626 :     print STDERR "$pair_set splits into $N partitions\n";
627 :    
628 : overbeek 1.3 my $keep = shift @partitioned;
629 : overbeek 1.8
630 : overbeek 1.7 if ((! $keep) || (@$keep < 5))
631 : overbeek 1.3 {
632 : overbeek 1.7 print STDERR "\tbut none of them are big enough, so they all get pitched\n";
633 : overbeek 1.3 &delete_PairSet($db,$pair_set);
634 :     }
635 :     else
636 :     {
637 : overbeek 1.8 my $sz = @$keep; print STDERR "\tkeeping partition of size=$sz\n";
638 : overbeek 1.3 my %to_keep = map { join(":",@$_) => 1 } @$keep;
639 :     foreach my $pair (@$set)
640 :     {
641 : overbeek 1.7 my $key = join(":",@$pair);
642 :     if (! $to_keep{ $key } )
643 : overbeek 1.3 {
644 : overbeek 1.7 print STDERR "\tdeleting $key from $pair_set\n";
645 : overbeek 1.3 &delete_pair_from_PairSet($db,$pair,$pair_set);
646 :     }
647 :     }
648 : overbeek 1.8 if (@$keep < @$set) { &check_and_rescore($db,$pair_set) }
649 :    
650 : overbeek 1.3 foreach my $new_set (@partitioned)
651 :     {
652 : overbeek 1.7 my $sz = @$new_set;
653 :     print STDERR "\tconsidering a partition of size $sz\n";
654 : overbeek 1.3 if (@$new_set >= 5)
655 :     {
656 :     my $next_set = &new_set($db);
657 : overbeek 1.7 print STDERR "\tcreating new PairSet $next_set\n";
658 : overbeek 1.3 foreach my $pair1 (@$new_set)
659 :     {
660 : overbeek 1.7 my $key = join(":",@$pair1);
661 :     print STDERR "\t\tadding $key\n";
662 : overbeek 1.3 &insert_pair($db,@$pair1,$next_set);
663 :     }
664 : overbeek 1.8 &check_and_rescore($db,$next_set);
665 : overbeek 1.3 }
666 :     }
667 :     }
668 :     }
669 :     }
670 :    
671 :     sub partition_pair_set {
672 : overbeek 1.8 my($fig,$db,$set) = @_;
673 : overbeek 1.3
674 :     my @split = ();
675 : overbeek 1.8 my @sets = &part1($fig,$set,0);
676 : overbeek 1.3 foreach my $part1 (@sets)
677 :     {
678 : overbeek 1.8 foreach my $part2 (&part1($fig,$part1,1))
679 : overbeek 1.3 {
680 :     push(@split,$part2);
681 :     }
682 :     }
683 :     return sort { @$b <=> @$a } @split;
684 :     }
685 :    
686 :     sub part1 {
687 : overbeek 1.8 my($fig,$set,$index) = @_;
688 : overbeek 1.3
689 :     my %entries = map { $_->[$index] => $_ } @$set;
690 : overbeek 1.8 my $ids = [grep { $fig->is_real_feature($_) } map { $_->[$index] } @$set];
691 :     my @partition = &PartitionSeqs::partition({ pegs => $ids, use => 'blast', identity => 0.4, coverage => 0.7 });
692 : overbeek 1.3 return map { [map { $entries{$_} } @$_] } @partition;
693 :     }
694 :    
695 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3