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

Annotation of /FigKernelPackages/FC.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (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.3 use FIG;
24 :     sub PartitionSeqs;
25 : overbeek 1.1
26 :     sub co_occurs {
27 :     my($db,$peg) = @_;
28 :    
29 :     my @tuples;
30 :     my $query = $db->Get("IsInPair Pairing Determines PairSet",
31 :     "IsInPair(from-link) = ?", [$peg]);
32 :     while (my $row = $query->Fetch())
33 :     {
34 :     my ($pair, $score) = $row->Values(['Determines(from-link)','PairSet(score)']);
35 :     my ($fid) = grep { $_ ne $peg } split(/:/, $pair);
36 :     push @tuples, [$fid, $score];
37 :     }
38 :     return @tuples;
39 :     }
40 :    
41 :     sub co_occurrence_and_evidence {
42 :     my($db,$peg) = @_;
43 :    
44 :     my @tuples;
45 :     my $query = $db->Get("IsInPair Pairing Determines PairSet IsDeterminedBy",
46 :     "IsInPair(from-link) = ?", [$peg]);
47 :     my @tuoples = ();
48 :     while (my $row = $query->Fetch())
49 :     {
50 :     push(@tuples,[$row->Values(['Determines(from-link)',
51 :     'PairSet(score)',
52 : overbeek 1.3 'IsDeterminedBy(inverted)',
53 : overbeek 1.1 'IsDeterminedBy(to-link)'])]);
54 :     }
55 : overbeek 1.2 @tuples = sort { $a->[0] <=> $b->[0] } @tuples;
56 :     my @co_occurs = ();
57 :    
58 :     my $x = shift @tuples;
59 :     while ($x)
60 :     {
61 :     my $curr = $x->[0];
62 :     my $sc = $x->[1];
63 :     my @set = ();
64 :     while ($x && ($x->[0] eq $curr))
65 :     {
66 :     my($peg3,$peg4) = split(/:/,$x->[3]);
67 :     push(@set,($x->[2]) ? [$peg3,$peg4] : [$peg4,$peg3]);
68 :     $x = shift @tuples;
69 :     }
70 :     my $i;
71 :     for ($i=0; ($i < @set) && ($peg ne $set[$i]->[0]); $i++) {}
72 :     if ($i == @set)
73 :     {
74 :     @set = map { [$_->[1],$_->[0]] } @set;
75 :     }
76 :     my($peg2) = grep { $_ ne $peg } split(/:/,$curr);
77 : overbeek 1.3 push(@co_occurs,[$sc,$peg2,[grep { $_->[0] ne $peg } @set]]);
78 : overbeek 1.2 }
79 :     return sort { $b->[0] <=> $a->[0] } @co_occurs;
80 : overbeek 1.1 }
81 :    
82 :     sub co_occurrence_set {
83 :     my($db,$set) = @_;
84 :    
85 :     my $query = $db->Get("PairSet IsDeterminedBy",
86 :     "PairSet(id) = ?", [$set]);
87 :     my $pairs;
88 :     my $sc;
89 :     while (my $row = $query->Fetch())
90 :     {
91 :     if (! $sc) { ($sc) = $row->Values(['PairSet(score)']); }
92 :     my($pairing,$inverted) = $row->Values(['IsDeterminedBy(to-link)','IsDeterminedBy(inverted)']);
93 :     my($peg1,$peg2) = split(/:/,$pairing);
94 :     push(@$pairs,$inverted ? [$peg2,$peg1] : [$peg1,$peg2]);
95 :     }
96 :     return ($sc,$pairs);
97 :     }
98 :    
99 : overbeek 1.3 sub all_co_occurrence_pair_sets {
100 :     my($db) = @_;
101 :    
102 :     my $query = $db->Get("PairSet","",[]);
103 :     my @PairSets = ();
104 :     while (my $row = $query->Fetch())
105 :     {
106 :     my($pair_set) = $row->Values(['PairSet(id)']);
107 :     push(@PairSets,$pair_set);
108 :     }
109 :     return @PairSets;
110 :     }
111 :    
112 :     sub all_co_occurrence_clusters {
113 :     my($db) = @_;
114 :    
115 :     my $query = $db->Get("Cluster","",[]);
116 :     my @Clusters = ();
117 :     while (my $row = $query->Fetch())
118 :     {
119 :     my($id) = $row->Values(['Cluster(id)']);
120 :     push(@Clusters,$id);
121 :     }
122 :     return @Clusters;
123 :     }
124 :    
125 : overbeek 1.1 sub co_occurrence_cluster {
126 :     my($db,$cluster) = @_;
127 :    
128 :     my $query = $db->Get("Cluster IsOccurrenceOf",
129 :     "Cluster(id) = ?", [$cluster]);
130 :     my @pegs = ();
131 :     while (my $row = $query->Fetch())
132 :     {
133 :     my($peg) = $row->Values(['IsOccurrenceOf(to-link)']);
134 :     push(@pegs,$peg);
135 :     }
136 :     return @pegs;
137 :     }
138 :    
139 :     sub co_occurrence_evidence {
140 :     my($db,$peg1,$peg2) = @_;
141 :    
142 :     my $key = ($peg1 lt $peg2) ? "$peg1:$peg2" : "$peg2:$peg1";
143 :     my $query = $db->Get("Determines PairSet",
144 :     "Determines(from-link) = ?", [$key]);
145 :     my($sc,$pairs);
146 :    
147 :     if ($query)
148 :     {
149 :     my $row = $query->Fetch();
150 :     if ($row)
151 :     {
152 :     my($set) = $row->Values(['PairSet(id)']);
153 :     ($sc,$pairs) = &co_occurrence_set($db,$set);
154 :     my $i;
155 :     for ($i=0; ($i < @$pairs) && ($pairs->[$i]->[0] ne $peg1); $i++) {}
156 :     if ($i == @$pairs)
157 :     {
158 :     $pairs = [map { [$_->[1],$_->[0]] } @$pairs];
159 :     for ($i=0; ($i < @$pairs) && ($pairs->[$i]->[0] ne $peg1); $i++) {}
160 :     if ($i == @$pairs)
161 :     {
162 :     print STDERR &Dumper($pairs,$i,$peg1,$peg2);
163 :     die "Something is screwed up";
164 :     }
165 :     }
166 :     splice(@$pairs,$i,1);
167 :     }
168 :     return $pairs;
169 :     }
170 :     return [];
171 :     }
172 : overbeek 1.3
173 :     sub in_pair_set {
174 :     my($db,$peg1,$peg2) = @_;
175 :    
176 :     my($key,$inv,$set);
177 : parrello 1.5 $key = ($peg1 lt $peg2) ? "$peg1:$peg2" : "$peg2:$peg1";
178 : overbeek 1.3
179 :     my $query = $db->Get("Pairing Determines",
180 :     "Pairing(id) = ?",[$key]);
181 :     if ($query)
182 :     {
183 :     my $row = $query->Fetch();
184 :     if ($row)
185 :     {
186 :     my($pair,$inverted) = $row->Values(['Determines(to-link)','Determines(inverted)']);
187 :     if ($peg1 gt $peg2)
188 :     {
189 :     $inverted = ! $inverted;
190 :     }
191 :     return ($pair,$inverted);
192 :     }
193 :     }
194 :     return undef;
195 :     }
196 :    
197 : overbeek 1.4 # This routine inserts a Pairing for $peg1 and $peg2
198 : overbeek 1.3 sub insert_pair {
199 :     my($db,$peg1,$peg2,$set) = @_;
200 : parrello 1.5 # Compute the key and the inversion flag.
201 :     my ($invertFlag, $key);
202 :     if ($peg1 lt $peg2) {
203 :     $key = "$peg1:$peg2";
204 :     $invertFlag = 0;
205 :     } else {
206 :     $key = "$peg2:$peg1";
207 :     $invertFlag = 1;
208 :     }
209 :     # Create the pairing.
210 :     $db->InsertObject('Pairing', id => $key);
211 :     # Connect it to its constituent features.
212 :     $db->InsertObject('IsInPair', from_link => $peg1, to_link => $key);
213 :     $db->InsertObject('IsInPair', from_link => $peg2, to_link => $key);
214 :     # Insert it into the pair set.
215 :     $db->InsertObject('IsDeterminedBy', from_link => $set, to_link => $key,
216 :     inverted => $invertFlag);
217 : overbeek 1.3 print STDERR "inserted $peg1:$peg2 into $set\n";
218 :     }
219 :    
220 : overbeek 1.4 # This routine creates a new PairSet, returning the ID
221 : overbeek 1.3 sub new_set {
222 :     my($db) = @_;
223 : parrello 1.5 # Create the new pair set.
224 :     my $retVal = $db->InsertNew('PairSet', score => 0);
225 :     # Return its ID.
226 :     return $retVal;
227 : overbeek 1.3 }
228 : overbeek 1.1
229 : overbeek 1.4
230 :     # This routine deletes a PairSet and relationships to Pairings
231 :     # (but not the Pairings).
232 :    
233 : overbeek 1.3 sub delete_PairSet {
234 :     my($db,$pair_set) = @_;
235 : parrello 1.5 # Disconnect the pair set from all of its pairings.
236 :     $db->Disconnect('IsDeterminedBy', PairSet => $pair_set);
237 :     # Delete the pair set. Because we've already disconnected, the pairings
238 :     # are safe.
239 :     $db->Delete(PairSet => $pair_set);
240 : overbeek 1.3 print STDERR "deleted PairSet $pair_set\n";
241 :     }
242 :    
243 : overbeek 1.4 # This routine deletes the relationship between a Pairing and a PairSet
244 :     # (but does not delete either entity).
245 :    
246 : overbeek 1.3 sub delete_pair_from_PairSet {
247 :     my($db,$pair,$PairSet) = @_;
248 : parrello 1.5 # Delete the relationship connecting this set to this pair.
249 :     $db->DeleteRow('IsDeterminedBy', $PairSet, $pair);
250 : overbeek 1.3 }
251 :    
252 : overbeek 1.4 # This routine updates the "score" field in a PairSet
253 : overbeek 1.3 sub set_pair_set_sc {
254 :     my($db,$pair_set,$sc) = @_;
255 : parrello 1.5 # Update the score in place.
256 :     $db->UpdateEntity(PairSet => $pair_set, score => $sc);
257 : overbeek 1.3 print STDERR "reset score of PairSet to $sc\n";
258 :     }
259 :    
260 :     sub rescore {
261 :     my($db,$pair_set) = @_;
262 :    
263 :     my $fig = new FIG;
264 :     my(undef,$pairs) = &co_occurrence_set($db,$pair_set);
265 :     my %genome_sets;
266 :     foreach my $pair (@$pairs)
267 :     {
268 : overbeek 1.6 $genome_sets{$fig->get_representative_genome(&FIG::genome_of($pair->[0]))} = 1;
269 : overbeek 1.3 }
270 :     my $sc = keys(%genome_sets);
271 :     &set_pair_set_sc($db,$pair_set,$sc);
272 :     }
273 :    
274 :     sub extend_pairs_and_pair_sets {
275 :     my($db,$pchF) = @_;
276 :    
277 :     die "NOT DEBUGGED: RAO March 8";
278 :    
279 :     open(PCHS,"sort $pchF |") || die "could not open $pchF";
280 :    
281 :     my $line = <PCHS>;
282 :     while ($line && ($line =~ /^((\S+)\t(\S+))\t(\S+)\t(\S+)/))
283 :     {
284 :     my $curr = $1;
285 :     my $set = [];
286 :     while ($line && ($line =~ /^((\S+)\t(\S+))\t(\S+)\t(\S+)/) && ($1 eq $curr))
287 :     {
288 :     push(@$set,[$2,$3,$4,$5]);
289 :     $line = <PCHS>;
290 :     }
291 :     &add_pch_set($db,$set);
292 :     }
293 :     close(PCHS);
294 :     }
295 :    
296 :     sub add_pch_set {
297 :     my($db,$set) = @_;
298 :    
299 :     my $fig = new FIG;
300 :     my @pairs = ([$set->[0]->[0],$set->[0]->[1]],map { [$_->[2],$_->[3]] } @$set);
301 :    
302 :     my($pair,@unplaced,%in);
303 :     foreach $pair (@pairs)
304 :     {
305 :     my($pair_set,$inverted) = &in_pair_set($db,@$pair);
306 :     if (! defined($pair_set))
307 :     {
308 :     push(@unplaced,$pair);
309 :     }
310 :     else
311 :     {
312 :     $in{"$pair_set:$inverted"}++;
313 :     }
314 :     }
315 :     my @in_sets = map { [split(/:/,$_)] } sort { $b->{$_} <=> $a->{$_} } keys(%in);
316 :     if (@in_sets > 1)
317 :     {
318 :     my $i;
319 :     for ($i=1; ($i < @in_sets); $i++)
320 :     {
321 :     my(undef,$in_old) = &co_occurence_set($db,$in_sets[$i]->[0]);
322 :     if ($in_sets[$i]->[1] ne $in_sets[0]->[1])
323 :     {
324 :     $in_old = [ map { [$_->[1],$_->[0]] } @$in_old];
325 :     }
326 :     foreach $pair (@$in_old)
327 :     {
328 :     &insert_pair($db,@$pair,$in_sets[0]->[0]);
329 :     }
330 :     &delete_set($db,$in_sets[$i]->[0]);
331 :     }
332 :     }
333 :    
334 :     if (@in_sets == 1)
335 :     {
336 :     foreach $pair (@unplaced)
337 :     {
338 :     my($peg3,$peg4) = $in_sets[0]->[0] ? ($pair->[1],$pair->[0]) : @$pair;
339 :     &insert_pair($db,$peg3,$peg4,$in_sets[0]->[0]);
340 :     }
341 :     }
342 :     elsif (@unplaced >= 5)
343 :     {
344 :     my %genome_sets;
345 :     foreach my $pair (@unplaced)
346 :     {
347 : overbeek 1.6 $genome_sets{$fig->get_representative_genome(&FIG::genome_of($pair->[0]))} = 1;
348 : overbeek 1.3 }
349 :     my $sc = keys(%genome_sets);
350 :     if ($sc >= 5)
351 :     {
352 :     my $new_set = &new_set($db);
353 :     foreach my $pair (@unplaced)
354 :     {
355 :     &insert_pair($db,@$pair,$new_set);
356 :     }
357 :     }
358 :     }
359 :     }
360 :    
361 :     sub cleanup_pair_sets {
362 :     my($db) = @_;
363 :    
364 :     die "NOT DEBUGGED: RAO March 8";
365 :    
366 :     foreach my $pair_set (&all_co_occurrence_pair_sets)
367 :     {
368 :     my(undef,$set) = &co_occurrence_set($db,$pair_set);
369 :     my @partitioned = sort { @$b <=> @$a } &partition_pair_set($db,$set);
370 :     my $keep = shift @partitioned;
371 :     if (@$keep < 5)
372 :     {
373 :     &delete_PairSet($db,$pair_set);
374 :     }
375 :     else
376 :     {
377 :     my %to_keep = map { join(":",@$_) => 1 } @$keep;
378 :     foreach my $pair (@$set)
379 :     {
380 :     if (! $to_keep{ join(":",@$pair) } )
381 :     {
382 :     &delete_pair_from_PairSet($db,$pair,$pair_set);
383 :     }
384 :     }
385 :    
386 :     foreach my $new_set (@partitioned)
387 :     {
388 :     if (@$new_set >= 5)
389 :     {
390 :     my $next_set = &new_set($db);
391 :     foreach my $pair1 (@$new_set)
392 :     {
393 :     &insert_pair($db,@$pair1,$next_set);
394 :     }
395 :     }
396 :     }
397 :     }
398 :     }
399 :     }
400 :    
401 :     sub partition_pair_set {
402 :     my($db,$set) = @_;
403 :    
404 :     my @split = ();
405 :     my @sets = &part1($set,0);
406 :     foreach my $part1 (@sets)
407 :     {
408 :     foreach my $part2 (&part1($part1,1))
409 :     {
410 :     push(@split,$part2);
411 :     }
412 :     }
413 :     return sort { @$b <=> @$a } @split;
414 :     }
415 :    
416 :     sub part1 {
417 :     my($set,$index) = @_;
418 :    
419 :     my %entries = map { $_->[$index] => $_ } @$set;
420 :     my $ids = [map { $_->[$index] } @$set];
421 :     my @partition = &PartitionSeqs::partition({ pegs => $ids, use => 'blast' });
422 :     my @tmp = map { [map { $entries{$_} } @$_] } @partition;
423 :     return map { [map { $entries{$_} } @$_] } @partition;
424 :     }
425 :    
426 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3