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

Annotation of /FigKernelScripts/CSA_predict_features_based_on_refs.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 use strict;
2 :     #
3 :     # This is a SAS Component
4 :     #
5 :     use Data::Dumper;
6 :     use SeedEnv;
7 :     use gjoseqlib;
8 :     use Carp;
9 :     use FS_RAST;
10 :    
11 :     my $usage = "usage: CSA_predict_features_based_on_refs [-install] Features ContigIndex CurrentDir > recommended.changes";
12 :    
13 :     my $install = $ARGV[0] eq "-install";
14 :     if ($install) { shift @ARGV }
15 :    
16 :     my($featuresF,$indexF,$currentD);
17 :     (
18 :     ($featuresF = shift @ARGV) &&
19 :     ($indexF = shift @ARGV) &&
20 :     ($currentD = shift @ARGV)
21 :     )
22 :     || die $usage;
23 :    
24 :     (-s $featuresF) || die "$featuresF does not exist";
25 :     (-s $indexF) || die "$indexF does not exist";
26 :     (-s "$currentD/Features/peg/tbl") || die "$currentD is not a valid Seed Directory";
27 :    
28 :     my $refH;
29 :     if ($install)
30 :     {
31 :     $refH = &load_dir($currentD);
32 :     }
33 :    
34 :     my $fids = [];
35 :     my $fidH = {};
36 :     &load_existing($currentD,'peg',$fids,$fidH);
37 :     &load_existing($currentD,'rna',$fids,$fidH);
38 :     my @predictions = &load_scored_predictions($featuresF,$indexF);
39 :    
40 :     my @kept;
41 :     foreach my $tuple (@predictions)
42 :     {
43 :     my($id,$contig,$beg,$end,$func,$seq,$errs,$template_peg,undef,undef,$template_translation,$sc) = @$tuple;
44 :     my $existing = 0;
45 :     my $i;
46 :     for ($i=0; ($i < @kept) && (&overlaps($contig,$beg,$end,$kept[$i]->[1],$kept[$i]->[2],$kept[$i]->[3]) < 30); $i++) {}
47 :     if ($i == @kept)
48 :     {
49 :     for ($i=0; ($i < @$fids); $i++)
50 :     {
51 :     my($id1,$contig1,$beg1,$end1) = @{$fids->[$i]};
52 :     if ($contig eq $contig1)
53 :     {
54 :     if (($beg != $beg1) || ($end != $end1))
55 :     {
56 :     if (($_ = &overlaps($contig,$beg,$end,$contig1,$beg1,$end1)) >= 150)
57 :     {
58 :     if ($sc >= 4)
59 :     {
60 :     if ($end == $end1)
61 :     {
62 :     print join("\t",('replace',$id1,$contig1,$beg1,$end1,$contig,$beg,$end,$func,$_,$sc)),"\n";
63 :     $fidH->{join(",",(&SeedUtils::type_of($id1),$contig,$beg,$end))} = 1;
64 :     if ($install)
65 :     {
66 :     &replace($refH,$id1,$contig,$beg,$end,$func,$seq,$errs,$template_peg,$template_translation);
67 :     }
68 :     }
69 :     else
70 :     {
71 :     print join("\t",('delete',$id1,$contig1,$beg1,$end1,$contig,$beg,$end,$func,$_,$sc)),"\n";
72 :     if ($install)
73 :     {
74 :     &delete_fid($refH,$id1);
75 :     }
76 :     }
77 :     }
78 :     else
79 :     {
80 :     $existing = 1;
81 :     }
82 :     }
83 :     }
84 :     }
85 :     }
86 :     my $type = &SeedUtils::type_of($id);
87 :     if ((! $fidH->{join(",",($type,$contig,$beg,$end))}) && (! $existing))
88 :     {
89 :     my $len = abs($end-$beg)+1;
90 :     if (($len > 150) || (! &SeedUtils::hypo($func)))
91 :     {
92 :     print join("\t",('add',$type,$contig,$beg,$end,$len,$func,$errs,$sc,$template_translation)),"\n";
93 :     push(@kept,$tuple);
94 :     if ($install)
95 :     {
96 :     &add($refH,$type,$contig,$beg,$end,$func,$errs,$seq,$template_peg,$template_translation);
97 :     }
98 :     }
99 :     }
100 :     }
101 :     }
102 :    
103 :     if ($install)
104 :     {
105 :     &update_dir($refH,$currentD);
106 :     }
107 :    
108 :     sub update_dir {
109 :     my($refH,$currentD) = @_;
110 :    
111 :     &update_tbl($currentD,$refH,'peg');
112 :     &update_tbl($currentD,$refH,'rna');
113 :     &update_fasta($currentD,$refH,'peg');
114 :     &update_fasta($currentD,$refH,'rna');
115 :     &update_func($currentD,$refH);
116 :     &update_annotations($currentD,$refH);
117 :     }
118 :    
119 :     sub overlaps {
120 :     my($c1,$b1,$e1,$c2,$b2,$e2) = @_;
121 :    
122 :     if ($b1 > $e1) { ($b1,$e1) = ($e1,$b1) }
123 :     if ($b2 > $e2) { ($b2,$e2) = ($e2,$b2) }
124 :    
125 :     if (($b1 <= $b2) && ($e1 >= $b2))
126 :     {
127 :     return (&min($e1,$e2) - $b2) + 1;
128 :     }
129 :     elsif (($b2 <= $b1) && ($e2 >= $b1))
130 :     {
131 :     return (&min($e1,$e2) - $b1) + 1;
132 :     }
133 :     else
134 :     {
135 :     return 0;
136 :     }
137 :     }
138 :    
139 :     sub load_scored_predictions {
140 :     my($features,$index) = @_;
141 :    
142 :     open(INDEX,"<$indexF") || die "could not open $indexF";
143 :     my %to_real = map { ($_ =~ /^(\S+)\t(\S+)$/) ? ($1 => $2) : () } <INDEX>;
144 :     close(INDEX);
145 :    
146 :     open(FEATURES,"<$featuresF") || die "could not open $featuresF";
147 :     my @scored = sort { $b->[11] <=> $a->[11] }
148 :     grep { $_->[11] >= 2 }
149 :     map { chop;
150 :     my @x = split(/\t/,$_);
151 :     $x[1] = $to_real{$x[1]};
152 :     push(@x,&score(&SeedUtils::type_of($x[0]),$x[4],length($x[5]),$x[6],$x[8],$x[9]));
153 :     \@x;
154 :     } <FEATURES>;
155 :     close(FEATURES);
156 :     return @scored;
157 :     }
158 :    
159 :     sub score {
160 :     my($type,$func,$seqlen,$errs,$mapped,$next_best) = @_;
161 :    
162 :     my $score = 2 * ($mapped - $next_best);
163 :     if (! &SeedUtils::hypo($func)) { $score += 2 }
164 :     if ($seqlen > 120) { $score += 2 }
165 :     if ($errs) { $score -= 10 }
166 :     return $score;
167 :     }
168 :    
169 :     sub load_existing {
170 :     my($seed_dir,$type,$fids,$fidH) = @_;
171 :     open(TBL,"<$seed_dir/Features/$type/tbl")
172 :     || die "$seed_dir/Features/$type/tbl does not exist";
173 :     foreach $_ (<TBL>)
174 :     {
175 :     if ($_ =~ /^(fig\|\d+\.\d+\.$type\.\d+)\t(\S+)_(\d+)_(\d+)/)
176 :     {
177 :     push(@$fids,[$1,$2,$3,$4]);
178 :     $fidH->{join(",",($type,$2,$3,$4))} = $1;
179 :     }
180 :     }
181 :     close(TBL);
182 :     }
183 :    
184 :     sub load_dir {
185 :     my($currentD) = @_;
186 :    
187 :     my $refH = {};
188 :     $refH->{'deleted'}->{peg} = &deleted_fids($currentD,'peg');
189 :     $refH->{'deleted'}->{rna} = &deleted_fids($currentD,'rna');
190 :     $refH->{annotations} = [];
191 :     &load_tbl($currentD,$refH,'rna');
192 :     &load_tbl($currentD,$refH,'peg');
193 :     $refH->{max}->{rna} = &max_fid($refH->{tbl}->{rna});
194 :     $refH->{max}->{peg} = &max_fid($refH->{tbl}->{peg});
195 :     &load_fasta($currentD,$refH,'rna');
196 :     &load_fasta($currentD,$refH,'peg');
197 :     &load_func($currentD,$refH);
198 :     $refH->{code} = 11;
199 :     return $refH;
200 :     }
201 :    
202 :     sub max_fid {
203 :     my($tblH) = @_;
204 :    
205 :     my @fids = sort { &SeedUtils::by_fig_id($a,$b) } keys(%$tblH);
206 :     return $fids[-1];
207 :     }
208 :    
209 :     sub deleted_fids {
210 :     my($currentD,$type) = @_;
211 :    
212 :     my $delH = {};
213 :     if (open(DEL,"<$currentD/Features/$type/delete.features"))
214 :     {
215 :     while (defined($_ = <DEL>))
216 :     {
217 :     if ($_ =~ /^(\S+)/) { $delH->{$1} = 1 }
218 :     }
219 :     close(DEL);
220 :     }
221 :     return $delH;
222 :     }
223 :    
224 :     sub load_tbl {
225 :     my($currentD,$refH,$type) = @_;
226 :    
227 :     my $file = "$currentD/Features/$type/tbl";
228 :     my $delH = $refH->{deleted}->{$type};
229 :     open(TBL,"<$file") || die "cannot open $file";
230 :     my %tbl = map { (($_ =~ /^(\S+)/) && (! $delH->{$1})) ? ($1 => $_) : () } <TBL>;
231 :     close(TBL);
232 :     $refH->{tbl}->{$type} = \%tbl;
233 :     }
234 :    
235 :     sub update_tbl {
236 :     my($currentD,$refH,$type) = @_;
237 :    
238 :     my $tbl = $refH->{tbl}->{$type};
239 :     my $file = "$currentD/Features/$type/tbl";
240 :     if (-s "$file~") { unlink("$file~") }
241 :     rename($file,"$file~");
242 :     open(TBL,">$file") || die "could not open $file";
243 :     foreach my $fid (sort { &SeedUtils::by_fig_id($a,$b) } keys(%$tbl))
244 :     {
245 :     print TBL $tbl->{$fid};
246 :     }
247 :     close(TBL);
248 :     my $delH = $refH->{deleted}->{$type};
249 :     if ($refH)
250 :     {
251 :     my @deleted = sort { &SeedUtils::by_fig_id($a,$b) } keys(%$delH);
252 :     if (@deleted > 0)
253 :     {
254 :     my $file = "$currentD/Features/$type/deleted.features";
255 :     open(DEL,">$file") || die "could not open $file";
256 :     foreach $_ (@deleted)
257 :     {
258 :     print DEL "$_\n";
259 :     }
260 :     close(DEL);
261 :     }
262 :     }
263 :     }
264 :    
265 :     sub load_fasta {
266 :     my($currentD,$refH,$type) = @_;
267 :    
268 :     my $delH = $refH->{deleted}->{$type};
269 :     my $file = "$currentD/Features/$type/fasta";
270 :     my @seqs = grep { ! $delH->{$_->[0]} } &gjoseqlib::read_fasta($file);
271 :     my %seqsH = map { $_->[0] => [$_->[1],$_->[2]] } @seqs;
272 :     $refH->{fasta}->{$type} = \%seqsH;
273 :     }
274 :    
275 :     sub update_fasta {
276 :     my($currentD,$refH,$type) = @_;
277 :    
278 :     my $file = "$currentD/Features/$type/fasta";
279 :     if (-s "$file~") { unlink("$file~") }
280 :     rename($file,"$file~");
281 :     my $seqH = $refH->{fasta}->{$type};
282 :     my @seqs = map { [$_,@{$seqH->{$_}}] }
283 :     sort { &SeedUtils::by_fig_id($a,$b) }
284 :     keys(%$seqH);
285 :     &gjoseqlib::print_alignment_as_fasta($file,\@seqs);
286 :     }
287 :    
288 :     sub load_func {
289 :     my($currentD,$refH) = @_;
290 :    
291 :     my $assignments = {};
292 :     &load_file($assignments,"$currentD/proposed_non_ff_functions");
293 :     &load_file($assignments,"$currentD/proposed_functions");
294 :     &load_file($assignments,"$currentD/assigned_functions");
295 :     $refH->{assignments} = $assignments;
296 :     }
297 :    
298 :     sub load_file {
299 :     my($assignments,$file) = @_;
300 :    
301 :     if (open(FILE,"<",$file))
302 :     {
303 :     while (defined($_ = <FILE>))
304 :     {
305 :     chop;
306 :     my($fid,$func) = split(/\t/,$_);
307 :     $assignments->{$fid} = $func;
308 :     }
309 :     close(FILE);
310 :     }
311 :     }
312 :    
313 :     sub update_func {
314 :     my($currentD,$refH) = @_;
315 :    
316 :     my $file = "$currentD/assigned_functions";
317 :     if (-s "$file~") { unlink("$file~") }
318 :     rename($file,"$file~");
319 :     open(ASSIGN,">$file") || die "could not open $file";
320 :     my $assignments = $refH->{assignments};
321 :     foreach my $fid (sort { &SeedUtils::by_fig_id($a,$b) } keys(%$assignments))
322 :     {
323 :     print ASSIGN join("\t",($fid,$assignments->{$fid})),"\n";
324 :     }
325 :     close(ASSIGN);
326 :     }
327 :    
328 :     sub update_annotations {
329 :     my($currentD,$refH) = @_;
330 :    
331 :     my $annotations = $refH->{annotations};
332 :     open(ANN,">>$currentD/annotations")
333 :     || confess "could not open $currentD/annotations";
334 :    
335 :     my $time_made = time;
336 :     foreach my $annotation (@$annotations) {
337 :     my ($fid,$text) = @$annotation;
338 :     print ANN "$fid\n$time_made\nchanges to sequence\n$text\n//\n";
339 :     }
340 :     close(ANN);
341 :     }
342 :    
343 :     sub replace {
344 :     my($refH,$fid,$contig,$beg,$end,$func,$seq,$errs,$template_peg,$template_translation) = @_;
345 :    
346 :     &delete_fid($refH,$fid);
347 :     my $type = &SeedUtils::type_of($fid);
348 :     &add($refH,$type,$contig,$beg,$end,$func,$errs,$seq,$template_peg,$template_translation,$fid);
349 :     }
350 :    
351 :     sub delete_fid {
352 :     my($refH,$fid) = @_;
353 :    
354 :     my $type = &SeedUtils::type_of($fid);
355 :     $refH->{deleted}->{$type}->{$fid} = 1;
356 :     delete $refH->{tbl}->{$type}->{$fid};
357 :     delete $refH->{fasta}->{$type}->{$fid};
358 :     }
359 :    
360 :     sub add {
361 :     my($refH,$type,$contig,$beg,$end,$func,$errs,$seq,$template_peg,$template_translation,$fid) = @_;
362 :    
363 :     if (! $fid)
364 :     {
365 :     my $last_fid = $refH->{max}->{$type};
366 :     ($last_fid =~ /^(fig\|\d+\.\d+\.[^\.]+\.)(\d+)/) || die "bad fid: $last_fid";
367 :     $fid = $1 . ($2 + 1);
368 :     $refH->{max}->{$type} = $fid;
369 :     }
370 :     if ($refH->{deleted}->{$type}->{$fid}) { delete $refH->{deleted}->{$type}->{$fid} }
371 :     my $tbl_line = join("\t",($fid,join("_",($contig,$beg,$end)),'')) . "\n";
372 :     $refH->{tbl}->{$type}->{$fid} = $tbl_line;
373 :     $func || ($func = '');
374 :     $refH->{assignments}->{$fid} = $func;
375 :    
376 :     my $prot_seq;
377 :     if ($type eq 'peg')
378 :     {
379 :     if (! $errs)
380 :     {
381 :     ### we need to construct the correct genetic code here, if it is not 11 ###
382 :     if ($refH->{code} == 11)
383 :     {
384 :     $prot_seq = &SeedUtils::translate($seq,undef,1);
385 :     $prot_seq =~ s/\*$//;
386 :     }
387 :     else
388 :     {
389 :     die "we do not support genetic code $refH->{code} yet";
390 :     }
391 :     }
392 :     else
393 :     {
394 :     my $params = {};
395 :     $params->{family} = [[ $template_peg, "", $template_translation ]];
396 :     $params->{code} = $refH->{code};
397 :     my ($new_loc, $new_translation, undef, $annotation) = &FS_RAST::best_match_in_family($params, [$contig,$beg,$end,$seq]);
398 :     if ($new_loc)
399 :     {
400 :     $prot_seq = $new_translation;
401 :     push(@{$refH->{annotations}},[$fid,$annotation]);
402 :     }
403 :     }
404 :     }
405 :     $refH->{fasta}->{$type}->{$fid} = ['',($type eq 'peg') ? $prot_seq : $seq];
406 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3