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

Annotation of /FigKernelPackages/AliTrees.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3