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

Annotation of /FigKernelPackages/fastDNAml.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3