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

Annotation of /FigKernelScripts/insert.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3