[Bio] / FigKernelScripts / insert.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/insert.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.2 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : overbeek 1.1 # usage: insert Ali Tree Id [Weights] > new_tree
19 :    
20 :     use tree_utilities;
21 :     use Carp;
22 :     use Data::Dumper;
23 :     use fastDNAml;
24 :    
25 :     $usage = "usage: insert Ali Tree Id Weights > new_tree";
26 :    
27 :     $size_rep = 50;
28 :    
29 :     ($ali = shift @ARGV)
30 :     || die "bad alignment: $usage";
31 :    
32 :     if ((! -s "$ali.index") || ((-M $ali) < (-M "$ali.index")))
33 :     {
34 :     (system("index_fasta $ali > $ali.index") == 0)
35 :     || die "could not make $ali.index";
36 :     }
37 :     open(INDEX,"<$ali.index")
38 :     || die "could not open $ali.index";
39 :     while (defined($_ = <INDEX>))
40 :     {
41 :     chop;
42 :     ($id,$seek,$ln) = split(/\t/,$_);
43 :     $seek_of{$id} = [$seek,$ln];
44 :     }
45 :     close(INDEX);
46 :     open(ALI,"<$ali") || die "could not open $ali";
47 :    
48 :     (($tree_file = shift @ARGV) && (-s $tree_file))
49 :     || die "missing tree: $usage";
50 :    
51 :     ($id_to_ins = shift @ARGV)
52 :     || die "missing id to insert: $usage";
53 :    
54 :     if (!($weights = shift @ARGV)) { $weights = ""; }
55 :    
56 :     ($tree = &uproot(&parse_newick_tree(join("",`cat \"$tree_file\"`))))
57 :     || die "could not parse the tree";
58 :     $big_indexes = &tree_index_tables($tree);
59 :    
60 :     ############## get the tree ##########
61 :     $tips = &tips_of_tree($tree);
62 :     @in_tree = ();
63 :     foreach $x (@$tips)
64 :     {
65 :     $in_tree{$x} = 1;
66 :     }
67 :    
68 :    
69 :     if ($in_tree{$id_to_ins})
70 :     {
71 :     die "$id_to_ins is already in the tree";
72 :     }
73 :    
74 :     $bad = 0;
75 :     foreach $id (@$tips)
76 :     {
77 :     if ((! &seq_of($id,\%seek_of,\*ALI)) && (! &seq_of("\'$id\'",\%seek_of,\*ALI)))
78 :     {
79 :     $bad = 1;
80 :     warn "The tree contains $id, but it is not in the alignment\n";
81 :     }
82 :     }
83 :     (! $bad) || die "Could not insert - alignment is missing sequences";
84 :    
85 :     print STDERR "insert: best_hits $id_to_ins $ali\n";
86 :     @closest = `best_hits $id_to_ins $ali`;
87 :     chop @closest;
88 :     foreach $x (@closest)
89 :     {
90 :     $x = &printable_to_label($x);
91 :     $in{$x} = 1;
92 :     }
93 :    
94 :     $rep_tree = &representative_by_size($tree,$size_rep);
95 :     $rep_tips = &tips_of_tree($rep_tree);
96 :     foreach $x (@$rep_tips)
97 :     {
98 :     $in{$x} = 1;
99 :     }
100 :    
101 :     $subtree = &subtree($tree,\%in);
102 :     $x = &display_tree($subtree,{});
103 :     print STDERR "Initial Representative Tree:\n\n$x\n";
104 :    
105 :     ##################################################
106 :     @seen = ();
107 :    
108 :     while (! &seen($subtree))
109 :     {
110 :     ($where,$len_to_add) = &sequence_ins_point($subtree,$id_to_ins,\%seek_of,\*ALI);
111 :     # print STDERR "node1=@{$where->[0]} node2=@{$where->[1]} frac=$where->[2] len_to_add=$len_to_add\n";
112 :     $subtree = &neighborhood_of_tree_point($big_indexes,$where,$size_rep);
113 :     $x = &display_tree($subtree,{});
114 :     # print STDERR "new neighborhood\n";
115 :     # print STDERR $x;
116 :     }
117 :     $new_tree = &to_newick(&add_branch($big_indexes,$where,$len_to_add,$id_to_ins));
118 :     print "$new_tree\n";
119 :     print STDERR "successfully inserted $id_to_ins\n";
120 :    
121 :     sub seen {
122 :     my($tree) = @_;
123 :     my($tips,@sorted,$key);
124 :    
125 :     $tips = &tips_of_tree($tree);
126 :     @sorted = sort @$tips;
127 :     $key = join(" ",@sorted);
128 :    
129 :     if (! $seen{$key})
130 :     {
131 :     $seen{$key} = 1;
132 :     return 0;
133 :     }
134 :     return 1;
135 :     }
136 :    
137 :     sub add_branch {
138 :     my($tabP,$where,$len_to_add,$id) = @_;
139 :     my($nd1,$nd2,$node1,$node2,$fract,$node,$desc);
140 :    
141 :     ($node1,$node2,$fract) = @$where;
142 :     $nd1 = &locate_node($node1,$tabP);
143 :     $nd2 = &locate_node($node2,$tabP);
144 :     $tree = &root_tree_between_nodes($nd1,$nd2,$fract,$tabP);
145 :     $node = [$id,$len_to_add,[$tree]];
146 :     $desc = $tree->[2];
147 :     push(@$desc,$node);
148 :     return $tree;
149 :     }
150 :    
151 :     sub sequence_ins_point {
152 :     my($tree,$id,$seek_of,$ali_fh) = @_;
153 :     my($id_seqs,$tips,$seq1,@weights_and_rates,@weights,$column_weights,$ncats);
154 :     my(@cats,$cats,@lines,$new_trees,$x,$id1);
155 :     my($indexes,$at_node,$sz1,$sz2,$len_to_add,$where);
156 :     my($tips1,$tips2,$node1,$node2,$d1,$d2,$frac);
157 :    
158 :     $id_seqs = [];
159 :     $tips = &tips_of_tree($tree);
160 :     foreach $id1 (@$tips)
161 :     {
162 :     $seq1 = &seq_of($id1,$seek_of,$ali_fh);
163 :     $seq1 || confess "could not locate sequence for $id1";
164 :     push(@$id_seqs,[$id1,$seq1]);
165 :     }
166 :     $seq1 = &seq_of($id,\%seek_of,\*ALI);
167 :     $seq1 || confess "could not locate sequence for $id";
168 :     push(@$id_seqs,[$id,$seq1]);
169 :    
170 :     if ($weights)
171 :     {
172 :     @weights_and_rates = `cat $weights`;
173 :     @weights = ();
174 :     while (($_ = shift @weights_and_rates) && ($_ !~ /^C/))
175 :     {
176 :     push(@weights,$_);
177 :     }
178 :     $column_weights = join("",@weights);
179 :     $column_weights =~ s/^Weights//;
180 :     $column_weights =~ s/\s//g;
181 :    
182 :     if ($_ =~ /^C\s+(\d+)\s+(\S.*\S)\s*$/)
183 :     {
184 :     $ncats = $1;
185 :     @cats = ($2);
186 :     while (($_ = shift @weights_and_rates) && ($_ !~ /^Categories/))
187 :     {
188 :     push(@cats,$_);
189 :     }
190 :    
191 :     $cats = join("",@cats);
192 :     $cats =~ s/\s+$//;
193 :     $cats =~ s/^\s+//;
194 :     @cats = split(/\s+/,$cats);
195 :     (@cats == $ncats) || die "category mismatch - BUG @cats";
196 :    
197 :     $_ =~ s/^Categories\s+//;
198 :     $_ =~ s/\s//g;
199 :     @lines = ($_);
200 :     while ($_ = shift @weights_and_rates)
201 :     {
202 :     $_ =~ s/\s//g;
203 :     push(@lines,$_);
204 :     }
205 :     $column_cats = join("",@lines);
206 :     }
207 :     else
208 :     {
209 :     die "bad weights and rates file: missing category data";
210 :     }
211 :    
212 :     $new_trees = &fastDNAml($id_seqs,
213 :     [
214 :     ["categories",\@cats,$column_cats],
215 :     ["weights",$column_weights],
216 :     ["fast_add"],
217 :     ["restart",$subtree],
218 :     ["global",0,0]
219 :     ]);
220 :     }
221 :     else
222 :     {
223 :     $new_trees = &fastDNAml($id_seqs,[["fast_add"],["restart",$subtree],["global",0,0]]);
224 :     }
225 :    
226 :     (@$new_trees == 1) || die "Bad value back from fastDNAml";
227 :     $new_tree = $new_trees->[0]->[1];
228 :     $x = &display_tree($new_tree,{});
229 :     print STDERR "Tree returned by fastDNAml\n";
230 :     print STDERR $x;
231 :    
232 :     $indexes = &tree_index_tables($new_tree);
233 :     $at_node = &root_tree_at_node($indexes,$id);
234 :    
235 :     # &print_tree($at_node);
236 :    
237 :     $sz1 = &tree_length($tree);
238 :     $sz2 = &tree_length($at_node->[2]->[1]) - $at_node->[2]->[1]->[1];
239 :     if (($sz1 >= 0) && ($sz2 > 0))
240 :     {
241 :     $len_to_add = $at_node->[2]->[1]->[1] * ($sz1 / $sz2);
242 :     }
243 :     else
244 :     {
245 :     $len_to_add = 0;
246 :     }
247 :     # print STDERR "sz1=$sz1 sz2=$sz2 len_to_add=$len_to_add\n";
248 :     $descendant1 = $at_node->[2]->[1]->[2]->[1];
249 :     $descendant2 = $at_node->[2]->[1]->[2]->[2];
250 :     $node1 = &rooted_identify($descendant1);
251 :     $node2 = &rooted_identify($descendant2);
252 :     if (@$node1 == 2) { push(@$node1,$node2->[0]); }
253 :     if (@$node2 == 2) { push(@$node2,$node1->[0]); }
254 :    
255 :     $d1 = $descendant1->[1];
256 :     $d2 = $descendant2->[1];
257 :     if (($d1 + $d2 ) > 0)
258 :     {
259 :     $frac = ($d1 / ($d1+$d2));
260 :     }
261 :     else
262 :     {
263 :     $frac = 0.5;
264 :     }
265 :     # print STDERR "node1=@$node1 node2=@$node2 frac=$frac\n";
266 :     $where = [$node1,$node2,$frac];
267 :     #
268 :     # Set $where
269 :     #
270 :     return ($where,$len_to_add);
271 :     }
272 :    
273 :     sub tree_length {
274 :     my($tree) = @_;
275 :     my($desc,$ln,$i);
276 :    
277 :     $desc = $tree->[2];
278 :     for ($ln=0,$i=1; ($i <= $#{$desc}); $i++)
279 :     {
280 :     $ln += &tree_length($desc->[$i]) + $desc->[$i]->[1];
281 :     }
282 :     return $ln;
283 :     }
284 :    
285 :     sub rooted_identify {
286 :     my($tree) = @_;
287 :    
288 :     $cc = $tree->[2];
289 :     if (@$cc == 1) { return [$tree->[0]]; }
290 :     if (@$cc == 2) { return ""; }
291 :     return [&node_of_tree($cc->[1]),&node_of_tree($cc->[2])];
292 :     }
293 :    
294 :     sub node_of_tree {
295 :     my($tree) = @_;
296 :    
297 :     while ($#{$tree->[2]} > 0)
298 :     {
299 :     $tree = $tree->[2]->[1];
300 :     }
301 :     return $tree->[0];
302 :     }
303 :    
304 :     sub seq_of {
305 :     my($id,$seek_of,$ali_fh) = @_;
306 :     my($x,$seek,$len,$seq);
307 :    
308 :     if ($x = $seek_of->{$id})
309 :     {
310 :     ($seek,$len) = @$x;
311 :     seek($ali_fh,$seek,0);
312 :     read($ali_fh,$seq,$len);
313 :     $seq =~ s/\s//g;
314 :     $seq =~ tr/a-z/A-Z/;
315 :     $seq =~ tr/U/T/;
316 :     $seq =~ s/[~.]/\?/g;
317 :     return $seq;
318 :     }
319 :     die "invalid seq id=$id";
320 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3