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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.18 # -*- Perl -*-
2 :    
3 : overbeek 1.14 ########################################################################
4 : olson 1.12 #
5 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
6 :     # for Interpretations of Genomes. All Rights Reserved.
7 :     #
8 :     # This file is part of the SEED Toolkit.
9 :     #
10 :     # The SEED Toolkit is free software. You can redistribute
11 :     # it and/or modify it under the terms of the SEED Toolkit
12 :     # Public License.
13 :     #
14 :     # You should have received a copy of the SEED Toolkit Public License
15 :     # along with this program; if not write to the University of Chicago
16 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
17 :     # Genomes at veronika@thefig.info or download a copy from
18 :     # http://www.theseed.org/LICENSE.TXT.
19 :     #
20 : gdpusch 1.18 ########################################################################
21 : olson 1.12
22 : efrank 1.1 use FIG;
23 : gdpusch 1.16 use Carp;
24 : gdpusch 1.18 use Data::Dumper;
25 :    
26 :     my $fig = FIG->new();
27 :     my %trans_table;
28 :    
29 :     my $usage = qq(usage: parse_genbank [-bioname=GENOME_BIONAME] [-taxonomy=TAXONOMY] [-source=DATA_SOURCE] Taxon_ID Org_Dir < genbank.entry);
30 : efrank 1.1
31 : gdpusch 1.18 my $trouble = 0;
32 :     my $force_bioname;
33 :     my $force_taxonomy;
34 : gdpusch 1.22 my $data_source = qq(Parsed GenBank File);
35 : gdpusch 1.9 while ($ARGV[0] =~ m/^-/)
36 :     {
37 : gdpusch 1.23 if ($ARGV[0] =~ m/^-{1,2}bioname=(\S+)/) {
38 :     $force_bioname = $1;
39 :     $force_bioname =~ s/^[\'\"]//o;
40 :     $force_bioname =~ s/[\'\"]$//o;
41 : gdpusch 1.18 $force_bioname =~ s/\s+/ /g;
42 :    
43 :     print STDERR qq(\nGenome bioname will be taken as: \"$force_bioname\"\n);
44 :     }
45 : gdpusch 1.23 elsif ($ARGV[0] =~ m/^-{1,2}taxonomy=(\S+)/) {
46 :     $force_taxonomy = $1;
47 :     $force_taxonomy =~ s/^[\'\"]//o;
48 :     $force_taxonomy =~ s/[\'\"]$//o;
49 : gdpusch 1.18 $force_taxonomy =~ s/\s+/ /g;
50 :     if ($force_taxonomy !~ m/\.$/) { $force_taxonomy .= qq(.); }
51 :    
52 :     print STDERR qq(\nTaxonomy will be taken as: \"$force_taxonomy\"\n);
53 :     }
54 : gdpusch 1.23 elsif ($ARGV[0] =~ m/^-{1,2}source=(\S+)/) {
55 :     $data_source = $1;
56 :     $data_source =~ s/^[\'\"]//o;
57 :     $data_source =~ s/[\'\"]$//o;
58 : gdpusch 1.18
59 :     print STDERR qq(Data source string is: \"$data_source\"\n);
60 :     }
61 :     else {
62 :     $trouble = 1;
63 :     die qq(Could not handle $ARGV[0]);
64 : gdpusch 1.9 }
65 : gdpusch 1.23 shift @ARGV;
66 : gdpusch 1.9 }
67 :    
68 : gdpusch 1.18 my $taxon_ID;
69 :     my $org_dir;
70 : efrank 1.2 (
71 : gdpusch 1.18 ($taxon_ID = shift(@ARGV)) &&
72 :     ($org_dir = shift(@ARGV))
73 : efrank 1.1 )
74 :     || die $usage;
75 :    
76 : gdpusch 1.18 $prefixP = qq(fig|$taxon_ID.peg.);
77 :     $prefixR = qq(fig|$taxon_ID.rna.);
78 :    
79 :     #...Create paths down to PEG and RNA dirs
80 :     &FIG::verify_dir("$org_dir/Features/peg");
81 :     &FIG::verify_dir("$org_dir/Features/rna");
82 :    
83 :     my $fh_contigs;
84 :     my $contigs_file = qq($org_dir/contigs);
85 :     open($fh_contigs, qq(>$contigs_file))
86 :     || die qq(could not open $contigs_file);
87 :    
88 :     my $fh_assigned_funcs;
89 :     my $assigned_funcs_file = qq($org_dir/assigned_functions);
90 :     open($fh_assigned_funcs, qq(>$assigned_funcs_file))
91 :     || die qq(could not write-open $assigned_funcs_file);
92 :    
93 :     my $fh_annotations;
94 :     my $annotations_file = qq($org_dir/annotations);
95 :     open($fh_annotations, qq(>$annotations_file))
96 :     || die qq(could not write-open $annotations_file);
97 :    
98 :     my $fh_ec_nums;
99 :     my $ec_nums_file = qq($org_dir/EC_numbers);
100 :     open($fh_ec_nums, qq(>$ec_nums_file))
101 :     || die" Could not write-open $ec_nums_file";
102 :    
103 :     my $fh_peg_tbl;
104 :     my $peg_tbl_file = qq($org_dir/Features/peg/tbl);
105 :     open($fh_peg_tbl, qq(>$peg_tbl_file))
106 :     || die qq(could not open $peg_tbl_file);
107 :    
108 :     my $fh_peg_fasta;
109 :     my $peg_fasta_file = qq($org_dir/Features/peg/fasta);
110 :     open($fh_peg_fasta, qq(>$peg_fasta_file))
111 :     || die qq(could not open $peg_fasta_file);
112 :    
113 :     my $fh_rna_tbl;
114 :     my $rna_tbl_file = qq($org_dir/Features/rna/tbl);
115 :     open($fh_rna_tbl, qq(>$rna_tbl_file))
116 :     || die qq(could not write-open $rna_tbl_file);
117 :    
118 :     my $fh_rna_fasta;
119 :     my $rna_fasta_file = qq($org_dir/Features/rna/fasta);
120 :     open($fh_rna_fasta, qq(>$rna_fasta_file))
121 :     || die qq(could not write-open $rna_fasta_file);
122 : efrank 1.2
123 : gdpusch 1.18
124 :     open( TAXONOMY, qq(>$org_dir/TAXONOMY)) || die qq(could not open $org_dir/TAXONOMY);
125 :     open( GENOME, qq(>$org_dir/GENOME)) || die qq(could not open $org_dir/GENOME);
126 :    
127 :     open( PROJECT, qq(>$org_dir/PROJECT)) || die qq(could not open $org_dir/PROJECT);
128 :     print PROJECT qq($data_source\n);
129 : overbeek 1.4 close(PROJECT);
130 :    
131 : gdpusch 1.18 my $idNp = 1;
132 :     my $idNr = 1;
133 : overbeek 1.4
134 : gdpusch 1.18 $/ = qq(\n//\n);
135 :     my $record;
136 :     my $contigs = {};
137 : gdpusch 1.16 while (defined($record = <STDIN>)) {
138 :     $record =~ s/^\s+//os;
139 :     last unless $record;
140 :    
141 :     if ($record !~ m/LOCUS\s+(\S+)/os) {
142 : gdpusch 1.18 print STDERR qq(No LOCUS line for record begining with:\n)
143 :     , substr($record, 0, 160), qq(\n\n);
144 : gdpusch 1.9 }
145 : gdpusch 1.16 else {
146 : gdpusch 1.18 my $id = $1;
147 : redwards 1.24 if ($record =~ m/\nORIGIN[^\n]*\n(.*?)(\/\/|LOCUS)/os) {
148 : gdpusch 1.18 my $seq = $1;
149 : gdpusch 1.16 $seq =~ s/\s//ogs;
150 :     $seq =~ s/\d//og;
151 : efrank 1.1 undef $contigs;
152 :     $contigs->{$id} = $seq;
153 : gdpusch 1.18 &FIG::display_id_and_seq($id, \$seq, $fh_contigs);
154 : efrank 1.1 }
155 : gdpusch 1.16 else {
156 : gdpusch 1.18 warn qq(could not find contig sequence for $id\n);
157 : overbeek 1.4 next;
158 :     }
159 : gdpusch 1.16
160 : gdpusch 1.18 my ($taxon_ID, $taxonomy, $written_genome);
161 : gdpusch 1.16 if ($record =~ /\n {0,4}ORGANISM\s+(\S[^\n]+(\n\s{10,14}\S[^\n]+)*)/os) {
162 : overbeek 1.14 my $block = $1;
163 :     my @lines = split(/\n/,$block);
164 : gdpusch 1.16
165 : overbeek 1.14 my @genome = ();
166 :     my @full_tax = ();
167 : gdpusch 1.16 for ($i=0; ($i < @lines) && ($lines[$i] !~ /;/); $i++) {
168 : overbeek 1.14 push(@genome,$lines[$i]);
169 :     }
170 : gdpusch 1.16
171 :     while ($i < @lines) {
172 : overbeek 1.14 push(@full_tax,$lines[$i]);
173 :     $i++;
174 :     }
175 : gdpusch 1.16
176 : gdpusch 1.18 $genome = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @genome);
177 :     $taxonomy = join(qq( ),map { $_ =~ s/^\s*(\S.*\S).*$/$1/; $1 } @full_tax);
178 : gdpusch 1.16
179 : gdpusch 1.18 $taxonomy =~ s/\n\s+//og;
180 :     $taxonomy =~ s/ {2,}/ /og;
181 :     $taxonomy =~ s/\.$//o;
182 :     $taxonomy = $taxonomy . qq(; $genome);
183 : overbeek 1.15
184 : gdpusch 1.16 if (! $written_genome) {
185 : gdpusch 1.20 if ($force_bioname) { $genome = $force_bioname; }
186 :     print GENOME qq($genome\n);
187 : overbeek 1.4 close(GENOME);
188 : gdpusch 1.9
189 : gdpusch 1.18 if ($force_taxonomy) { $taxonomy = $force_taxonomy; }
190 :     print TAXONOMY qq($taxonomy\n);
191 :     close(TAXONOMY);
192 :    
193 : gdpusch 1.21 $written_genome = qq($taxonomy\t$genome);
194 : overbeek 1.4 }
195 : gdpusch 1.18 elsif (($written_genome ne qq($taxonomy\t$bioname)) && (!$force_bioname || !$force_taxonomy)) {
196 :     print STDERR qq(serious mismatch in GENOME/TAX for $id\n$written_genome\n$taxonomy\t$bioname\n\n);
197 : overbeek 1.4 }
198 : efrank 1.1 }
199 : gdpusch 1.16
200 :     while ($record =~ m/\n\s{4,6}CDS\s+([^\n]+(\n {20,}[^\n]*)+)/ogs) {
201 : gdpusch 1.18 my $cds = $1;
202 : gdpusch 1.16 if (($cds !~ m/\/pseudo/o) &&
203 :     (($cds !~ m/\/exception/o) || ($cds =~ m/\/translation/o))
204 :     ) {
205 : gdpusch 1.18 &process_cds($id, \$cds, $prefixP, \$idNp, $contigs, $fh_peg_tbl, $fh_peg_fasta, $fh_assigned_funcs, $fh_annotations, $fh_ec_nums);
206 : overbeek 1.7 }
207 : efrank 1.2 }
208 :    
209 : gdpusch 1.16 while ($record =~ m/\n\s{3,6}[tr]RNA\s+([^\n]+(\n\s{20,22}\S[^\n]*)+)/ogs)
210 : efrank 1.2 {
211 :     $rna = $1;
212 : gdpusch 1.18 &process_rna($id, \$rna, $prefixR, \$idNr, $contigs, $fh_rna_tbl, $fh_rna_fasta);
213 : efrank 1.1 }
214 :     }
215 :     }
216 : gdpusch 1.18 close($fh_contigs);
217 :     close($fh_peg_tbl);
218 :     close($fh_peg_fasta);
219 :     close($fh_rna_tbl);
220 :     close($fh_rna_fasta);
221 :     close($fh_assigned_funcs);
222 :     close($fh_annotations);
223 :     close($fh_ec_nums);
224 :    
225 :     if (!-s $assigned_funcs_file) {
226 : gdpusch 1.20 warn qq(WARNING: no assigned_functions in $org_dir --- deleting\n);
227 : gdpusch 1.18 unlink("$assigned_funcs_file") unless $ENV{DEBUG};
228 :     }
229 :    
230 : gdpusch 1.22 if (!-s $annotations_file) {
231 :     warn qq(WARNING: no annotations in $org_dir --- deleting\n);
232 :     unlink("$annotations_file") unless $ENV{DEBUG};
233 :     }
234 :    
235 : gdpusch 1.18 if (!-s $ec_nums_file) {
236 : gdpusch 1.20 warn qq(WARNING: no EC_numbers in $org_dir --- deleting\n);
237 : gdpusch 1.18 unlink("$ec_nums_file") unless $ENV{DEBUG};
238 :     }
239 : efrank 1.2
240 : gdpusch 1.20 if ((!-s $rna_tbl_file) || (!-s $rna_fasta_file)) {
241 :     warn qq(WARNING: no RNAs in $org_dir --- deleting\n);
242 :     system(qq(rm -rf $org_dir/Features/rna)) unless $ENV{DEBUG};;
243 :     }
244 :    
245 : gdpusch 1.18 if ((!-s $peg_tbl_file) || (!-s $peg_fasta_file)) {
246 : gdpusch 1.20 warn qq(WARNING: no PEGs in $org_dir --- deleting\n);
247 : gdpusch 1.18 system(qq(rm -rf $org_dir/Features/peg)) unless $ENV{DEBUG};
248 :     }
249 :    
250 : gdpusch 1.20 if ((!-d qq($org_dir/Features/peg)) && (!-d qq($org_dir/Features/peg))) {
251 :     warn qq(WARNING: no Features in $org_dir --- deleting\n);
252 :     system(qq(rmdir $org_dir/Features)) unless $ENV{DEBUG};;
253 : gdpusch 1.18 }
254 : gdpusch 1.20
255 : gdpusch 1.18 if (!-s $contigs_file) {
256 :     $trouble = 1;
257 :     unlink($contigs_file);
258 : gdpusch 1.20 print STDERR qq(ERROR: No contigs in $org_dir\n);
259 : gdpusch 1.18 }
260 :    
261 :     if ($trouble) {
262 :     warn qq(Genome directory $org_dir is corrupt --- deleting\n\n);
263 :     system(qq(rm -rf $org_dir)) unless $ENV{DEBUG};
264 :     }
265 :    
266 :     exit($trouble);
267 : efrank 1.1
268 :     sub process_cds {
269 : gdpusch 1.20 my ($contigID, $cdsP, $prefix, $idNp, $contigs, $fh_tbl, $fh_fasta, $fh_assign, $fh_annot, $fh_ec_nums) = @_;
270 : gdpusch 1.18 my ($id, $prot);
271 : gdpusch 1.9
272 : overbeek 1.13 ++$recnum;
273 : gdpusch 1.18 my ($loc, $precise) = &get_loc($contigID,$cdsP);
274 :     my @aliases = &get_aliases($cdsP);
275 :     my @ec_nums = &get_ec_numbers($cdsP);
276 :     my ($func, $notes) = &get_func($cdsP);
277 :     my $trans = &get_trans($cdsP);
278 :    
279 : gdpusch 1.16 if ($loc || $trans) {
280 : gdpusch 1.18 my $dna = &FIG::extract_seq($contigs,$loc);
281 :     if ($dna) {
282 : gdpusch 1.16 if ($precise) {
283 : gdpusch 1.18 if (not $trans) {
284 :     warn qq(WARNING: Translation missing for CDS $recnum; attempting to generate translation\n);
285 :    
286 :     my $genetic_code = 11;
287 :     if ($$cdsP =~ m/\/transl_table=(\d+)/o) {
288 :     $genetic_code = $1;
289 :     }
290 :    
291 :     if (not defined($trans_table{$genetic_code})) {
292 :     $trans_table{$genetic_code} = &FIG::genetic_code($genetic_code);
293 :     }
294 :    
295 :     my $prot = &FIG::translate($dna, $trans_table{$genetic_code}, 1);
296 :    
297 :     if ($prot !~ m/\*/o) {
298 :     $trans = $prot;
299 :     }
300 :     else {
301 :     warn qq(Translation contains STOPs, changing to \'x\'s, for CDS $recnum:\n$$cdsP\n);
302 :     }
303 : gdpusch 1.9 }
304 : efrank 1.1 }
305 :     }
306 : gdpusch 1.16 else {
307 : gdpusch 1.18 warn qq(Could not get DNA sequence for CDS at $loc for entry:\n$$cdsP\nof record begining with:\n)
308 :     , substr($record, 0, 160), qq(\n\n);
309 : gdpusch 1.9 }
310 :    
311 : gdpusch 1.18 $id = $prefix . qq($$idNp);
312 :     ++$$idNp;
313 : olson 1.19 print $fh_tbl qq($id\t$loc\t) . join("\t",@aliases) . "\n";
314 : gdpusch 1.16
315 :     if ($func) {
316 : gdpusch 1.20 print $fh_assign qq($id\t$func\n);
317 : efrank 1.1 }
318 : gdpusch 1.16
319 : gdpusch 1.18 if (@ec_nums) {
320 :     print $fh_ec_nums (join(qq(\t), ($id, @ec_nums)), qq(\n));
321 :     }
322 :    
323 :     foreach my $note (@$notes) {
324 :     &make_annotation($id, $note, $fh_annot);
325 :     }
326 :    
327 : gdpusch 1.16 if ($trans) {
328 : efrank 1.1 &FIG::display_id_and_seq($id,\$trans,$fh_fasta);
329 :     }
330 : gdpusch 1.16 else {
331 : gdpusch 1.18 die qq(No translation for $id);
332 : gdpusch 1.9 }
333 : efrank 1.1 }
334 : gdpusch 1.16 else {
335 : gdpusch 1.18 print &Dumper($cdsP); die qq(aborted);
336 : efrank 1.1 }
337 :     }
338 :    
339 :     sub get_loc {
340 :     my($contigID,$cdsP) = @_;
341 :     my($beg,$end,$loc,$locS,$iter,$n,$args,@pieces);
342 :     my($precise);
343 : gdpusch 1.16
344 :     if ($$cdsP =~ m/^(([^\n]+)((\n\s{21,23}[^\/ ][^\n]+)+)?)/os) {
345 : efrank 1.1 $locS = $1;
346 :     $locS =~ s/\s//g;
347 : gdpusch 1.16 $precise = ($locS !~ m/[<>]/o);
348 : gdpusch 1.10
349 : efrank 1.1 @pieces = ();
350 :     $n = 0;
351 :     $iter = 500;
352 : gdpusch 1.16 while (($locS !~ m/^\[\d+\]$/o) && $iter)
353 : efrank 1.1 {
354 :     $iter--;
355 : gdpusch 1.16 if ($locS =~ s/[<>]?(\d+)\.\.[<>]?(\d+)/\[$n\]/o) {
356 : efrank 1.1 push(@pieces,["loc",$1,$2]);
357 :     $n++;
358 :     }
359 : gdpusch 1.10
360 : gdpusch 1.16 if ($locS =~ s/([,\(])(\d+)([,\)])/$1\[$n\]$3/o) {
361 : efrank 1.1 push(@pieces,["loc",$2,$2]);
362 :     $n++;
363 :     }
364 : gdpusch 1.10
365 : gdpusch 1.16 if ($locS =~ s/join\((\[\d+\](,\[\d+\])*)\)/\[$n\]/o) {
366 : efrank 1.1 $args = $1;
367 : gdpusch 1.16 push(@pieces,["join",map { $_ =~ /^\[(\d+)\]$/; $1 } split(m/,/,$args)]);
368 : efrank 1.1 $n++;
369 :     }
370 : gdpusch 1.10
371 : gdpusch 1.16 if ($locS =~ s/complement\((\[\d+\](,\[\d+\])*)\)/\[$n\]/og) {
372 : efrank 1.1 $args = $1;
373 : gdpusch 1.16 push(@pieces,["complement", map { $_ =~ m/^\[(\d+)\]$/o; $1 } split(m/,/o, $args)]);
374 : efrank 1.1 $n++;
375 :     }
376 :     }
377 : gdpusch 1.10
378 : gdpusch 1.16 if ($locS =~ m/^\[(\d+)\]$/o) {
379 : efrank 1.1 $loc = &conv(\@pieces,$1,$contigID);
380 :     }
381 : gdpusch 1.16 else {
382 : efrank 1.1 print STDERR &Dumper(["could not parse $locS $iter",$cdsP]);
383 : gdpusch 1.18 die qq(aborted);
384 : efrank 1.1 }
385 : gdpusch 1.10
386 : gdpusch 1.16 my @locs = split(m/,/o, $loc);
387 : gdpusch 1.10 #...STOP is now included, so don't trim it off...
388 :     # &trim_stop(\@locs);
389 : efrank 1.1 $loc = join(",",@locs);
390 :     }
391 : gdpusch 1.16 else {
392 : gdpusch 1.18 print STDERR &Dumper($cdsP); die qq(could not parse location);
393 :     die qq(aborted);
394 : efrank 1.1 }
395 :     return ($loc,$precise);
396 :     }
397 :    
398 :     sub trim_stop {
399 :     my($locs) = @_;
400 :     my($to_trim,$n);
401 : gdpusch 1.16
402 : efrank 1.1 $to_trim = 3;
403 : gdpusch 1.16 while ($to_trim && (@$locs > 0)) {
404 : efrank 1.1 $n = @$locs - 1;
405 : gdpusch 1.16 if ($locs->[$n] =~ m/^(\S+)_(\d+)_(\d+)$/o) {
406 :     if ($2 <= ($3-$to_trim)) {
407 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3-$to_trim));
408 :     $to_trim = 0;
409 :     }
410 : gdpusch 1.16 elsif ($2 >= ($3 + $to_trim)) {
411 : efrank 1.1 $locs->[$n] = join("_",($1,$2,$3+$to_trim));
412 :     $to_trim = 0;
413 :     }
414 : gdpusch 1.16 else {
415 : efrank 1.1 $to_trim -= abs($3-$2) + 1;
416 :     pop @$locs;
417 :     }
418 :     }
419 : gdpusch 1.16 else {
420 : gdpusch 1.18 die qq(could not parse $locs->[$n]);
421 : efrank 1.1 }
422 :     }
423 :     }
424 :    
425 :    
426 :     sub conv {
427 :     my($pieces,$n,$contigID) = @_;
428 : gdpusch 1.16
429 : gdpusch 1.18 if ($pieces->[$n]->[0] eq qq(loc)) {
430 : efrank 1.1 return join("_",$contigID,@{$pieces->[$n]}[1..2]);
431 :     }
432 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(join)) {
433 : efrank 1.1 return join(",",map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}]);
434 :     }
435 : gdpusch 1.18 elsif ($pieces->[$n]->[0] eq qq(complement)) {
436 : gdpusch 1.16 return join(",",&complement(join(",", map { &conv($pieces,$_,$contigID) } @{$pieces->[$n]}[1..$#{$pieces->[$n]}])));;
437 : efrank 1.1 }
438 :     }
439 :    
440 :     sub complement {
441 :     my($loc) = @_;
442 :     my(@locs);
443 :    
444 :     @locs = reverse split(/,/,$loc);
445 : gdpusch 1.16 foreach $loc (@locs) {
446 :     if ($loc =~ m/^(\S+)_(\d+)_(\d+)$/o) {
447 :     $loc = join("_",($1,$3,$2));
448 :     }
449 :     else {
450 : gdpusch 1.18 confess qq(Bad location: $loc);
451 : gdpusch 1.16 }
452 : efrank 1.1 }
453 :     return join(",",@locs);
454 :     }
455 :    
456 :     sub get_aliases {
457 :     my($cdsP) = @_;
458 : gdpusch 1.18 my($id, $prefix, $alias, $db_ref);
459 : gdpusch 1.16
460 : efrank 1.1 my @aliases = ();
461 : gdpusch 1.16 while ($$cdsP =~ m/\/(protein_id|gene|locus_tag)=\"([^\"]+)\"/ogs) {
462 : gdpusch 1.20 ($type, $alias) = ($1, $2);
463 : gdpusch 1.18
464 :     # define prefixes for different types of ids
465 :     if ($type eq qq(locus_tag)){
466 :     $id = qq(locus|$alias);
467 :     }
468 :     elsif ( $type eq qq(protein_id) ) {
469 :     $id = qq(protein_id|$alias);
470 :     }
471 :     elsif ( $type eq qq(gene) ){
472 :     $id = qq(gene_name|$alias);
473 :     }
474 :     else{
475 :     $id = $alias;
476 :     }
477 :    
478 :     push(@aliases,$id);
479 : efrank 1.1 }
480 :    
481 : gdpusch 1.16 while ($$cdsP =~ m/\/db_xref=\"([^\"]+)\"/ogs) {
482 : efrank 1.1 $db_ref = $1;
483 : gdpusch 1.16 $db_ref =~ s/^GI:/gi\|/o;
484 : wilke 1.17 $db_ref =~ s/^GeneID:/geneID\|/o;
485 : gdpusch 1.16 $db_ref =~ s/SWISS-PROT:/sp\|/o;
486 : efrank 1.1 push(@aliases,$db_ref);
487 :     }
488 : gdpusch 1.16
489 : efrank 1.1 return @aliases;
490 :     }
491 :    
492 : gdpusch 1.18 sub get_ec_numbers {
493 :     my ($cdsP) = @_;
494 :     my @ec_numbers = ($$cdsP =~ m{/EC_number=\"([0-9]+\.[0-9\-]+\.[0-9\-]+\.[0-9\-]+)\"}ogs);
495 :     return @ec_numbers;
496 :     }
497 :    
498 : efrank 1.1 sub get_trans {
499 :     my($cdsP) = @_;
500 : gdpusch 1.18 my $tran = qq();
501 : efrank 1.1
502 : gdpusch 1.16 if ($$cdsP =~ m/\/translation=\"([^"]*)\"/os) {
503 : efrank 1.1 $tran = $1;
504 :     $tran =~ s/\s//gs;
505 :     }
506 : gdpusch 1.16 elsif ($$cdsP =~ m/\/protein_id=\"([^"]+)\"/o) {
507 : overbeek 1.8 $tran = $1;
508 : gdpusch 1.16 $tran =~ s/\s//ogs;
509 : overbeek 1.8 }
510 : efrank 1.1 return $tran;
511 :     }
512 :    
513 :     sub get_func {
514 :     my($cdsP) = @_;
515 : overbeek 1.11
516 : gdpusch 1.18 my $functions = [];
517 :     my $products = [];
518 :     my $prot_descs = [];
519 :     my $notes = [];
520 :    
521 :     print STDERR qq(\nRecord $recnum:\n$$cdsP\n) if $ENV{VERBOSE};
522 :    
523 : gdpusch 1.20 @$functions = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/function=\"([^\"]*)\"/ogs);
524 : gdpusch 1.18 if ($ENV{VERBOSE} && @$functions) {
525 :     print STDERR (qq(Functions: ), ((@$functions > 1) ? qq(\n) : qq()));
526 :     print STDERR (join(qq(\n), @$functions), qq(\n));
527 :     }
528 :    
529 : gdpusch 1.20 @$products = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/product=\"([^\"]*)\"/ogs);
530 : gdpusch 1.18 if ($ENV{VERBOSE} && @$products) {
531 :     print STDERR (qq(Products: ), ((@$products > 1) ? qq(\n) : qq()));
532 :     print STDERR (join(qq(\n), @$products), qq(\n));
533 : overbeek 1.4 }
534 : gdpusch 1.18
535 : gdpusch 1.20 @$prot_descs = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/prot_desc=\"([^\"]*)\"/ogs);
536 : gdpusch 1.18 if ($ENV{VERBOSE} && @$prot_descs) {
537 :     print STDERR (qq(Prot_Descs: ), ((@$prot_descs > 1) ? qq(\n) : qq()));
538 :     print STDERR (join(qq(\n), @$prot_descs), qq(\n));
539 :     }
540 :    
541 :     @$notes = map { s/[\s\n]+/ /ogs; $_ } ($$cdsP =~ m/\/note=\"([^\"]*)\"/ogs);
542 :     if ($ENV{VERBOSE} && @$notes) {
543 :     print STDERR (qq(Notes: ), ((@$notes > 1) ? qq(\n) : qq()));
544 :     print STDERR (join(qq(\n), @$notes), qq(\n));
545 :     }
546 :    
547 :     @$products = grep { !m/^hypothetical\s+protein$/io } @$products;
548 :    
549 :     my $func = qq();
550 :     my $annotations = [];
551 :     if (@$products) {
552 :     $func = join(q( / ), @$products);
553 :     @$annotations = (@$functions, @$prot_descs, @$notes);
554 :     }
555 :     elsif (@$functions) {
556 :     $func =join(q( / ), @$functions);
557 :     @$annotations =(@$prot_descs, @$notes);
558 :     }
559 :     elsif (@$prot_descs) {
560 :     $func =join(q( / ), @$prot_descs);
561 :     @$annotations = @$notes;
562 : overbeek 1.4 }
563 : gdpusch 1.16 else {
564 : gdpusch 1.18 $func = qq(hypothetical protein);
565 :     @$annotations = @$notes;
566 : efrank 1.2 }
567 : gdpusch 1.18
568 :     return ($func, $annotations);
569 : overbeek 1.4 }
570 :    
571 :     sub fixup_func {
572 :     my($func) = @_;
573 : gdpusch 1.16
574 :     $func =~ s/^COG\d+:\s+//oi;
575 :     $func =~ s/^ORFID:\S+\s+//oi;
576 : efrank 1.2 return $func;
577 :     }
578 :    
579 : overbeek 1.4 sub ok_func {
580 :     my($func) = @_;
581 : gdpusch 1.16
582 : overbeek 1.4 return (
583 : gdpusch 1.16 ($func !~ m/^[a-zA-Z]{1,3}\d+[a-zA-Z]{0,3}(\.\d+([a-zA-Z])?)?$/o) &&
584 :     ($func !~ m/^\d+$/o) &&
585 :     ($func !~ m/ensangp/oi) &&
586 :     ($func !~ m/agr_/oi) &&
587 :     ($func !~ m/^SC.*:/o) &&
588 :     ($func !~ m/^RIKEN/o) &&
589 :     ($func !~ m/\slen: \d+/o) &&
590 :     ($func !~ m/^similar to/oi) &&
591 :     ($func !~ m/^CPX/oi) &&
592 :     ($func !~ m/^\S{3,4}( or \S+ protein)?$/oi) &&
593 :     ($func !~ m/^putative$/o) &&
594 :     ($func !~ m/\scomes from/oi) &&
595 :     ($func !~ m/^unknown( function)?$/oi) &&
596 :     ($func !~ m/^hypothetical( (function|protein))?$/oi) &&
597 :     ($func !~ m/^orf /oi) &&
598 :     ($func !~ m/^ORF\s?\d+[lr]?$/oi) &&
599 :     ($func !~ m/^[a-z]{1,3}\s?\d+$/oi)
600 :     );
601 : overbeek 1.4 }
602 :    
603 : efrank 1.2 sub process_rna {
604 : overbeek 1.4 my($contigID,$cdsP,$prefix,$idNr,$contigs,$fh_tbl,$fh_dna) = @_;
605 : gdpusch 1.9 my($loc,@aliases,$func,$trans,$id,$precise,$DNA);
606 :    
607 : gdpusch 1.16 ($loc, $precise) = &get_loc($contigID,$cdsP);
608 :     $func = &get_rna_func($cdsP);
609 : efrank 1.2
610 : gdpusch 1.16 if ($loc) {
611 :     if ($dna = &FIG::extract_seq($contigs,$loc)) {
612 : gdpusch 1.18 $id = $prefix . qq($$idNr);
613 : gdpusch 1.9 $$idNr++;
614 : gdpusch 1.16
615 :     if (! $func ) {
616 : gdpusch 1.18 warn qq(WARNING: could not get func $$cdsP\n) if $ENV{VERBOSE};
617 : gdpusch 1.9 }
618 : gdpusch 1.16 else {
619 : gdpusch 1.18 print $fh_tbl qq($id\t$loc\t$func\n);
620 : gdpusch 1.9 &FIG::display_id_and_seq($id,\$dna,$fh_dna);
621 :     }
622 : efrank 1.2 }
623 : gdpusch 1.16 else {
624 : gdpusch 1.18 warn qq(Could not get DNA sequence for RNA at $loc for entry\n$$cdsP\nof record begining with:\n)
625 :     , substr($record, 0, 160), qq(\n\n);
626 : efrank 1.2 }
627 :     }
628 : gdpusch 1.16 else {
629 : gdpusch 1.18 print &Dumper($cdsP); die qq(aborted);
630 : efrank 1.2 }
631 :     }
632 :    
633 :     sub get_rna_func {
634 :     my($cdsP) = @_;
635 : gdpusch 1.18 my $func = qq();
636 : gdpusch 1.16
637 :     if ($$cdsP =~ m/\/product=\"([^"]*)\"/os) {
638 : efrank 1.2 $func = $1;
639 :     $func =~ s/\s+/ /gs;
640 :     }
641 : gdpusch 1.16 elsif ($$cdsP =~ m/\/gene=\"([^"]*)\"/os) {
642 : efrank 1.2 $func = $1;
643 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
644 : efrank 1.2 }
645 : gdpusch 1.16 elsif ($$cdsP =~ m/\/note=\"([^"]*)\"/os) {
646 : efrank 1.1 $func = $1;
647 : gdpusch 1.16 $func =~ s/\s+/ /ogs;
648 : efrank 1.1 }
649 :     return $func;
650 :     }
651 : gdpusch 1.18
652 :     sub make_annotation {
653 :     my ($fid, $note, $fh_ann) = @_;
654 :     # print STDERR Dumper($fid, $note, $fh_ann);
655 :    
656 :     print $fh_ann ($fid, qq(\n));
657 :     print $fh_ann (time, qq(\n));
658 :     print $fh_ann qq(parse_genbank\n);
659 :     print $fh_ann ($note, qq(\n));
660 :     print $fh_ann qq(//\n);
661 :    
662 :     return 1;
663 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3