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

Annotation of /FigKernelPackages/fastDNAml.pm

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 package fastDNAml;
19 :    
20 :     use tree_utilities;
21 :     use Carp;
22 :     require Exporter;
23 :     @ISA = (Exporter);
24 :     @EXPORT = qw(
25 :     fastDNAml
26 :     );
27 :    
28 :     # Example:
29 :     # $ans = fastDNAml([["id1","acgtacgt"],["id2","acgt-cgt"]],[[jumble,13],
30 :     # [weights,"01010101"],
31 :     # ["categories","ABCD1234"]],
32 :     #
33 :     # The answer is a list of [LogLikelihood,Tree] pairs
34 :     #
35 :    
36 :     %codes = (
37 :     "bootstrap" => "B",
38 :     "categories" => "C",
39 :     "fast_add" => "Q",
40 :     "frequencies" => "F",
41 :     "global" => "G",
42 :     "jumble" => "J",
43 :     "restart" => "R",
44 :     "treefile" => "Y",
45 :     "transition" => "T",
46 :     "userlengths" => "U",
47 :     "weights" => "W",
48 :     );
49 :    
50 :    
51 :     sub fastDNAml {
52 :     my($id_seqs,$options) = @_;
53 :     my($tmp,$functor,$arity,$opt,$code,$num_seqs,$len_seqs);
54 :     my($n,$rate,$rates,$i,$x,$tree);
55 :    
56 :     $added_n = 1;
57 :    
58 :     @translated = ();
59 :     foreach $pair (@$id_seqs)
60 :     {
61 :     $id = $pair->[0];
62 :     if ((index($id,"_") >= 0) || (length($id) > 10))
63 :     {
64 :     $id1 = "ZZadded$added_n";
65 :     $added_n++;
66 :     $pair->[0] = $id1;
67 :     push(@translated,[$id,$id1]);
68 :     # print STDERR "$id -> $id1\n";
69 :     }
70 :     }
71 :    
72 :     &fix_options($options);
73 :     mkdir("Tmp$$",0777);
74 :     $num_seqs = @$id_seqs;
75 :     $tmp = $id_seqs->[0]->[1];
76 :     $tmp =~ s/\s//g;
77 :     $len_seqs = length($tmp);
78 :     $tmp = "Tmp$$/tmp$$";
79 :    
80 :     open(TMP,">$tmp") || die "could not open $tmp";
81 :     print TMP "$num_seqs $len_seqs ";
82 :     foreach $opt (@$options)
83 :     {
84 :     $functor = $opt->[0];
85 :     $arity = $#{$opt} - 1;
86 :     if (($code = $codes{$functor}) &&
87 :     (! (($functor eq "frequencies") && ($arity == 4))))
88 :     {
89 :     print TMP " $code";
90 :     }
91 :     }
92 :     print TMP "\n";
93 :    
94 :     foreach $opt (@$options)
95 :     {
96 :     $functor = $opt->[0];
97 :     $arity = $#{$opt};
98 :     if (($functor eq "bootstrap") && ($arity == 1))
99 :     {
100 :     print TMP "B $opt->[1]\n";
101 :     }
102 :     elsif (($functor eq "global") && ($arity == 1))
103 :     {
104 :     print TMP "G $opt->[1]\n";
105 :     }
106 :     elsif (($functor eq "global") && ($arity == 2))
107 :     {
108 :     print TMP "G $opt->[1] $opt->[2]\n";
109 :     }
110 :     elsif (($functor eq "jumble") && ($arity == 1))
111 :     {
112 :     print TMP "J $opt->[1]\n";
113 :     }
114 :     elsif (($functor eq "transition") && ($arity == 1))
115 :     {
116 :     print TMP "T $opt->[1]\n";
117 :     }
118 :     elsif (($functor eq "treefile") && ($arity == 1))
119 :     {
120 :     print TMP "Y $opt->[1]\n";
121 :     }
122 :     elsif (($functor eq "categories") && ($arity == 2))
123 :     {
124 :     $n = $#{$opt->[1]} + 1;
125 :     print TMP "C $n";
126 :     $rates = $opt->[1];
127 :     foreach $rate (@$rates)
128 :     {
129 :     print TMP " $rate";
130 :     }
131 :     print TMP "\n";
132 :     print TMP "Categories $opt->[2]\n";
133 :     }
134 :     elsif (($functor eq "weights") && ($arity == 1))
135 :     {
136 :     print TMP "Weights $opt->[1]\n";
137 :     }
138 :     elsif (($functor eq "frequencies") && ($arity == 4))
139 :     {
140 :     print TMP "$opt->[1] $opt->[2] $opt->[3] $opt->[4]\n";
141 :     }
142 :     }
143 :     foreach $pair (@$id_seqs)
144 :     {
145 :     ($id,$seq) = @$pair;
146 :     $pad = " " x (13 - length($id));
147 :     print TMP "$id$pad$seq\n";
148 :     # print STDERR "seq: $id\n";
149 :     }
150 :    
151 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne "usertree"); $i++) {}
152 :     if ($i <= $#{$options})
153 :     {
154 :     $trees = $options->[$i]->[1];
155 :     $n = $#{$trees} + 1;
156 :     print TMP "$n\n";
157 :     foreach $tree (@$trees)
158 :     {
159 :     $x = &to_translated_newick($tree);
160 :     print TMP "$x\n";
161 :     }
162 :     }
163 :    
164 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne "restart"); $i++) {}
165 :     if ($i <= $#{$options})
166 :     {
167 :     $tree = $options->[$i]->[1];
168 :     $x = &to_translated_newick($tree);
169 :     print TMP "$x\n";
170 :     }
171 :     close(TMP);
172 :     system "cd Tmp$$; fastDNAml < tmp$$ > output; cd ..";
173 :    
174 :     @trees = glob("Tmp$$/treefile*");
175 :     if (@trees != 1)
176 :     {
177 :     confess "no tree produced by fastDNAml - look in Tmp$$";
178 :     }
179 :     system "cd Tmp$$; mv treefile.* fastDNAml_tree; cd ..";
180 :     @trees = `cat Tmp$$/fastDNAml_tree`;
181 :     $trees = [];
182 :     foreach $tree (@trees)
183 :     {
184 :     if ($tree =~ /^.*\[.*likelihood\s*=\s*([0-9.,e-]+),.*\]\s*(\S.*\S)\s*$/)
185 :     {
186 :     $lk = $1; $just_tree = $2;
187 :     $just_tree =~ s/\'//g;
188 :     foreach $pair (@translated)
189 :     {
190 :     $just_tree =~ s/\b$pair->[1]\b/$pair->[0]/g;
191 :     }
192 :     $parsed = &parse_newick_tree($just_tree);
193 :     push(@$trees,[$lk,$parsed]);
194 :     }
195 :     }
196 :     system "/bin/rm -r Tmp$$";
197 :     return $trees;
198 :     }
199 :    
200 :     sub to_translated_newick {
201 :     my($tree) = @_;
202 :     my($x,$pair,$y);
203 :    
204 :     $x = &to_newick($tree);
205 :     foreach $pair (@translated)
206 :     {
207 :     $y = quotemeta($pair->[0]);
208 :     $x =~ s/\b$y\b/$pair->[1]/g;
209 :     }
210 :     return $x;
211 :     }
212 :    
213 :     sub memberchk {
214 :     my($opt,$options) = @_;
215 :     my($i);
216 :    
217 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne $opt); $i++) {}
218 :     return ($i <= $#{$options});
219 :     }
220 :    
221 :     sub fix_options {
222 :     my($options) = @_;
223 :     my($i,$x);
224 :    
225 :     if (! (&memberchk("bootstrap",$options) ||
226 :     &memberchk("jumble",$options) ||
227 :     &memberchk("transition",$options) ||
228 :     &memberchk("weights",$options) ||
229 :     &memberchk("categories",$options)))
230 :     {
231 :     unshift(@$options,["transition",2.0]);
232 :     }
233 :    
234 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne "frequencies"); $i++) {}
235 :     if ($i <= $#{$options})
236 :     {
237 :     $x = $options->[$i];
238 :     splice(@$options,$i,1);
239 :     push(@$options,$x);
240 :     }
241 :     else
242 :     {
243 :     push(@$options,["frequencies"]);
244 :     }
245 :    
246 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne "global"); $i++) {}
247 :     if ($i <= $#{$options})
248 :     {
249 :     $x = $options->[$i];
250 :     splice(@$options,$i,1);
251 :     unshift(@$options,$x);
252 :     }
253 :    
254 :     for ($i=0; ($i <= $#{$options}) && ($options->[$i]->[0] ne "treefile"); $i++) {}
255 :     if ($i <= $#{$options})
256 :     {
257 :     $x = $options->[$i];
258 :     splice(@$options,$i,1);
259 :     unshift(@$options,$x);
260 :     }
261 :     else
262 :     {
263 :     unshift(@$options,["treefile"]);
264 :     }
265 :     }
266 :    
267 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3