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

Annotation of /FigKernelPackages/AliTrees.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3