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

Annotation of /FigKernelPackages/AliTrees.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
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 :    
19 :     package AliTrees;
20 :    
21 :     use strict;
22 :     use FIG;
23 :     use Data::Dumper;
24 :     use DBrtns;
25 :     use Carp;
26 :     use AliTree;
27 :     use FIG;
28 :    
29 :     # This is the constructor. Presumably, $class is 'AliTrees'.
30 :     #
31 :    
32 :     sub new {
33 :     my($class,$fig,$data) = @_;
34 :    
35 :     my $ali_trees = {};
36 :    
37 :     $ali_trees->{fig} = defined($fig->{_fig}) ? $fig->{_fig} : $fig;
38 :    
39 :     if (! defined($data))
40 :     {
41 :     if (-d "$FIG_Config::data/AlignmentsAndTrees")
42 :     {
43 :     $data = "$FIG_Config::data/AlignmentsAndTrees";
44 :     }
45 :     else
46 :     {
47 : olson 1.3 warn "Where are the alignments and trees?";
48 :     return undef;
49 : overbeek 1.1 }
50 :     }
51 :    
52 :     $ali_trees->{data} = $data;
53 :    
54 :     bless $ali_trees,$class;
55 :     return $ali_trees;
56 :     }
57 :    
58 :     sub all_alignments {
59 :     my($self) = @_;
60 :    
61 :     my $fig = $self->{fig};
62 :     my $dbh = $fig->db_handle;
63 :    
64 :     my $res = $dbh->SQL("SELECT DISTINCT ali_id from alignments");
65 :     return sort map { $_->[0] } @$res;
66 :     }
67 :    
68 :     sub alignments_containing_peg {
69 :     my($self,$peg) = @_;
70 :    
71 :     my $fig = $self->{fig};
72 :     my $dbh = $fig->db_handle;
73 :    
74 :     my $res = $dbh->SQL("SELECT DISTINCT ali_id from alignments WHERE peg = '$peg'");
75 :     return sort map { $_->[0] } @$res;
76 :     }
77 :    
78 :     sub delete_ali {
79 :     my($self,$ali) = @_;
80 :    
81 :     if (! -d "$self->{data}\/Library/$ali") { return 0 }
82 :     my $data = $self->{data};
83 :     my $fig = $self->{fig};
84 :     my $dir = "$data/Library/$ali";
85 :     &FIG::verify_dir("$data/Deleted");
86 :     my $dbh = $fig->db_handle;
87 :     if (-d "$data/Deleted/$ali") { system "/bin/rm -rf $data/Deleted/$ali" }
88 :     return ($dbh->SQL("DELETE FROM alignments WHERE ali_id = '$ali'") &&
89 :     (system("mv $dir $data/Deleted") == 0));
90 :     }
91 :    
92 :     sub add_ali {
93 :     my($self,$dir_to_add) = @_;
94 :    
95 :     my $data = $self->{data};
96 :     my $fig = $self->{fig};
97 :     my $dbh = $fig->db_handle;
98 :    
99 :     my $ali;
100 :     if (($dir_to_add =~ /([^\/]+)$/) &&
101 :     ($ali = $1) &&
102 :     (! -d "$data/Library/$ali") &&
103 :     (-s "$dir_to_add/full.ali") &&
104 :     (system("cp -r $dir_to_add $data/Library/$ali") == 0))
105 :     {
106 :     my $ali_tree = new AliTree($ali,$fig,$data);
107 :     my $h = $ali_tree->pegs_in_alignment($ali_tree->id);
108 :     foreach my $peg (sort { &FIG::by_fig_id($a,$b) } keys(%$h))
109 :     {
110 :     $dbh->SQL("INSERT INTO alignments (ali_id,peg) VALUES ('$ali','$peg')");
111 :     }
112 :     return 1;
113 :     }
114 :     return 0;
115 :     }
116 :    
117 :     sub merge_ali {
118 :     my($self,$id1,$id2) = @_;
119 :    
120 :     ($id1 =~ /^\d+\.\d+\.peg\.\d+(-(\d+))?$/)
121 :     || return undef;
122 :     my $next_id = &next_id($self,$id1);
123 :    
124 :     my $data = $self->{data};
125 :     my $ali1 = &gjoseqlib::read_fasta("$data/Library/$id1/full.ali");
126 :     my $ali2 = &gjoseqlib::read_fasta("$data/Library/$id2/full.ali");
127 :     my $ali3;
128 :     if ($ali3 = &merge($ali1,$ali2))
129 :     {
130 :     my $tmpdir = "$FIG_Config::var/Temp";
131 :     &FIG::verify_dir($tmpdir);
132 :     mkdir("$tmpdir/$next_id")
133 :     || confess "could not make $tmpdir/$next_id";
134 :     open(ALI,">$tmpdir/$next_id/full.ali")
135 :     || confess "could not open $tmpdir/$next_id/full.ali";
136 :     &gjoseqlib::print_alignment_as_fasta(\*ALI,$ali3);
137 :     close(ALI);
138 :     if ($self->add_ali("$tmpdir/$next_id"))
139 :     {
140 :     system "/bin/rm -r $tmpdir/$next_id";
141 :     return $next_id;
142 :     }
143 :     }
144 :     return undef;
145 :     }
146 :    
147 :     sub merge {
148 :     my($ali1,$ali2) = @_;
149 :    
150 :     my $ln1 = length($ali1->[0]->[2]);
151 :     my $ln2 = length($ali2->[0]->[2]);
152 :    
153 :     my $s;
154 :     my %in1 = map { $s = $_->[2]; $_->[0] => $_->[2] } @$ali1;
155 :     my %in2 = map { $s = $_->[2]; $_->[0] => $_->[2] } @$ali2;
156 :    
157 :     my %common = map { $_->[0] => 1 } grep { $in2{$_->[0]} } @$ali1;
158 :    
159 :     if (keys(%common) == 0)
160 :     {
161 :     return undef;
162 :     }
163 :    
164 :     my($id,$maps_to,%offsets);
165 :     foreach $id (keys(%common))
166 :     {
167 :     my($off1,$off2) = &offsets($in1{$id},$in2{$id});
168 :     if (! defined($off1))
169 :     {
170 :     my $s1 = $in1{$id}; $s1 =~ s/-//g;
171 :     my $s2 = $in2{$id}; $s2 =~ s/-//g;
172 :     print STDERR &Dumper($off1,$off2,$s1,$s2);
173 :     print STDERR "Cannot merge these: probably multidomains that were trimmed: check $id\n";
174 :     return undef;
175 :     }
176 :     $offsets{$id} = [$off1,$off2];
177 :     }
178 :    
179 :     foreach $id (keys(%common))
180 :     {
181 :     my $seq1 = $in1{$id};
182 :     my $seq2 = $in2{$id};
183 :    
184 :     my($off1,$off2) = @{$offsets{$id}};
185 :     my $i1 = &start($off1->[0],$seq1);
186 :     my $i2 = &start($off2->[0],$seq2);
187 :    
188 :     while (($i1 < $ln1) && ($i2 < $ln2))
189 :     {
190 :     while (($i1 < $ln1) && (substr($seq1,$i1,1) eq "-")) { $i1++ }
191 :     while (($i2 < $ln2) && (substr($seq2,$i2,1) eq "-")) { $i2++ }
192 :     if (($i1 < $ln1) && ($i2 < $ln2))
193 :     {
194 :     if (($i1 > $off1->[0]) &&
195 :     ($i2 > $off2->[0]) &&
196 :     (uc substr($seq1,$i1,1) ne uc substr($seq2,$i2,1)))
197 :     {
198 :     print STDERR &Dumper($seq1,$seq2,$i1,$i2,$off1,$off2);
199 :     die "Something is seriously wrong";
200 :     }
201 :     $maps_to->{$i1}->{$i2}++;
202 :     }
203 :     $i1++;
204 :     $i2++;
205 :     }
206 :     }
207 :     # At this point maps_to will map columns in alignment1 to columns in alignment2.
208 :     # Or, more precisely, it gives votes on how to map the columns.
209 :    
210 :     my $whole_map = &build_map($ln1,$ln2,$maps_to);
211 :     # whole_map has mapps from alignment1 to columns in the new alignment, and the same
212 :     # for alignment2
213 :    
214 :     my $ali3 = &build_ali($ali1,$ali2,$whole_map,\%common);
215 :     return $ali3;
216 :     }
217 :    
218 :     sub offsets {
219 :     my($seq1,$seq2) = @_;
220 :    
221 :     $seq1 =~ s/-//g;
222 :     $seq2 =~ s/-//g;
223 :    
224 :     my($beg1,$beg2) = &find_offset($seq1,$seq2);
225 :     if (! defined($beg1)) { return undef }
226 :    
227 :     my $rev_seq1 = reverse $seq1;
228 :     my $rev_seq2 = reverse $seq2;
229 :     my($end1,$end2) = &find_offset($rev_seq1,$rev_seq2);
230 :     if (! defined($end1)) { return undef }
231 :    
232 :     ($end1,$end2) = (length($seq1) - ($end1+1),length($seq2) - ($end2+1));
233 :     return ([$beg1,$end1],[$beg2,$end2]);
234 :     }
235 :    
236 :     sub find_offset {
237 :     my($seq1,$seq2) = @_;
238 :    
239 :     my $off1 = index($seq2,substr($seq1,1,20));
240 :     my $off2 = index($seq1,substr($seq2,1,20));
241 :     if (($off1 <= 0) && ($off2 <= 0))
242 :     {
243 :     return undef;
244 :     }
245 :     elsif (($off1 > 0) && ($off2 <= 0))
246 :     {
247 :     return (0,$off1-1);
248 :     }
249 :     elsif (($off2 > 0) && ($off1 <= 0))
250 :     {
251 :     return ($off2-1,0);
252 :     }
253 :     elsif (($off2 == 1) && ($off1 == 1))
254 :     {
255 :     return (0,0);
256 :     }
257 :     return undef;
258 :     }
259 :    
260 :     sub start {
261 :     my($off,$seq) = @_;
262 :    
263 :     my $i=0;
264 :     while (($i < length($seq)) && (($off > 0) || (substr($seq,$i,1) eq "-")))
265 :     {
266 :     if (substr($seq,$i,1) ne "-") { $off-- }
267 :     $i++;
268 :     }
269 :     return $i;
270 :     }
271 :    
272 :     sub build_map {
273 :     my($ln1,$ln2,$maps_to) = @_;
274 :    
275 :     my $expanded_map = {};
276 :     my @cols_connected = sort { $a <=> $b } keys(%$maps_to);
277 :    
278 :     my($n,$colI,$col,$h,@tuples,$tuple,$to);
279 :     foreach $col (@cols_connected)
280 :     {
281 :     $h = $maps_to->{$col};
282 :     my @poss = sort { ($h->{$b} <=> $h->{$a}) or ($a <=> $b) } keys(%$h);
283 :     push(@tuples,[$col,$poss[0]]);
284 :     }
285 :    
286 :     @tuples = &pins(\@tuples,10);
287 :     foreach $tuple (@tuples)
288 :     {
289 :     ($col,$to) = @$tuple;
290 :     $expanded_map->{$col} = $to;
291 :     }
292 :     my @cols_connected = sort { $a <=> $b } keys(%$expanded_map);
293 :    
294 :     for ($colI=0; ($colI < (@cols_connected - 1)) &&
295 :     ($expanded_map->{$cols_connected[$colI]} < $expanded_map->{$cols_connected[$colI+1]});
296 :     $colI++) {}
297 :    
298 :    
299 :     if ($colI < (@cols_connected - 1))
300 :     {
301 :     print STDERR &Dumper($colI,$cols_connected[$colI],$colI+1,$cols_connected[$colI+1],
302 :     $expanded_map->{$cols_connected[$colI]},$expanded_map->{$cols_connected[$colI+1]}
303 :     );
304 :     confess "mapping is inconsistent";
305 :     }
306 :    
307 :     for ($colI=0; ($colI < (@cols_connected - 1)); $colI++)
308 :     {
309 :     if ((($n = ($cols_connected[$colI+1] - $cols_connected[$colI])) < 5) &&
310 :     (($expanded_map->{$cols_connected[$colI+1]} - $expanded_map->{$cols_connected[$colI]}) == $n))
311 :     {
312 :     $n--;
313 :     while ($n > 0)
314 :     {
315 :     $expanded_map->{$cols_connected[$colI]+$n} = $expanded_map->{$cols_connected[$colI]} + $n;
316 :     $n--;
317 :     }
318 :     }
319 :     }
320 :    
321 :     my($map1,$map2,$i1,$i2,$i3);
322 :     @cols_connected = sort { $a <=> $b } keys(%$expanded_map);
323 :     $map1 = [];
324 :     $map2 = [];
325 :     $i1 = $cols_connected[0];
326 :     $i2 = $expanded_map->{$cols_connected[0]};
327 :     $i3 = 0;
328 :    
329 :     push(@$map1,[$i1,$i3]); ### We begin processing adjacent pairs in mapping
330 :     push(@$map2,[$i2,$i3]); ### with the first coordinate of each pair "already processed"
331 :     $i3++;
332 :    
333 :     my $j;
334 :     for ($colI=0; ($colI < (@cols_connected-1)); $colI++)
335 :     {
336 :     $i1 = $cols_connected[$colI] + 1;
337 :     $n = $cols_connected[$colI+1] - $i1;
338 :     for ($j=0; ($j < $n); $j++)
339 :     {
340 :     push(@$map1,[$i1+$j,$i3]);
341 :     $i3++;
342 :     }
343 :    
344 :     $i2 = $expanded_map->{$cols_connected[$colI]} + 1;
345 :     $n = $expanded_map->{$cols_connected[$colI+1]} - $i2;
346 :     for ($j=0; ($j < $n); $j++)
347 :     {
348 :     push(@$map2,[$i2+$j,$i3]);
349 :     $i3++;
350 :     }
351 :    
352 :     push(@$map1,[$cols_connected[$colI+1],$i3]);
353 :     push(@$map2,[$expanded_map->{$cols_connected[$colI+1]},$i3]);
354 :     $i3++;
355 :     }
356 :     return [$i3,$map1,$map2]; # [length-of-merge,map-for-ali1,map-for-ali2]
357 :     }
358 :    
359 :     sub build_ali {
360 :     my($ali1,$ali2,$map,$common) = @_;
361 :    
362 :     my $ali3 = [];
363 :     my($ln,$map1,$map2) = @$map;
364 :    
365 :     my($x,$entry);
366 :     foreach $x (@$ali1)
367 :     {
368 :     $entry = &expand_entry($ln,$map1,$x);
369 :     push(@$ali3,$entry);
370 :     }
371 :    
372 :     foreach $x (@$ali2)
373 :     {
374 :     if (! $common->{$x->[0]})
375 :     {
376 :     $entry = &expand_entry($ln,$map2,$x);
377 :     push(@$ali3,$entry);
378 :     }
379 :     }
380 :     return $ali3;
381 :     }
382 :    
383 :     sub expand_entry {
384 :     my($ln,$map,$old_entry) = @_;
385 :    
386 :     my $seq = "-" x $ln;
387 :     my $old_seq = $old_entry->[2];
388 :     my $tuple;
389 :     foreach $tuple (@$map)
390 :     {
391 :     my($from,$to) = @$tuple;
392 :     substr($seq,$to,1) = substr($old_seq,$from,1);
393 :     }
394 :     return [$old_entry->[0],"",$seq];
395 :     }
396 :    
397 :     sub set_all {
398 :     my($self) = @_;
399 :    
400 :     my $data = $self->{data};
401 :     my $fig = $self->{fig};
402 :    
403 :     my @all = $self->all_alignments;
404 :     my($id,$ali);
405 :     foreach $id (@all)
406 :     {
407 :     $ali = new AliTree($id,$fig,$data);
408 :     undef $ali;
409 :     }
410 :     }
411 :    
412 :     sub load_db {
413 :     my($self) = @_;
414 :    
415 :     my $data = $self->{data};
416 :     my $fig = $self->{fig};
417 :    
418 :     my $dbf = $fig->db_handle;
419 :     $dbf->drop_table( tbl => 'alignments' );
420 :     $dbf->create_table( tbl => 'alignments',
421 :     flds => 'ali_id varchar(32), peg varchar(32)'
422 :     );
423 :    
424 :     my $tmpdir = "$FIG_Config::var/Temp";
425 :     &FIG::verify_dir($tmpdir);
426 :     open(OUT,">$tmpdir/alignments.$$") || die "could not open $tmpdir/alignments.$$";
427 : overbeek 1.2 opendir(DIR,"$data/Library") || die "could not open $data";
428 : overbeek 1.1 my @alis = grep { $_ !~ /^\./ } readdir(DIR);
429 :     closedir(DIR);
430 :     foreach my $ali (@alis)
431 :     {
432 :     my @pegs = map { ($_ =~ /^(fig\|\d+\.\d+\.peg\.\d+)/) ? $1 : () } `cut -f1 $data/Library/$ali/coords`;
433 :     foreach my $peg (@pegs)
434 :     {
435 :     print OUT "$ali\t$peg\n";
436 :     }
437 :     }
438 :     close(OUT);
439 :     $dbf->load_table( tbl => "alignments", file => "$tmpdir/alignments.$$");
440 :     unlink("$tmpdir/alignments.$$");
441 :    
442 :     $dbf->create_index( idx => "alignments_id_ix",
443 :     tbl => "alignments",
444 :     type => "btree",
445 :     flds => "ali_id" );
446 :     $dbf->create_index( idx => "alignments_peg_ix",
447 :     tbl => "alignments",
448 :     type => "btree",
449 :     flds => "peg" );
450 :     }
451 :    
452 :     sub next_id {
453 :     my($self,$id) = @_;
454 :    
455 :     my $x;
456 :     $id =~ s/\-\d+$//;
457 :     my @id_derivatives = grep { ($_ =~ /^(\d+\.\d+\.peg\.\d+)(-(\d+))?$/) && ($1 eq $id) } $self->all_alignments;
458 :     my $n = 1;
459 :     foreach $x (@id_derivatives)
460 :     {
461 :     if (($x =~ /\-(\d+)$/) && ($n <= $1)) { $n = $1 + 1; }
462 :     }
463 :     return "$id-$n";
464 :     }
465 :    
466 :     sub pins {
467 :     my($tuples,$run_dist) = @_;
468 :    
469 :     my @grouped = &group($tuples,$run_dist);
470 :     my $i = $#grouped;
471 :     while ($i > 0)
472 :     {
473 :     while (($i > 0) &&
474 :     ($i < @grouped) &&
475 :     &cross($grouped[$i-1],$grouped[$i]))
476 :     {
477 :     my $n1 = @{$grouped[$i-1]};
478 :     my $n2 = @{$grouped[$i]};
479 :     if ($n1 < $n2)
480 :     {
481 :     pop(@{$grouped[$i-1]});
482 :     if ($n1 == 1)
483 :     {
484 :     splice(@grouped,$i-1,1);
485 :     $i--;
486 :     }
487 :     }
488 :     else
489 :     {
490 :     pop(@{$grouped[$i]});
491 :     if ($n2 == 1)
492 :     {
493 :     splice(@grouped,$i,1);
494 :     }
495 :     }
496 :     }
497 :     $i--;
498 :     }
499 :     my @left = ();
500 :     my $group;
501 :     foreach $group (@grouped)
502 :     {
503 :     push(@left,@$group);
504 :     }
505 :     return @left;
506 :     }
507 :    
508 :     sub cross {
509 :     my($g1,$g2) = @_;
510 :    
511 :     my($x,$y,$b1,$e1,$b2,$e2);
512 :     if ((@$g1 > 0) && (@$g2 > 0))
513 :     {
514 :     $x = $g1->[@$g1 - 1];
515 :     $y = $g2->[0];
516 :     ($b1,$e1) = @$x;
517 :     ($b2,$e2) = @$y;
518 :     return ((($b1 <= $b2) && ($e1 >= $e2)) ||
519 :     (($b1 >= $b2) && ($e1 <= $e2)));
520 :     }
521 :     return 0;
522 :     }
523 :    
524 :     sub group {
525 :     my($tuples,$run_dist) = @_;
526 :    
527 :     my @groups = ();
528 :     my($i,$group);
529 :     $i = 0;
530 :     while ($i < @$tuples)
531 :     {
532 :     my $group = [$tuples->[$i]];
533 :     $i++;
534 :     while (($i < @$tuples) &&
535 :     ($tuples->[$i]->[0] <= ($group->[$#{$group}]->[0] + $run_dist)) &&
536 :     ($tuples->[$i]->[1] <= ($group->[$#{$group}]->[1] + $run_dist)) &&
537 :     ($tuples->[$i]->[1] > $group->[$#{$group}]->[1]))
538 :     {
539 :     push(@$group,$tuples->[$i]);
540 :     $i++;
541 :     }
542 :     push(@groups,$group);
543 :     }
544 :     return @groups;
545 :     }
546 :    
547 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3