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

Annotation of /FigKernelPackages/FC.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3