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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3