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

Annotation of /FigKernelPackages/FC.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3