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

Annotation of /FigKernelPackages/SeedSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 :     package SeedSims;
19 :    
20 :     use FIG;
21 :     use Data::Dumper;
22 :    
23 :     sub load_seed_sims {
24 :    
25 :     my $dir = "$FIG_Config::data/SimPart";
26 :    
27 :     if (! -d "$dir")
28 :     {
29 :     system "$FIG_Config::bin/partition_prots 20000 2> /dev/null | make_sim_part_data $dir";
30 :     }
31 :    
32 :     ((-s "$dir/SimilarityPartitions") &&
33 :     (-s "$dir/ToSimilarityPartition") &&
34 :     (-d "$dir/Partitions"))
35 :     || die "$dir did not get built properly";
36 :    
37 :     my $fig = new FIG;
38 :     my $dbf = $fig->db_handle;
39 :    
40 :     my @tmp = `tail -n 1 $dir/ToSimilarityPartition`;
41 :    
42 :     my $rebuild = 1;
43 :     if ($fig->table_exists('ToSimilarityPartition') &&
44 :     (@tmp == 1) && ($tmp[0] =~ /^(\d+)\t(\S+)/))
45 :     {
46 :     my $part = $1;
47 :     my $md5 = $2;
48 :     my $db_loc = $dbf->SQL("SELECT partition FROM ToSimilarityPartition WHERE md5_hash = '$md5'");
49 :     my $x = $db_loc->[0]->[0];
50 :     $rebuild = ($x ne $part);
51 :     }
52 :    
53 :     if ($rebuild)
54 :     {
55 :     print STDERR "loading similarity partition data\n";
56 :    
57 :     $dbf->drop_table( tbl => "fc_scores" );
58 :     $dbf->create_table( tbl => "fc_scores",
59 :     flds => "partition1 integer, partition2 integer, score integer"
60 :     );
61 :     $dbf->load_table( tbl => "fc_scores",
62 :     file => "$dir/PartitionFC"
63 :     );
64 :    
65 :     $dbf->create_index( idx => "fc_scores_peg1_ix",
66 :     tbl => "fc_scores",
67 :     type => "btree",
68 :     flds => "partition1" );
69 :    
70 :     $dbf->create_index( idx => "fc_scores_sc_ix",
71 :     tbl => "fc_scores",
72 :     type => "btree",
73 :     flds => "score" );
74 :    
75 :     $dbf->drop_table( tbl => "SimilarityPartitions" );
76 :    
77 :     $dbf->create_table( tbl => "SimilarityPartitions",
78 :     flds => "partition integer, md5_hash varchar(32)"
79 :     );
80 :     $dbf->load_table( tbl => "SimilarityPartitions",
81 :     file => "$dir/SimilarityPartitions"
82 :     );
83 :    
84 :     $dbf->create_index( idx => "partition_ix",
85 :     tbl => "SimilarityPartitions",
86 :     type => "btree",
87 :     flds => "partition" );
88 :    
89 :     $dbf->create_index( idx => "hash_ix",
90 :     tbl => "SimilarityPartitions",
91 :     type => "btree",
92 :     flds => "md5_hash" );
93 :    
94 :    
95 :     $dbf->drop_table( tbl => "ToSimilarityPartition" );
96 :     $dbf->create_table( tbl => "ToSimilarityPartition",
97 :     flds => "partition integer, md5_hash varchar(32)"
98 :     );
99 :     $dbf->load_table( tbl => "ToSimilarityPartition",
100 :     file => "$dir/ToSimilarityPartition"
101 :     );
102 :    
103 :     $dbf->create_index( idx => "partition_ix",
104 :     tbl => "ToSimilarityPartition",
105 :     type => "btree",
106 :     flds => "partition" );
107 :    
108 :     $dbf->create_index( idx => "hash_ix",
109 :     tbl => "ToSimilarityPartition",
110 :     type => "btree",
111 :     flds => "md5_hash" );
112 :    
113 :     $dbf->drop_table( tbl => "PartitionReps" );
114 :     $dbf->create_table( tbl => "PartitionReps",
115 :     flds => "partition integer, md5_hash1 varchar(32), md5_hash2 varchar(32)"
116 :     );
117 :     $dbf->load_table( tbl => "PartitionReps",
118 :     file => "$dir/representatives"
119 :     );
120 :    
121 :     $dbf->create_index( idx => "representatives_partition_ix",
122 :     tbl => "PartitionReps",
123 :     type => "btree",
124 :     flds => "partition" );
125 :    
126 :     $dbf->create_index( idx => "representatives_md5_ix",
127 :     tbl => "PartitionReps",
128 :     type => "btree",
129 :     flds => "md5_hash1" );
130 :    
131 :     print STDERR "finished loading similarity partition data\n";
132 :     }
133 :     }
134 :    
135 :     sub representative_of {
136 :     my($fig,$partition,$md5) = @_;
137 :    
138 :     my $dbf = $fig->db_handle;
139 :     my $rep = $dbf->SQL("SELECT md5_hash2 FROM ToSimilarityPartition,PartitionReps WHERE
140 :     (PartitionReps.md5_hash1 = '$md5') AND
141 :     (ToSimilarityPartition.md5_hash = '$md5') AND
142 :     (ToSimilarityPartition.partition = $partition)");
143 :    
144 :     return ($rep && (@$rep == 1)) ? $rep->[0]->[0] : undef;
145 :     }
146 :    
147 :     sub seed_sims {
148 :     my($peg,$parms) = @_;
149 :    
150 :     &load_seed_sims;
151 :     my $fig = new FIG;
152 :     my $md5 = $fig->md5_of_peg($peg);
153 :     # print STDERR "before md5 sims\n";
154 :     my @simsH = &md5_sims($md5,$parms);
155 :     # print STDERR "after md5 sims\n";
156 :     my $sim;
157 :     my @sims = ();
158 :     foreach $sim (@simsH)
159 :     {
160 :     $sim->[0] = $peg;
161 :     my $md5 = $sim->[1];
162 :     my @pegs = $fig->pegs_with_md5($md5);
163 :     foreach my $peg1 (@pegs)
164 :     {
165 :     my $sim1 = bless([@$sim],'Sim');
166 :     if ($peg1 ne $peg)
167 :     {
168 :     $sim1->[1] = $peg1;
169 :     push(@sims,$sim1);
170 :     }
171 :     }
172 :     }
173 :     return @sims;
174 :     }
175 :    
176 :     sub md5_sims {
177 :     my($md5,$parms) = @_;
178 :    
179 :     my @sims = ();
180 :     my $fig = new FIG;
181 :     my $dbf = $fig->db_handle;
182 :     if (! $md5)
183 :     {
184 :     print STDERR "MISSING MD5 HASH for $peg\n";
185 :     return ();
186 :     }
187 :     my $seq1 = &seq_with_md5($fig,$md5);
188 :     my @blastout;
189 :     if ($seq1)
190 :     {
191 :     my $ln1 = length($seq1);
192 :     my $tmpF = "$FIG_Config::temp/md5$$.fasta";
193 :     open(TMP,">$tmpF") || die "could not open $tmpF";
194 :     print TMP ">$md5\n$seq1\n";
195 :     close(TMP);
196 :    
197 :     my $db_loc = $dbf->SQL("SELECT partition FROM ToSimilarityPartition WHERE md5_hash = '$md5'");
198 :     my $x = $db_loc->[0]->[0];
199 :     if ($x)
200 :     {
201 :     my $dir1 = $x % 1000;
202 :     my $partF = "$FIG_Config::data/SimPart/Partitions/$dir1/$x";
203 :     if (! ((-s "$partF.psq") && ((-M $partF) > (-M "$partF.psq"))))
204 :     {
205 :     system "$FIG_Config::ext_bin/formatdb -i $partF -p T";
206 :     }
207 :    
208 :     # print STDERR "before blast\n";
209 :     @blastout = map { chop; bless([split(/\s+/,$_),$ln1],'Sim') }
210 :     `blastall -i $tmpF -d $partF -p blastp -m 8 -e 0.01 $parms`;
211 :     my %ln2 = map { $_ =~ /^(\S+)\t(\S+)/; $1 => $2 } `cut -f1,2 $partF.summary`;
212 :    
213 :     my $sim;
214 :     foreach $sim (@blastout)
215 :     {
216 :     my $ln2 = $ln2{$sim->id2};
217 :     if ($ln2)
218 :     {
219 :     push(@$sim,$ln2);
220 :     push(@sims,$sim);
221 :     }
222 :     }
223 :     }
224 :     unlink $tmpF;
225 :     }
226 :     # return @sims;
227 :     return @blastout;
228 :     }
229 :    
230 :     sub seq_with_md5 {
231 :     my($fig,$md5) = @_;
232 :    
233 :     my @pegs = $fig->pegs_with_md5($md5);
234 :    
235 :     my($i,$x);
236 :     for ($i=0; ($i < @pegs) && (! ($x = $fig->get_translation($pegs[$i]))); $i++) {}
237 :     if ($i == @pegs)
238 :     {
239 :     return undef;
240 :     }
241 :     $x =~ s/\*$//;
242 :     return $x;
243 :     }
244 :    
245 :     sub length_of_seq_with_md5 {
246 :     my($fig,$md5) = @_;
247 :    
248 :     my @pegs = $fig->pegs_with_md5($md5);
249 :    
250 :     my($i,$x);
251 :     for ($i=0; ($i < @pegs) && (! ($x = $fig->translation_length($pegs[$i]))); $i++) {}
252 :     if ($i == @pegs)
253 :     {
254 :     return undef;
255 :     }
256 :     return $x;
257 :     }
258 :    
259 :     sub md5_to_sim_set {
260 :     my($md5) = @_;
261 :    
262 :     my $fig = new FIG;
263 :     my $dbf = $fig->db_handle;
264 :    
265 :     my $db_loc = $dbf->SQL("SELECT partition FROM ToSimilarityPartition WHERE md5_hash = '$md5'");
266 :     my $x = defined($db_loc->[0]->[0]) ? $db_loc->[0]->[0] : undef;
267 :     }
268 :    
269 :     sub partition_coupling_evidence {
270 :     my($fig,$ss1,$ss2) = @_;
271 :    
272 :     if ($ss1 && $ss2)
273 :     {
274 :     my $in1 = &load_locs($ss1);
275 :     my $in2 = &load_locs($ss2);
276 :     my @pairs = ();
277 :     foreach my $genome (sort keys(%$in1))
278 :     {
279 :     my($xH1,$xH2);
280 :     if ($xH1 = $in1->{$genome})
281 :     {
282 :     $xH2 = $in2->{$genome};
283 :    
284 :     my($contig,$YL1,$YL2);
285 :     foreach $contig (keys(%$xH1))
286 :     {
287 :     if ($YL2 = $xH2->{$contig})
288 :     {
289 :     $YL1 = $xH1->{$contig};
290 :     my $i1;
291 :     my $i2;
292 :     for ($i1 =0; ($i1 < @$YL1); $i1++)
293 :     {
294 :     for ($i2=0; ($i2 < @$YL2); $i2++)
295 :     {
296 :     if (($YL1->[$i1]->[2] ne $YL2->[$i2]->[2]) &&
297 :     &close_enough($YL1->[$i1],$YL2->[$i2]))
298 :     {
299 :     push(@pairs,[$YL1->[$i1],$YL2->[$i2]]);
300 :     }
301 :     }
302 :     }
303 :     }
304 :     }
305 :     }
306 :     }
307 :     # @pairs now contains the relevant pairs of PEGs as evidence of coupling
308 :     my $score = &score_pairs($fig,\@pairs);
309 :     return ($score,\@pairs);
310 :     }
311 :     }
312 :    
313 :     sub close_enough {
314 :     my($x,$y) = @_;
315 :    
316 :     return abs($x->[0] - $y->[0]) < 5000;
317 :     }
318 :    
319 :     sub score_pairs {
320 :     my($fig,$pairs) = @_;
321 :    
322 :     my %md5_vals;
323 :     foreach my $pair (@$pairs)
324 :     {
325 :     my $rep = $pair->[0]->[3];
326 :     if (my $rep1 = &representative_of($fig,$pair->[0]->[4],$rep))
327 :     {
328 :     $rep = $rep1;
329 :     }
330 :     $md5_vals{$rep}++;
331 :     }
332 :     return scalar keys(%md5_vals);
333 :     }
334 :    
335 :     sub load_locs {
336 :     my($ss) = @_;
337 :    
338 :     my $file = &summary_file($ss);
339 :     my $locs = {};
340 :     open(TMP,"cut -f1,3 $file |") || die "could not open $file";
341 :     while (defined($_ = <TMP>))
342 :     {
343 :     if ($_ =~ /^(\S+)\t(\S+)/)
344 :     {
345 :     my $md5 = $1;
346 :     my $peg_info = $2;
347 :     my @peg_locs = split(/;/,$peg_info);
348 :     foreach my $peg_loc (@peg_locs)
349 :     {
350 :     if ($peg_loc =~ /^(fig\|(\d+\.\d+)\.peg\.\d+),(\S+)/)
351 :     {
352 :     my $peg = $1;
353 :     my $genome = $2;
354 :     my $loc = $3;
355 :     my($contig,$beg,$end) = &FIG::boundaries_of($loc);
356 :     push(@{$locs->{$genome}->{$contig}},[&FIG::min($beg,$end),&FIG::max($beg,$end),$peg,$md5,$ss]);
357 :     }
358 :     }
359 :     }
360 :     }
361 :     close(TMP);
362 :     return $locs;
363 :     }
364 :    
365 :     sub summary_file {
366 :     my($x) = @_;
367 :    
368 :     my $dir1 = $x % 1000;
369 :     my $summF = "$FIG_Config::data/SimPart/Partitions/$dir1/$x.summary";
370 :     return (-s $summF) ? $summF : undef;
371 :     }
372 :    
373 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3