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

Annotation of /FigKernelPackages/AliTree.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 AliTree;
20 :    
21 :     use strict;
22 :     use FIG;
23 :     use SameFunc;
24 :     use gjoseqlib;
25 :     use gjoalignment;
26 :     use tree_utilities;
27 :     use PHOB;
28 :     use Tracer;
29 :     use AliTrees;
30 :    
31 :     use Carp;
32 :     use Data::Dumper;
33 :    
34 :    
35 :     #### To turn on tracing, at the command line type
36 :     ###
37 :     ### export TRACING=Ross
38 :     ### trace 3 AliTree
39 :     ###
40 :    
41 :     sub new {
42 :     my($class,$id,$fig,$data) = @_;
43 :    
44 :     ($id =~ /^\d+\.\d+\.peg\.\d+/) || confess "invalid id: $id";
45 :     if (! defined($fig)) { $fig = new FIG }
46 :     if (! defined($data))
47 :     {
48 :     if (-d "$FIG_Config::data/AlignmentsAndTrees")
49 :     {
50 :     $data = "$FIG_Config::data/AlignmentsAndTrees";
51 :     }
52 :     else
53 :     {
54 :     confess "Where are the alignments and trees?";
55 :     }
56 :     }
57 :    
58 :     my $ali_tree = {};
59 :     $ali_tree->{id} = $id;
60 :     $ali_tree->{fig} = $fig;
61 :     $ali_tree->{data} = $data;
62 :     my $dir = $ali_tree->{dir} = "$data/Library/$id";
63 :    
64 :     if (! defined("$dir/full.ali")) { return undef }
65 :    
66 :     my($ali,$core_ali);
67 : overbeek 1.3 open(ALI,"<$dir/full.ali") || die "could not open $dir/full.ali";
68 :     ($ali = &gjoseqlib::read_fasta(\*ALI)) || die "$dir/full.ali is not well-formatted fasta";
69 : overbeek 1.1 close(ALI);
70 :    
71 :     if (open(ALI,"<$dir/core.ali") && ($core_ali = &gjoseqlib::read_fasta(\*ALI)))
72 :     {
73 :     close(ALI);
74 :     }
75 :     else
76 :     {
77 :     $core_ali = undef;
78 :     }
79 :     my @deleted = grep {$fig->is_deleted_fid($_->[0]) } @$ali;
80 :    
81 :     my %deleted;
82 :     if (@deleted > 0)
83 :     {
84 :     %deleted = map { $_->[0] => 1 } @deleted;
85 :     $ali = [grep { ! $deleted{$_->[0]} } @$ali];
86 :     $core_ali = $core_ali ? [grep { ! $deleted{$_->[0]} } @$core_ali] : undef;
87 :     }
88 : overbeek 1.3
89 : overbeek 1.1 if (@$ali < 2) { return undef }
90 :    
91 :     my @coords = ();
92 :     if (open(COORDS,"<$dir/coords"))
93 :     {
94 :     @coords = map { chomp; [split(/\t/,$_)] } <COORDS>;
95 :     close(COORDS);
96 :     }
97 :    
98 :     if (@coords > @$ali)
99 :     {
100 :     @coords = grep { ! $deleted{$_->[0]} } @coords;
101 :     open(COORDS,">$dir/coords") || die "could not update $dir/coords";
102 :     foreach my $x (@coords)
103 :     {
104 :     print COORDS join("\t",@$x),"\n";
105 :     }
106 :     close(COORDS);
107 :     }
108 :     elsif (@coords < @$ali)
109 :     {
110 :     open(COORDS,">$dir/coords") || die "could not update $dir/coords";
111 :     for (my $tupleI = (@$ali - 1); ($tupleI >= 0); $tupleI--)
112 :     {
113 :     my $tuple = $ali->[$tupleI];
114 :     my($peg,undef,$seq) = @$tuple;
115 :     $seq =~ s/-//g;
116 :     $seq = lc $seq;
117 :     my $full_seq = lc $fig->get_translation($peg);
118 :     my($offset,$n);
119 :     if (($offset = index($full_seq,$seq)) >= 0)
120 :     {
121 :     print COORDS join("\t",($peg,$offset+1,$offset + length($seq),length($full_seq))),"\n";
122 :     push(@coords,[$peg,$offset+1,$offset + length($seq),length($full_seq)]);
123 :     }
124 :     elsif (($n = int(0.4 * length($full_seq))) &&
125 :     (($offset = index($full_seq,substr($seq,$n))) >= 0))
126 :     {
127 :     $offset -= $n;
128 :     print COORDS join("\t",($peg,$offset+1,$offset + length($seq),length($full_seq))),"\n";
129 :     push(@coords,[$peg,$offset+1,$offset + length($seq),length($full_seq)]);
130 :     }
131 :     else
132 :     {
133 :     push(@deleted,$peg);
134 :     print STDERR "SERIOUS ERROR: ",&Dumper($peg,$seq,$full_seq);
135 :     splice(@$ali,$tupleI,1);
136 :     }
137 :     }
138 :     close(COORDS);
139 :     }
140 :    
141 :     if (@deleted > 0)
142 :     {
143 :     open(ALI,">$dir/full.ali") || die "could not update $dir/full.ali";
144 :     &gjoseqlib::print_alignment_as_fasta(\*ALI,$ali);
145 :     close(ALI);
146 :    
147 :     if (open(ALI,">$dir/core.ali"))
148 :     {
149 :     &gjoseqlib::print_alignment_as_fasta(\*ALI,$core_ali);
150 :     close(ALI);
151 :     }
152 :     }
153 :     $ali_tree->{ali} = $ali;
154 :     ############# we now have an up-to-date alignment ##############
155 :    
156 :     my %coords = map { $_->[0] => [$_->[1],$_->[2],$_->[3]] } @coords;
157 :     $ali_tree->{coords} = \%coords;
158 :     ############# we now have up-to-date coords ##############
159 :    
160 :     my $tree = "";
161 :     # if ((! -s "$dir/tree.newick") && (@$ali > 3))
162 :     # {
163 :     # &FIG::run("make_neigh_joining_tree -a 2 $dir/full.ali > $dir/tree.newick");
164 :     # }
165 :    
166 :     if ((@$ali > 3) && (-s "$dir/tree.newick"))
167 :     {
168 :     $tree = &tree_utilities::parse_newick_tree(join("",`cat $dir/tree.newick`));
169 :     if (@deleted > 0)
170 :     {
171 :     my $ids = &tree_utilities::tips_of_tree($tree);
172 :    
173 :     if (@$ids > @$ali)
174 :     {
175 :     my %good = map {$_->[0] => 1} @$ali;
176 :     $tree = &tree_utilities::subtree($tree,\%good);
177 :     open(TREE,">$dir/tree.newick")
178 :     || die "could not update $dir/tree.newick";
179 :     print TREE &tree_utilities::to_newick($tree),"\n";
180 :     close(TREE);
181 :     }
182 :     }
183 :     }
184 :     $ali_tree->{tree} = $tree;
185 :     ############# we now have an up-to-date tree ##############
186 :    
187 :     bless $ali_tree,$class;
188 :     return $ali_tree;
189 :     }
190 :    
191 :     sub tree {
192 :     my($self) = @_;
193 :    
194 :     if (! $self->{tree})
195 :     {
196 :     my $id = $self->{id};
197 :     my $data = $self->{data};
198 : overbeek 1.2 my $dir = $self->{dir};
199 : overbeek 1.1 my $ali = $self->{ali};
200 :     if ((! -s "$dir/tree.newick") && (@$ali > 3))
201 :     {
202 :     &FIG::run("make_neigh_joining_tree -a 2 $dir/full.ali > $dir/tree.newick");
203 :     $self->{tree} = &tree_utilities::parse_newick_tree(join("",`cat $dir/tree.newick`));
204 :     }
205 :     }
206 :     return $self->{tree};
207 :     }
208 :    
209 :     sub id {
210 :     my($self) = @_;
211 :    
212 :     return $self->{id};
213 :     }
214 :    
215 :     sub ali {
216 :     my($self) = @_;
217 :    
218 :     return $self->{ali};
219 :     }
220 :    
221 :     sub pegs_in_alignment {
222 :     my($self) = @_;
223 :    
224 :     return $self->{coords};
225 :     }
226 :    
227 : overbeek 1.4 sub functions_in {
228 :     my($self) = @_;
229 :     my $fig = $self->{fig};
230 :    
231 :     my($pegH,$peg,%funcs);
232 :     $pegH = $self->{coords};
233 :     foreach $peg (keys(%$pegH))
234 :     {
235 :     my $func = $fig->function_of($peg);
236 :     $funcs{$func}++;
237 :     }
238 :     return map { [$funcs{$_},$_] } sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
239 :     }
240 :    
241 : overbeek 1.1 sub display {
242 :     my($self) = @_;
243 :    
244 :     my @lines = ();
245 :     my $ali = $self->{ali};
246 :     foreach my $tuple (@$ali)
247 :     {
248 :     push(@lines,sprintf("%32s",$tuple->[0]) . "\t" . $tuple->[2] . "\n");
249 :     }
250 :     return join("",@lines);
251 :     }
252 :    
253 :     sub indel_ratio {
254 :     my($self) = @_;
255 :    
256 :     my $seq = join("",map { $_->[2] } @{$self->{ali}});
257 :     my $tot = length($seq);
258 :     my $indels = ($seq =~ tr/-//);
259 :     my $inrat = sprintf "%.3f", $indels/$tot;
260 :     return $inrat;
261 :     }
262 :    
263 :     sub overlaps {
264 :     my($self) = @_;
265 :    
266 :     my $id = $self->{id};
267 :     my $fig = $self->{fig};
268 :     my $data = $self->{data};
269 :     my $coords = $self->{coords};
270 :     my @pegs = keys(%$coords);
271 :    
272 :     my $ali_trees = new AliTrees($fig,$data);
273 :     my $overlaps = {};
274 :     my($peg,@poss,$ali2,$id2,$coords2,%seen,@overlap);
275 :     my($peg1,$c1,$c2,$b1,$e1,$b2,$e2,$ln1,$ln2,$ov);
276 :     foreach $peg (@pegs)
277 :     {
278 :     @poss = $ali_trees->alignments_containing_peg($peg);
279 :     foreach $id2 (@poss)
280 :     {
281 :     if (($id2 ne $id) && (! $seen{$id2}))
282 :     {
283 :     $ali2 = new AliTree($id2,$fig,$data);
284 :     $coords2 = $ali2->coords_of;
285 :     @overlap = ();
286 :     foreach $peg1 (@pegs)
287 :     {
288 :     if (($c1 = $coords->{$peg1}) && ($c2 = $coords2->{$peg1}))
289 :     {
290 :     ($b1,$e1) = @$c1;
291 :     ($b2,$e2) = @$c2;
292 :     $ln1 = ($e1-$b1)+1;
293 :     $ln2 = ($e2-$b2)+1;
294 :     $ov = &FIG::min($e1,$e2) - &FIG::max($b1,$b2);
295 :     if (($ov > (0.9 * $ln1)) && ($ov > (0.9 * $ln2)))
296 :     {
297 :     push(@overlap,[$peg1,$b1,$e1,$b2,$e2]);
298 :     }
299 :     }
300 :     }
301 :    
302 :     if (@overlap > 0)
303 :     {
304 :     $overlaps->{$id2} = [sort { &FIG::by_fig_id($a,$b) } @overlap];
305 :     }
306 :     $seen {$id2} = 1;
307 :     }
308 :     }
309 :     }
310 :     return $overlaps;
311 :     }
312 :    
313 :     sub coords_of {
314 :     my($self,$peg) = @_;
315 :    
316 :     return $peg ? $self->{coords}->{$peg} : $self->{coords};
317 :     }
318 :    
319 :     sub sz {
320 :     my($self) = @_;
321 :     my $n = @{$self->{ali}};
322 :     return $n;
323 :     }
324 :    
325 :    
326 :     sub html {
327 :     my($self) = @_;
328 :    
329 :     my $id = $self->{id};
330 :     my $data = $self->{data};
331 :     my $dir = "$data/$id";
332 :    
333 :     if ((! -s "$dir/full.html") && (-s "$dir/full.ali"))
334 :     {
335 :     &FIG::run("alignment_to_html < $dir/full.ali > $dir/full.html");
336 :     }
337 :     return (-s "$dir/full.html") ? join("",`cat $dir/full.html`) : undef;
338 :     }
339 :    
340 : overbeek 1.2 sub phob_dir {
341 :     my($self) = @_;
342 :    
343 :     my $dir = $self->{dir};
344 :     if (! -d "$dir/PHOB")
345 :     {
346 :     if ($self->make_phob_dir("$dir/PHOB"))
347 :     {
348 :     return "$dir/PHOB";
349 :     }
350 :     else
351 :     {
352 :     return "";
353 :     }
354 :     }
355 :     else
356 :     {
357 :     return "$dir/PHOB";
358 :     }
359 :     }
360 :    
361 :     sub make_phob_dir {
362 :     my($self,$dir) = @_;
363 :    
364 :     if (-d $dir)
365 :     {
366 :     print STDERR "Attempt to make $dir failed: it already exists\n";
367 :     return 0;
368 :     }
369 :     &FIG::verify_dir($dir);
370 :    
371 :     my($ali,$tree);
372 :     if (($ali = $self->ali) && ($tree = $self->tree))
373 :     {
374 :     my($n,$to,$id,$seq,$idN,$tuple);
375 :     open(IDS,">$dir/ids") || die "could not open $dir/ids";
376 :     open(ALI,">$dir/aln.fasta") || die "could not open $dir/aln.fasta";
377 :     open(TREE,">$dir/tree.dnd") || die "could not open $dir/tree.dnd";
378 :     $n = 1;
379 :     $to = {};
380 :     foreach $tuple (@$ali)
381 :     {
382 :     ($id,undef,$seq) = @$tuple;
383 :     $idN = "id$n";
384 :     $n++;
385 :     print IDS "$idN\t$id\n";
386 :     $to->{$id} = $idN;
387 :     $seq =~ s/\s//gs;
388 :     $seq =~ s/[uU]/x/g;
389 :     &FIG::display_id_and_seq($idN,\$seq,\*ALI);
390 :     }
391 :     &tree_utilities::relabel_nodes($tree,$to);
392 :     print TREE &tree_utilities::to_newick($tree),"\n";
393 :     close(IDS);
394 :     close(ALI);
395 :     close(TREE);
396 :     return 1;
397 :     }
398 :     else
399 :     {
400 :     return 0;
401 :     }
402 :     }
403 :    
404 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3