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

Annotation of /FigKernelScripts/parse_genbank.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3