[Bio] / FigKernelScripts / svr_update_corr_maps.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/svr_update_corr_maps.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #########################################################################
2 :     use SeedEnv;
3 :     use gjoseqlib;
4 :    
5 :     use strict;
6 :     use Data::Dumper;
7 :     use Carp;
8 :     use CorrTableEntry;
9 :    
10 :     =head1 svr_update_corr_maps
11 :    
12 :     Update correspondence maps after PEG insertions/deletions/function-changes
13 :    
14 :     Note only updates to the given genome are reflected in the updates.
15 :    
16 :     ------
17 :     Example: svr_update_corr_maps -g GenomeDir
18 :    
19 :     would update the 20-column correspondence tables kept in Dir/CorrToReferenceGenomes/*
20 :    
21 :     ------
22 :    
23 :     =head2 Command-Line Options
24 :    
25 :     The program is invoked using
26 :    
27 :     svr_update_corr_maps -g GenomeDir
28 :    
29 :     =over 4
30 :    
31 :     =item -g GenomeDir
32 :    
33 :     This is used to specify a SEED genome directory. If there is a subdirectory
34 :     CorrToReferenceGenomes containing 20-column maps to a set of reference genomes,
35 :     and if PEG insertions or deletions have occurred, this program will update
36 :     the maps (overwriting the proevious set of maps in CorrToReferenceGenomes.
37 :    
38 :     =head2 Output Format
39 :    
40 :     The files in CorrToReferenceGenomes (each containing a map) are updated.
41 :     Otherwise, no changes are made.
42 :    
43 :     =back
44 :    
45 :     =cut
46 :    
47 :     use SeedEnv;
48 :     use SeedUtils;
49 :     use SAPserver;
50 :     use SeedAware;
51 :     use Getopt::Long;
52 :     use ProtSims;
53 :    
54 :     my $usage = "usage: svr_update_corr_maps -g SEEDdirectory\n";
55 :    
56 :     my $genome_dir;
57 :    
58 :     my $rc = GetOptions("g=s" => \$genome_dir );
59 :     if ((! $rc) ||
60 : overbeek 1.2 (! $genome_dir) || (! -s "$genome_dir/Features/peg/tbl")
61 : overbeek 1.1 )
62 :     {
63 :     print STDERR "$usage\n"; exit
64 :     }
65 :     my ($genome1_id) = ($genome_dir =~ /(\d+\.\d+)$/);
66 :    
67 :     if (! opendir(CORR,"$genome_dir/CorrToReferenceGenomes"))
68 :     {
69 :     print STDERR "There is no $genome_dir/CorrToReferenceGenomes directory to update\n";
70 :     exit;
71 :     }
72 :    
73 :     my @genome2s = grep { $_ =~ /^(\d+\.\d+)$/ } readdir(CORR);
74 :     closedir(CORR);
75 :     if (@genome2s == 0)
76 :     {
77 :     print STDERR "There are no $genome_dir/CorrToReferenceGenomes/* to update\n";
78 :     exit;
79 :     }
80 :    
81 :     my $sapObject = SAPserver->new;
82 :    
83 :     my $genome1 = make_genome_source($genome_dir, $sapObject);
84 :     $genome1 or die "Cannot load genome data from $genome1\n";
85 :    
86 :     $genome1->init_data();
87 :     my $functions1 = {};
88 :     $genome1->get_functions($functions1);
89 :     my $aliases1 = {};
90 :     $genome1->get_aliases($aliases1);
91 :    
92 :     foreach my $genome2_name (@genome2s)
93 :     {
94 :     print STDERR "processing $genome2_name\n";
95 :    
96 :     my $genome2 = make_genome_source($genome2_name, $sapObject);
97 :     $genome2 or die "Cannot load genome data from $genome2_name\n";
98 :    
99 :     $genome2->init_data();
100 :     my $functions2 = {};
101 :     $genome2->get_functions($functions2);
102 :     my $aliases2 = {};
103 :     $genome2->get_aliases($aliases2);
104 :    
105 :     my @old_corr = map { chop; [split(/\t/,$_)] } `cat $genome_dir/CorrToReferenceGenomes/$genome2_name`;
106 :     my $corr = &update_correspondence(\@old_corr,$genome1,$functions1,$aliases1,$genome2,$functions2,$aliases2);
107 :     &print_corr($corr,$genome2_name,$genome_dir);
108 :     }
109 :    
110 :     sub update_correspondence {
111 :     my($corr,$genome1,$functions1,$aliases1,$genome2,$functions2,$aliases2) = @_;
112 :    
113 :     my $deleted1 = $genome1->{deleted} ? $genome1->{deleted} : {};
114 :     my $deleted2 = $genome2->{deleted} ? $genome2->{deleted} : {};
115 :     my $sims1 = {};
116 :     my $sims2 = {};
117 :     my $not_clear = {};
118 :     my $in_corr1 = {};
119 :     foreach $_ (@$corr)
120 :     {
121 :     my($id1,$id2,$context_count,$context,$function1,$function2,$aliases1,$aliases2,$bbh,$iden,$psc,
122 :     $b1,$e1,$ln1,$b2,$e2,$ln2,$bit_sc,$mcount) = @$_;
123 :     if ((! $deleted1->{$id1}) && (! $deleted2->{$id2}))
124 :     {
125 :     $in_corr1->{$id1} = 1;
126 :     &update_best($sims1,$id1,$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$ln1,$ln2);
127 :     &update_best($sims2,$id2,$id1,$iden,$psc,$bit_sc,$b2,$e2,$b1,$e1,$ln2,$ln1);
128 :     }
129 :     }
130 :     my @possibly_inserted_pegs = grep { ! $in_corr1->{$_} } keys(%$aliases1);
131 :     if (@possibly_inserted_pegs > 0)
132 :     {
133 :     my @sims_new = &get_new_sims($genome1,$genome2,\@possibly_inserted_pegs);
134 :     foreach my $sim (@sims_new)
135 :     {
136 :     &update_best($sims1,$sim->id1,$sim->id2,$sim->iden,$sim->psc,$sim->bitsc,$sim->b1,$sim->e1,$sim->b2,$sim->e2,$sim->ln1,$sim->ln2);
137 :     &update_best($sims2,$sim->id2,$sim->id1,$sim->iden,$sim->psc,$sim->bitsc,$sim->b2,$sim->e2,$sim->b1,$sim->e1,$sim->ln2,$sim->ln1);
138 :     }
139 :     }
140 :    
141 :     &set_best_and_not_clear($sims1,$not_clear);
142 :     &set_best_and_not_clear($sims2,$not_clear);
143 :     $corr = &build_corr($genome1,$sims1,$functions1,$aliases1,$genome2,$sims2,$functions2,$aliases2,5,1000000,$not_clear);
144 :     return $corr;
145 :     }
146 :    
147 :     sub get_new_sims {
148 :     my($genome1,$genome2,$inserted_pegs) = @_;
149 :    
150 :     my %ins_pegs = map { $_ => 1 } @$inserted_pegs;
151 :     my $tmp_dir = SeedAware::location_of_tmp();
152 :     my $tmp1 = "$tmp_dir/tmp1_$$.fasta";
153 :     my $tmp2 = "$tmp_dir/tmp2_$$.fasta";
154 :     my $lens1 = $genome1->get_fasta($tmp1);
155 :     my $lens2 = $genome2->get_fasta($tmp2);
156 :     my @seqs1 = grep { $ins_pegs{$_->[0]} } &gjoseqlib::read_fasta($tmp1);
157 :     my @sims = &ProtSims::blastP(\@seqs1,$tmp2,1,1); # this last argument forces the use of blast, bypassing blat
158 :     return @sims;
159 :     }
160 :    
161 :     sub print_corr {
162 :     my($corr,$genome2_name,$genome_dir) = @_;
163 :    
164 :     my $file = "$genome_dir/CorrToReferenceGenomes/$genome2_name";
165 :     rename($file,$file . "~") || die "could not rename $file";
166 :     open(CORR,">$file") || die "could not open $file";
167 :     foreach $_ (sort { &SeedUtils::by_fig_id($a->[0],$b->[0]) } @$corr)
168 :     {
169 :     print CORR join("\t",@$_),"\n";
170 :     }
171 :     close(CORR);
172 :     }
173 :    
174 :     sub build_corr {
175 :     my($genome1,$sims1,$functions1,$aliases1,$genome2,$sims2,$functions2,$aliases2,$sz_context,$ignore_ov,$not_clear) = @_;
176 :    
177 :     my $corr = [];
178 :    
179 :     my($matching_context, $matching_count) =
180 :     &matching_neighbors($genome1,$sims1,$functions1,$genome2,$sims2,$functions2,$sz_context,$ignore_ov);
181 :    
182 :     foreach my $peg1 (keys(%$sims1))
183 :     {
184 :     my $peg2 = $sims1->{$peg1}->[0];
185 :     if ($peg2)
186 :     {
187 :     my $context = "";
188 :     my $context_count = 0;
189 :     my $function2 = "";
190 :     my $aliases_peg2 = "";
191 :    
192 :     my $function1 = $functions1->{$peg1} ? $functions1->{$peg1} : "";
193 :     my $aliases1 = $aliases1->{$peg1} ? $aliases1->{$peg1} : "";
194 :     my $peg3 = $sims2->{$peg2}->[0];
195 :     my $bbh = ($peg3 && ($peg3 eq $peg1)) ? "<=>" : "->";
196 :    
197 :     if ($_ = $matching_context->{"$peg1,$peg2"})
198 :     {
199 :     $context = $_;
200 :     $context_count = ($context =~ tr/,//) + 1;
201 :     }
202 :    
203 :     my($iden,$psc,$b1,$e1,$b2,$e2,$ln1,$ln2,$bitsc);
204 :     (undef,$iden,$psc,$bitsc,$b1,$e1,$b2,$e2,$ln1,$ln2) = @{$sims1->{$peg1}};
205 :     if ($functions2->{$peg2}) { $function2 = $functions2->{$peg2} }
206 :     if ($aliases2->{$peg2}) { $aliases_peg2 = $aliases2->{$peg2} }
207 :     my $mcount = $matching_count->{"$peg1,$peg2"};
208 :     $mcount = 0 unless defined($mcount);
209 :     push(@$corr,[$peg1,$peg2,$context_count,$context,$function1,$function2,
210 :     $aliases1,$aliases_peg2,$bbh,$iden,$psc,
211 :     $b1,$e1,$ln1,$b2,$e2,$ln2,$bitsc,$mcount,
212 :     (($bbh ne "<=>") || $not_clear->{$peg1}) ? 0 : 1]);
213 :     }
214 :     }
215 :    
216 :     foreach my $peg2 (keys(%$sims2))
217 :     {
218 :     my $peg1 = $sims2->{$peg2}->[0];
219 :     if ($peg1)
220 :     {
221 :     my $context = "";
222 :     my $context_count = 0;
223 :     my $function1 = "";
224 :     my $aliases1 = "";
225 :    
226 :     my $function2 = $functions2->{$peg2} ? $functions2->{$peg2} : "";
227 :     my $aliases_peg2 = $aliases2->{$peg2} ? $aliases2->{$peg2} : "";
228 :     my $peg3 = $sims1->{$peg1}->[0];
229 :     if ($peg3 ne $peg2)
230 :     {
231 :     if ($_ = $matching_context->{"$peg1,$peg2"})
232 :     {
233 :     $context = $_;
234 :     $context_count = ($context =~ tr/,//) + 1;
235 :     }
236 :     my $mcount = $matching_count->{"$peg1,$peg2"};
237 :     $mcount = 0 unless defined($mcount);
238 :     my($iden,$psc,$b1,$e1,$b2,$e2,$ln1,$ln2,$bitsc);
239 :     (undef,$iden,$psc,$bitsc,$b2,$e2,$b1,$e1,$ln2,$ln1) = @{$sims2->{$peg2}};
240 :     if ($functions1->{$peg1}) { $function1 = $functions1->{$peg1} }
241 :     if ($aliases1->{$peg1}) { $aliases1 = $aliases1->{$peg1} }
242 :    
243 :     push(@$corr,[$peg1,$peg2,$context_count,$context,$function1,$function2,
244 :     $aliases1,$aliases_peg2,"<-",$iden,$psc,
245 :     $b1,$e1,$ln1,$b2,$e2,$ln2,$bitsc,$mcount,0]);
246 :     }
247 :     }
248 :     }
249 :     return $corr;
250 :     }
251 :    
252 :     sub merge {
253 :     my($b1a,$e1a,$b2a,$e2a,$b1b,$e1b,$b2b,$e2b,$bsc1,$bsc2) = @_;
254 :    
255 :     if (($b1a < $b1b) && (abs($b1b - $e1a) < 10) &&
256 :     ($b2a < $b2b) && (abs($b2b - $e2a) < 10))
257 :     {
258 :     return ($b1a,$e1b,$b2a,$e2b);
259 :     }
260 :     elsif (($b1b < $b1a) && (abs($b1a - $e1b) < 10) &&
261 :     ($b2b < $b2a) && (abs($b2a - $e2b) < 10))
262 :     {
263 :     return ($b1b,$e1a,$b2b,$e2a);
264 :     }
265 :     elsif ($bsc1 >= $bsc2)
266 :     {
267 :     return ($b1a,$e1a,$b2a,$e2a);
268 :     }
269 :     else
270 :     {
271 :     return ($b1b,$e1b,$b2b,$e2b);
272 :     }
273 :     }
274 :    
275 :     sub matching_neighbors {
276 :     my($genome1,$sims1,$functions1,$genome2,$sims2,$functions2,$sz_context,$ignore_ov) = @_;
277 :    
278 :     my %by_genome;
279 :     my @peg_loc_tuples_in_genome;
280 :    
281 :     @peg_loc_tuples_in_genome = $genome1->get_peg_loc_tuples();
282 :     $by_genome{$genome1} = &set_neighbors(\@peg_loc_tuples_in_genome, $sz_context, $ignore_ov);
283 :    
284 :     @peg_loc_tuples_in_genome = $genome2->get_peg_loc_tuples();
285 :     $by_genome{$genome2} = &set_neighbors(\@peg_loc_tuples_in_genome, $sz_context, $ignore_ov);
286 :    
287 :     my %matched_pairs;
288 :     my %matching_count;
289 :     foreach my $peg1 (keys(%$sims1))
290 :     {
291 :     my $peg2 = $sims1->{$peg1}->[0];
292 :     my $neigh1 = $by_genome{$genome1}->{$peg1};
293 :     my $neigh2 = $by_genome{$genome2}->{$peg2};
294 :     my %neigh2H = map { $_ => 1 } @$neigh2;
295 :     my @pairs = ();
296 :     my $matching_count = 0;
297 :     foreach my $n1 (@$neigh1)
298 :     {
299 :     my $maps_to = $sims1->{$n1}->[0];
300 :     if ($maps_to && $neigh2H{$maps_to})
301 :     {
302 :     push(@pairs,"$n1:$maps_to");
303 :    
304 :     if ($functions1->{$n1} eq $functions2->{$maps_to})
305 :     {
306 :     $matching_count++;
307 :     }
308 :     }
309 :     }
310 :     $matched_pairs{"$peg1,$peg2"} = join(",",@pairs);
311 :     $matching_count{"$peg1,$peg2"} = $matching_count;
312 :     }
313 :     return \%matched_pairs, \%matching_count;
314 :     }
315 :    
316 :     sub set_best_and_not_clear {
317 :     my($sims,$not_clear) = @_;
318 :    
319 :     foreach my $id (keys(%$sims))
320 :     {
321 :     my $bestL = $sims->{$id};
322 :     my $best = $bestL->[0];
323 :     if ($best->[2] > 1.0e-20)
324 :     {
325 :     $not_clear->{$id} = 1 ;
326 :     }
327 :    
328 :     if (! &ok_len($best))
329 :     {
330 :     $not_clear->{$id} = 1; # poor coverage -> not-clear
331 :     }
332 :    
333 :     if (@$bestL > 1)
334 :     {
335 :     if (&ok_len($bestL->[1]) &&
336 :     (abs($best->[1] - $bestL->[1]->[1]) < 5)) # diff in identity is < 5 -> not-clear
337 :     {
338 :     $not_clear->{$id} = 1;
339 :     }
340 :     }
341 :     $sims->{$id} = $best;
342 :     }
343 :     }
344 :    
345 :     sub ok_len {
346 :     my($sim) = @_;
347 :    
348 :     my($b1,$e1,$b2,$e2,$ln1,$ln2) = @{$sim}[4..9];
349 :     return ((($e1-$b1) >= (0.8 * $ln1)) && (($e2-$b2) >= (0.8 * $ln2)));
350 :     }
351 :    
352 :     sub better {
353 :     my($psc1,$psc2,$bsc1,$bsc2) = @_;
354 :    
355 :     return (($psc1 < $psc2) || (($psc1 == $psc2) && ($bsc1 > $bsc2)));
356 :     }
357 :    
358 :     sub update_best {
359 :     my($sims1,$id1,$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$len1,$len2) = @_;
360 :    
361 :     my $x = $sims1->{$id1};
362 :     if (! $x)
363 :     {
364 :     $sims1->{$id1} = [[$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$len1,$len2]];
365 :     }
366 :     elsif ($x && (($x->[0]->[0] ne $id2) && &better($psc,$x->[0]->[2],$bit_sc,$x->[0]->[3])))
367 :     {
368 :     splice(@$x,0,0,[$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$len1,$len2]);
369 :     $#{$x} = 1;
370 :     }
371 :     elsif ($x && ($x->[0]->[0] eq $id2))
372 :     {
373 :     ($b1,$e1,$b2,$e2) = &merge($b1,$e1,$b2,$e2,$x->[0]->[4],$x->[0]->[5],$x->[0]->[6],$x->[0]->[7],$bit_sc,$x->[0]->[3]);
374 :     $x->[0]->[4] = $b1;
375 :     $x->[0]->[5] = $e1;
376 :     $x->[0]->[6] = $b2;
377 :     $x->[0]->[7] = $e2;
378 :     }
379 :     elsif ((@$x == 1) && (($x->[0]->[0] ne $id2) && (! &better($psc,$x->[0]->[2],$bit_sc,$x->[0]->[3]))))
380 :     {
381 :     $x->[1] = [$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$len1,$len2];
382 :     }
383 :     elsif ((@$x == 2) && (($x->[1]->[0] ne $id2) && ($psc < $x->[1]->[2])))
384 :     {
385 :     $x->[1] = [$id2,$iden,$psc,$bit_sc,$b1,$e1,$b2,$e2,$len1,$len2];
386 :     }
387 :     }
388 :    
389 :     sub compare_locs {
390 :     my($loc1,$loc2) = @_;
391 :    
392 :     my($contig1,$min1,$max1) = &SeedUtils::boundaries_of($loc1);
393 :     my($contig2,$min2,$max2) = &SeedUtils::boundaries_of($loc2);
394 :     return (($contig1 cmp $contig2) or (($min1+$max1) <=> ($min2+$max2)));
395 :     }
396 :    
397 :     sub set_neighbors {
398 :     my($peg_loc_tuples,$N,$ignore_ov) = @_;
399 :    
400 :     my $peg_to_neighbors = {};
401 :     my $i;
402 :    
403 :     for ($i=0; ($i < @$peg_loc_tuples); $i++)
404 :     {
405 :     my($contigI,$minI,$maxI) = &SeedUtils::boundaries_of($peg_loc_tuples->[$i]->[1]);
406 :     $contigI || confess "BAD";
407 :     my $neighbors = [];
408 :     my $j = $i-1;
409 :     my $to_add = $N;
410 :     while (($j >= 0) && ($to_add > 0) &&
411 :     &same_contig($peg_loc_tuples->[$j]->[1],$contigI))
412 :     {
413 :     $j--;
414 :     if (&distinct($peg_loc_tuples->[$j]->[1],$peg_loc_tuples->[$j+1]->[1],$ignore_ov))
415 :     {
416 :     $to_add--;
417 :     }
418 :     }
419 :     $j++;
420 :     while ($j < $i) { push(@$neighbors,$peg_loc_tuples->[$j]->[0]); $j++ }
421 :    
422 :     $j = $i+1;
423 :     $to_add = $N;
424 :     while (($j < @$peg_loc_tuples) && ($to_add > 0) &&
425 :     &same_contig($peg_loc_tuples->[$j]->[1],$contigI))
426 :     {
427 :     push(@$neighbors,$peg_loc_tuples->[$j]->[0]);
428 :     if (&distinct($peg_loc_tuples->[$j]->[1],$peg_loc_tuples->[$j-1]->[1],$ignore_ov))
429 :     {
430 :     $to_add--;
431 :     }
432 :     $j++;
433 :     }
434 :     $peg_to_neighbors->{$peg_loc_tuples->[$i]->[0]} = $neighbors;
435 :     }
436 :     return $peg_to_neighbors;
437 :     }
438 :    
439 :     sub distinct {
440 :     my($x,$y,$ignore_ov) = @_;
441 :    
442 :     return ($ignore_ov > &overlap($x,$y));
443 :     }
444 :    
445 :     sub overlap {
446 :     my($x,$y) = @_;
447 :    
448 :     my($contig1,$min1,$max1) = &SeedUtils::boundaries_of($x);
449 :     my($contig2,$min2,$max2) = &SeedUtils::boundaries_of($y);
450 :     if ($contig1 ne $contig2) { return 0 }
451 :     if (&SeedUtils::between($min1,$min2,$max1)) { return ($max1 - $min2)+1 }
452 :     if (&SeedUtils::between($min2,$min1,$max2)) { return ($max2 - $min1)+1 }
453 :     return 0;
454 :     }
455 :    
456 :     sub same_contig {
457 :     my($entry,$contig) = @_;
458 :    
459 :     $contig || confess "BAD";
460 :     my($contig1,$minI,$maxI) = &SeedUtils::boundaries_of($entry);
461 :     return ($contig eq $contig1);
462 :     }
463 :    
464 :     sub make_genome_source
465 :     {
466 :     my($name, $sap) = @_;
467 :    
468 :     if ($name =~ /^\d+\.\d+$/)
469 :     {
470 :     return SapGenomeSource->new($name, $sap);
471 :     }
472 :     elsif (-d $name)
473 :     {
474 :     my %deleted;
475 :     if (-s "$name/Features/peg.deleted.features")
476 :     {
477 :     %deleted = map { $_ =~ /^(\S+)/; ($1 => 1) } `cut -f1 $name/Features/peg/deleted.features`;
478 :     }
479 :    
480 :     my $fasta = "$name/Features/peg/fasta";
481 :    
482 :     my $tbl = "$name/Features/peg/tbl";
483 :    
484 :     if (! -f $fasta)
485 :     {
486 :     die "No fasta found in $fasta\n";
487 :     }
488 :     if (! -f $tbl)
489 :     {
490 :     die "No tbl file found in $tbl\n";
491 :     }
492 :    
493 :     #
494 :     # We might have multiple functions files. If there are,
495 :     # collapse into one based on the usual RAST
496 :     # rules (old assigned funcs overwritten by
497 :     # auto-assign funcs overwritten by FIGfams funcs).
498 :     #
499 :    
500 :     my @files;
501 :     foreach my $file (qw(assigned_functions proposed_non_ff_functions proposed_functions))
502 :     {
503 :     if (-f "$name/$file")
504 :     {
505 :     push(@files, "$name/$file")
506 :     }
507 :     }
508 :     if (@files == 0)
509 :     {
510 :     die "No functions file found for $name\n";
511 :     }
512 :     return SapFileSource->new($fasta, $tbl, \@files, \%deleted);
513 :     }
514 :     else
515 :     {
516 :     #
517 :     # Must be a comma-sep triple
518 :     #
519 :     my @x = split(/,/, $name);
520 :     if (@x != 3)
521 :     {
522 :     die "Invalid genome specifier: $name\n";
523 :     }
524 :     my($fasta, $tbl, $func) = @x;
525 :     if (! -f $fasta)
526 :     {
527 :     die "No fasta found in $fasta\n";
528 :     }
529 :     if (! -f $tbl)
530 :     {
531 :     die "No tbl file found in $tbl\n";
532 :     }
533 :     if (! -f $func)
534 :     {
535 :     die "No function file found in $func\n";
536 :     }
537 :    
538 :     return SapFileSource->new($fasta, $tbl, [$func], {});
539 :     }
540 :     }
541 :    
542 :     package SapFileSource;
543 :     use strict;
544 :    
545 :     sub new
546 :     {
547 :     my($class, $fasta, $tbl, $func_files,$deleted) = @_;
548 :     my $self = {
549 :     fasta => $fasta,
550 :     tbl => $tbl,
551 :     deleted => $deleted,
552 :     func_files => $func_files,
553 :     };
554 :     bless $self, $class;
555 :     return $self;
556 :     }
557 :    
558 :     sub init_data
559 :     {
560 :     my($self) = @_;
561 :    
562 :     my $deleted = $self->{deleted} ? $self->{deleted} : {};
563 :     open(TBL, "<", $self->{tbl}) or die "Cannot read $self->{tbl}: $!";
564 :    
565 :     while (<TBL>)
566 :     {
567 :     chomp;
568 :     my($id, $loc, @aliases) = split(/\t/);
569 :     if (! $deleted->{$id})
570 :     {
571 :     $self->{aliases}->{$id} = [@aliases];
572 :     }
573 :     }
574 :     close(TBL);
575 :     }
576 :    
577 :     sub get_fasta
578 :     {
579 :     my($self, $file) = @_;
580 :    
581 :     my $deleted = $self->{deleted} ? $self->{deleted} : {};
582 :     my @seqs = grep { ! $deleted->{&SeedUtils::genome_of($_->[0])} } &gjoseqlib::read_fasta($self->{fasta});
583 :     &gjoseqlib::print_alignment_as_fasta($file,\@seqs);
584 :     my $lens = {};
585 :     foreach $_ (@seqs) { $lens->{$_->[0]} = length($_->[2]) }
586 :     return $lens;
587 :     }
588 :    
589 :     sub get_functions
590 :     {
591 :     my($self, $hash) = @_;
592 :     my $deleted = $self->{deleted} ? $self->{deleted} : {};
593 :    
594 :     for my $file (@{$self->{func_files}})
595 :     {
596 :     open(FFILE, "<", $file) or die "Cannot read $file: $!";
597 :     while (<FFILE>)
598 :     {
599 :     chomp;
600 :     my($id, $fn) = split(/\t/);
601 :     if (! $deleted->{$id})
602 :     {
603 :     $hash->{$id} = $fn;
604 :     }
605 :     }
606 :     close(FFILE);
607 :     }
608 :     }
609 :    
610 :     sub get_aliases
611 :     {
612 :     my($self, $hash) = @_;
613 :     my $deleted = $self->{deleted} ? $self->{deleted} : {};
614 :     for my $key (keys % {$self->{aliases}})
615 :     {
616 :     if (! $deleted->{$key})
617 :     {
618 :     $hash->{$key} = join(",", @{$self->{aliases}->{$key}});
619 :     }
620 :     }
621 :     }
622 :    
623 :     sub get_peg_loc_tuples
624 :     {
625 :     my($self) = @_;
626 :     my $deleted = $self->{deleted} ? $self->{deleted} : {};
627 :     my @all_fids = grep { ! $deleted->{$_} } keys(%{$self->{loc}});;
628 :    
629 :     my @peg_loc_tuples_in_genome =
630 :     sort { &main::compare_locs($a->[1],$b->[1]) }
631 :     map { [$_, [split(/,/,$self->{loc}->{$_})]] }
632 :     @all_fids;
633 :     return @peg_loc_tuples_in_genome;
634 :     }
635 :    
636 :     package SapGenomeSource;
637 :    
638 :     sub new
639 :     {
640 :     my($class, $genome, $sap) = @_;
641 :     my $self = {
642 :     genome => $genome,
643 :     sap => $sap,
644 :     };
645 :    
646 :     return bless $self, $class;
647 :     }
648 :    
649 :     sub init_data
650 :     {
651 :     my($self) = @_;
652 :     my $genome = $self->{genome};
653 :     my $sap = $self->{sap};
654 :    
655 :     my $fidHash = $sap->all_features(-ids => $genome, -type => 'peg');
656 :     $self->{all_fids} = $fidHash->{$genome};
657 :    
658 :     if (@{$self->{all_fids}} == 0)
659 :     {
660 :     die "Could not load pegs for $genome\n";
661 :     }
662 :    
663 :     my $locHash = $sap->fid_locations(-ids => $self->{all_fids});
664 :     }
665 :    
666 :     sub get_functions
667 :     {
668 :     my($self, $hash) = @_;
669 :    
670 :     my $fns = $self->{sap}->ids_to_functions(-ids => $self->{all_fids});
671 :     $hash->{$_} = $fns->{$_} for keys %$fns;
672 :     }
673 :    
674 :     sub get_aliases {
675 :     my($self, $aliases) = @_;
676 :    
677 :     my $aliasHash = $sapObject->fids_to_ids(-ids => $self->{all_fids});
678 :    
679 :     foreach my $peg (keys %$aliasHash)
680 :     {
681 :     my $typeH = $aliasHash->{$peg} ? $aliasHash->{$peg} : {};
682 :     my @all_aliases = map { @{$typeH->{$_}} } keys(%$typeH);
683 :     my $aliasStr = (@all_aliases > 0) ? join(",",@all_aliases) : "";
684 :     $aliases->{$peg} = $aliasStr;
685 :     }
686 :     return $aliases;
687 :     }
688 :    
689 :     sub get_peg_loc_tuples
690 :     {
691 :     my($self) = @_;
692 :    
693 :     my $locHash = $sapObject->fid_locations(-ids => $self->{all_fids});
694 :     my @peg_loc_tuples_in_genome =
695 :     sort { &main::compare_locs($a->[1],$b->[1]) }
696 :     map { [$_,$locHash->{$_}] }
697 :     keys(%$locHash);
698 :     return @peg_loc_tuples_in_genome;
699 :     }
700 :    
701 :     sub get_fasta {
702 :     my($self, $file) = @_;
703 :    
704 :     my $lens = {};
705 :     my $fastaHash = $sapObject->ids_to_sequences(-ids => $self->{all_fids},
706 :     -protein => 1);
707 :    
708 :     open(FASTA,">$file") || die "could not open $file";
709 :     foreach my $peg (keys(%$fastaHash))
710 :     {
711 :     my $seq = $fastaHash->{$peg};
712 :     if ($seq)
713 :     {
714 :     print FASTA ">$peg\n$seq\n";
715 :     $lens->{$peg} = length($seq);
716 :     }
717 :     }
718 :     close(FASTA);
719 :     return $lens;
720 :     }
721 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3