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

Annotation of /FigKernelPackages/gjogenbank.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : golsen 1.3 #
2 : gdpusch 1.1 # This is a SAS component.
3 : golsen 1.3 #
4 : gdpusch 1.1 package gjogenbank;
5 :    
6 :     #===============================================================================
7 : golsen 1.3 # Parse one or more GenBank entries in a file into perl structures.
8 :     # All of the entries in the file as a list:
9 : gdpusch 1.1 #
10 : golsen 1.3 # @entries = parse_genbank( ) # \*STDIN
11 :     # \@entries = parse_genbank( ) # \*STDIN
12 :     # @entries = parse_genbank( \*FH )
13 :     # \@entries = parse_genbank( \*FH )
14 :     # @entries = parse_genbank( $file )
15 :     # \@entries = parse_genbank( $file )
16 :     #
17 :     # One entry per call:
18 :     #
19 :     # $entry = parse_next_genbank( ) # STDIN
20 :     # $entry = parse_next_genbank( \*FH )
21 :     # $entry = parse_next_genbank( $file )
22 :     #
23 :     # Error or end-of-file returns undef.
24 :     #
25 :     # Each entry is a hash with key value pairs, some of which have special
26 :     # processing.
27 :     #
28 :     # ACCESSION => [ Accession, ... ]
29 : golsen 1.7 # COMMENT => [ Comment_line, ... ] # Original wrap is preserved
30 :     # DBLINK => [ Field, ... ]
31 :     # DBSOURCE => [ Field, ... ]
32 : golsen 1.3 # DEFINITION => Definition
33 : golsen 1.7 # KEYWORDS => [ Keyword, ... ]
34 :     # LOCUS => Locus_id
35 : golsen 1.3 # ORGANISM => Organism_name
36 : golsen 1.7 # ORIGIN => Description_of_sequence_origin
37 : golsen 1.3 # REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]
38 :     # SEQUENCE => Sequence
39 :     # SOURCE => Source_string
40 :     # TAXONOMY => [ Taxon, Taxon, ... ]
41 :     # VERSION => [ Version, Other_information ]
42 :     #
43 : golsen 1.7 # Data that are more derived or parsed:
44 :     #
45 :     # date => Date
46 :     # division => Division
47 :     # ftr_list => [ [ type, loc, quals ], ... ]
48 :     # geometry => linear | circular
49 :     # gi => Entry_gi_number
50 :     # is_protein => Bool
51 :     # key_index => { Keyword => 1, ... }
52 :     # locus_id => Locus_id
53 :     # mol_type => Mol_type
54 :     #
55 :     # Feature records are merged by type. Slash is removed from qualifier name.
56 : golsen 1.3 # Surrounding quotation marks are stripped from qualifier values.
57 :     #
58 :     # FEATURES => { Type => [ [ Location, { Qualifier => \@values } ],
59 :     # [ Location, { Qualifier => \@values } ],
60 :     # ...
61 :     # ],
62 :     # Type => ...
63 :     # }
64 : gdpusch 1.1 #
65 :     #
66 : golsen 1.3 # Access functions to parts of structure:
67 : gdpusch 1.1 #
68 : golsen 1.4 # @types = feature_types( $entry );
69 :     # \@types = feature_types( $entry );
70 :     #
71 :     # @ftrs = features_of_type( $entry, @types );
72 :     # \@ftrs = features_of_type( $entry, @types );
73 :     # @ftrs = features_of_type( $entry, \@types );
74 :     # \@ftrs = features_of_type( $entry, \@types );
75 :     #
76 : gdpusch 1.1 # Sequence of a feature, optionally including information on partial ends.
77 :     #
78 :     # $seq = ftr_dna( $dna, $ftr )
79 :     # $seq = ftr_dna( \$dna, $ftr )
80 : golsen 1.3 # ( $seq, $partial_5, $partial_3 ) = ftr_dna( $dna, $ftr ) # boolean of > or < in location
81 : gdpusch 1.1 # ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna, $ftr )
82 :     #
83 :     # $ftr_location = location( $ftr ) # Returns empty string on failure.
84 :     #
85 : golsen 1.3 # Identify features with partial 5' or 3' ends.
86 : gdpusch 1.1 #
87 : golsen 1.3 # $partial_5_prime = partial_5_prime( $ftr )
88 :     # $partial_3_prime = partial_3_prime( $ftr )
89 : gdpusch 1.1 #
90 : golsen 1.3 # \%ftr_qualifiers = qualifiers( $ftr ) # Returns empty hash reference on failure.
91 : gdpusch 1.1 #
92 : golsen 1.3 # $gene = gene( $ftr )
93 :     # @gene_and_synonyms = gene( $ftr )
94 : gdpusch 1.1 #
95 : golsen 1.3 # $id = CDS_id( $ftr ) # Prefer protein_id as id:
96 :     # $id = CDS_gi_or_id( $ftr ) # Prefer gi number as id:
97 :     # $gi = CDS_gi( $ftr ) # gi number or nothing:
98 : gdpusch 1.1 #
99 : golsen 1.3 # $product = product( $ftr )
100 :     #
101 :     # @EC_number = EC_number( $ftr )
102 :     # \@EC_number = EC_number( $ftr )
103 :     #
104 :     # $translation = CDS_translation( $ftr ) # Uses in situ if found
105 :     # $translation = CDS_translation( $ftr, $dna ) # If not in feature, translate
106 :     # $translation = CDS_translation( $ftr, \$dna )
107 :     # $translation = CDS_translation( $ftr, $entry )
108 : gdpusch 1.1 #
109 :     # Convert GenBank location to [ [ $contig, $begin, $dir, $length ], ... ]
110 :     #
111 :     # \@cbdl = genbank_loc_2_cbdl( $loc, $contig_id )
112 :     #
113 :     # Convert GenBank location to a SEED or Sapling location string.
114 :     #
115 :     # $loc = genbank_loc_2_seed( $acc, $loc )
116 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_seed( $acc, $loc )
117 :     #
118 :     # $loc = genbank_loc_2_sapling( $acc, $loc )
119 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_sapling( $acc, $loc )
120 :     #
121 :     # Convert a [ [ contig, begin, dir, length ], ... ] location to GenBank.
122 :     #
123 :     # $gb_location = cbdl_2_genbank( \@cbdl )
124 :     # ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )
125 :     #
126 :     #===============================================================================
127 : golsen 1.7 # Write a GenBank entry:
128 :     #
129 :     # write_genbank( \%entry );
130 :     # write_genbank( $file_name, \%entry );
131 :     # write_genbank( \*FH, \%entry );
132 :     #
133 :     #===============================================================================
134 : gdpusch 1.1
135 :     use strict;
136 :    
137 : golsen 1.7 use Data::Dumper;
138 :    
139 : gdpusch 1.1 require Exporter;
140 : golsen 1.3 our @ISA = qw( Exporter );
141 :     our @EXPORT = qw( parse_genbank
142 :     parse_next_genbank
143 : golsen 1.4 feature_types
144 :     features_of_type
145 : golsen 1.7 feature_list
146 :     ftr_dna
147 :    
148 :     genbank_loc_2_seed
149 :     genbank_loc_2_sapling
150 :     genbank_loc_2_cbdl
151 :     cbdl_2_genbank
152 :    
153 :     write_genbank
154 : golsen 1.3 );
155 : gdpusch 1.1
156 : golsen 1.7
157 :     # An approximate ordering of the common qualifiers in GenBank feature table
158 :     # entries:
159 :    
160 :     my %qualifier_order;
161 :     {
162 :     my $tmp_n = 1;
163 :     %qualifier_order = map { $_ => $tmp_n++ }
164 :     qw( organism
165 :     sub_species
166 :     sex
167 :     mating_type
168 :     chromosome
169 :     macronuclear
170 :     plasmid
171 :     organelle
172 :     strain
173 :     sub_strain
174 :     cultivar
175 :     serotype
176 :     serovar
177 :     variety
178 :     isolate
179 :     tissue_type
180 :     cell_type
181 :     dev_stage
182 :     cell_line
183 :     rearranged
184 :     clone_lib
185 :     clone
186 :     sub_clone
187 :     host
188 :     lab_host
189 :     mol_type
190 :     direction
191 :    
192 :     rpt_type
193 :     rpt_family
194 :    
195 :     operon
196 :     gene
197 :     gene_synonym
198 :     allele
199 :     locus_tag
200 :     codon_start
201 :     transl_table
202 :    
203 :     anticodon
204 :     product
205 :     function
206 :     EC_number
207 :     protein_id
208 :     locus_peptide
209 :     pseudo
210 :    
211 :     note
212 :     db_xref
213 :    
214 :     translation
215 :     );
216 :     }
217 :    
218 : golsen 1.5
219 : golsen 1.8 # Qualifiers that do not require values:
220 :    
221 :     my %valueless_qual = map { $_ => 1 }
222 :     qw( artificial_location
223 :     environmental_sample
224 :     focus
225 :     germline
226 :     macronuclear
227 :     partial
228 :     proviral
229 :     pseudo
230 :     rearranged
231 :     ribosomal_slippage
232 :     transgenic
233 :     trans_splicing
234 :     );
235 :    
236 :    
237 : gdpusch 1.1 #===============================================================================
238 :     #
239 :     # @entries = parse_genbank( ) # \*STDIN
240 :     # \@entries = parse_genbank( ) # \*STDIN
241 :     # @entries = parse_genbank( \*FH )
242 :     # \@entries = parse_genbank( \*FH )
243 :     # @entries = parse_genbank( $file )
244 :     # \@entries = parse_genbank( $file )
245 :     #
246 : golsen 1.3 #===============================================================================
247 :     my %genbank_streams;
248 :    
249 :     sub parse_genbank
250 :     {
251 :     my $file = shift;
252 :    
253 :     my ( $fh, $close ) = &input_filehandle( $file );
254 :     $fh or return wantarray ? () : [];
255 :    
256 :     my @entries;
257 :     while ( my $entry = parse_one_genbank_entry( $fh ) ) { push @entries, $entry }
258 :    
259 :     close $fh if $close;
260 :    
261 :     wantarray ? @entries : \@entries;
262 :     }
263 :    
264 :    
265 :     #-------------------------------------------------------------------------------
266 :     # Read and parse a GenBank file, on entry at a time. Successive calls with
267 :     # same parameter will return successive entries. Calls to different files
268 :     # can be interlaced.
269 : gdpusch 1.1 #
270 : golsen 1.3 # $entry = parse_next_genbank( ) # STDIN
271 :     # $entry = parse_next_genbank( \*FH )
272 :     # $entry = parse_next_genbank( $file )
273 : gdpusch 1.1 #
274 : golsen 1.3 # Error or end-of-file returns undef.
275 :     #-------------------------------------------------------------------------------
276 :     sub parse_next_genbank
277 :     {
278 :     my $file = shift;
279 :    
280 :     my $stream = $genbank_streams{ $file || '' };
281 :     if ( ! $stream )
282 :     {
283 :     $stream = [ &input_filehandle( $file ) ]; # Value is ( $fh, $close )
284 :     $stream->[0] or return undef; # Got a file handle?
285 :     $genbank_streams{ $file || '' } = $stream;
286 :     }
287 :    
288 :     my ( $fh, $close ) = @$stream;
289 :     my $entry = parse_one_genbank_entry( $fh );
290 :    
291 :     if ( ! $entry ) { close $fh if $close; delete $genbank_streams{ $file || '' }; }
292 :    
293 :     return $entry;
294 :     }
295 :    
296 :     #-------------------------------------------------------------------------------
297 :     # If it should be necessary to close a stream openned by parse_next_genbank()
298 :     # before it reaches the end-of-file, this will do it.
299 : gdpusch 1.1 #
300 : golsen 1.3 # close_next_genbank( ) # does nothing
301 :     # close_next_genbank( \*FH )
302 :     # close_next_genbank( $file )
303 : gdpusch 1.1 #
304 : golsen 1.3 #-------------------------------------------------------------------------------
305 :     sub close_next_genbank
306 : gdpusch 1.1 {
307 :     my $file = shift;
308 : golsen 1.3 my $stream = $genbank_streams{ $file || '' };
309 :     close $stream->[0] if $stream && ref $stream eq 'ARRAY' && $stream->[1];
310 :     }
311 :    
312 : gdpusch 1.1
313 : golsen 1.3 #-------------------------------------------------------------------------------
314 :     # Parse the next GenBank format entry read from an open file handle. This is
315 :     # primarily intended as an internal function called through parse_genbank()
316 :     # or parse_next_genbank().
317 :     #
318 :     # \%entry = parse_one_genbank_entry( \*FH )
319 :     #
320 :     # Error or end-of-file returns undef
321 :     #-------------------------------------------------------------------------------
322 :     sub parse_one_genbank_entry
323 :     {
324 :     my $fh = shift;
325 : golsen 1.7 local $_;
326 : gdpusch 1.1
327 :     my $state = 0;
328 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
329 : golsen 1.3
330 : gdpusch 1.1 # 0 = Looking for LOCUS
331 :     # 1 = Header information
332 :     # 2 = Features
333 :     # 3 = Sequence
334 :     # -1 = Error
335 :    
336 : golsen 1.3 my %entry = ();
337 : gdpusch 1.1
338 :     # 1 2 3 4 5 6 7 8
339 :     # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
340 :     # LOCUS NC_000909 1664970 bp DNA circular BCT 03-DEC-2005
341 :     # LOCUS DGRINCAD_6 9696 BP DS-DNA SYN 22-AUG-2006
342 :     #
343 : golsen 1.3 while ( $state == 0 )
344 :     {
345 : golsen 1.5 if ( s/^LOCUS\s+// )
346 : golsen 1.3 {
347 : golsen 1.5 my @parts = split;
348 : golsen 1.7 $entry{ locus_id } = $entry{ LOCUS } = shift @parts;
349 :     $entry{ date } = pop @parts if $parts[-1] =~ m/^\d+-[A-Z][A-Z][A-Z]-\d+$/i;
350 :     $entry{ division } = pop @parts if $parts[-1] =~ m/^\S\S\S$/;
351 :     $entry{ geometry } = pop @parts if $parts[-1] =~ m/^(lin|circ)/i;
352 :     $entry{ mol_type } = pop @parts if $parts[-1] =~ m/na$/i;
353 : golsen 1.3 $state = 1;
354 :     }
355 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
356 :     }
357 :    
358 :     # Reading the header requires merging continuations, then dealing
359 :     # with the data:
360 :    
361 :     while ( $state == 1 )
362 :     {
363 :     if ( /^FEATURES / )
364 : gdpusch 1.1 {
365 : golsen 1.3 $state = 2;
366 : gdpusch 1.1 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
367 :     }
368 :    
369 : golsen 1.3 elsif ( /^ORIGIN / )
370 :     {
371 :     $state = 3;
372 : golsen 1.7 last;
373 : golsen 1.3 }
374 : gdpusch 1.1
375 : golsen 1.3 elsif ( /^REFERENCE / )
376 : gdpusch 1.1 {
377 : golsen 1.3 my $ref;
378 :     ( $ref, $state, $_ ) = read_ref( $_, $fh );
379 :     push @{ $entry{ REFERENCES } }, $ref if $ref;
380 :     defined() or $state = -1;
381 :     }
382 :    
383 : golsen 1.7 elsif ( /^(.{10}) +(.*\S)\s*$/ ) # Any other keyword
384 : golsen 1.3 {
385 : golsen 1.7 my ( $tag, @value ) = ( uc( $1 ), $2 );
386 :     $tag =~ s/^ +| +$//g;
387 : golsen 1.3
388 :     # Merge continuations:
389 :    
390 : golsen 1.7 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
391 : golsen 1.3
392 : golsen 1.7 while ( $state >= 0 && s/^ {12}\s*// )
393 : gdpusch 1.1 {
394 : golsen 1.7 s/\s+$//;
395 :     push @value, $_;
396 : gdpusch 1.1 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
397 :     }
398 :    
399 : golsen 1.3 # Special case formats
400 :    
401 : golsen 1.7 if ( $tag eq 'COMMENT' )
402 : gdpusch 1.1 {
403 : golsen 1.7 push @{ $entry{ COMMENT } }, @value;
404 : gdpusch 1.1 }
405 : golsen 1.7
406 :     elsif ( $tag eq 'DBLINK' )
407 : golsen 1.3 {
408 : golsen 1.7 my $data = '';
409 :     my @dbs;
410 :     foreach ( @value )
411 :     {
412 :     if ( $data && /: / )
413 :     {
414 :     push @dbs, $data;
415 :     $data = '';
416 :     }
417 :     $data .= $_;
418 :     }
419 :     push @dbs, $data if length $data;
420 :    
421 :     $entry{ DBLINK } = \@dbs if @dbs;
422 : golsen 1.3 }
423 : golsen 1.7
424 :     elsif ( $tag eq 'DBSOURCE' )
425 : golsen 1.3 {
426 : golsen 1.7 my $data = '';
427 :     my @dbs;
428 :     foreach ( @value )
429 :     {
430 :     if ( $data && /: / )
431 :     {
432 :     push @dbs, $data;
433 :     $data = '';
434 :     }
435 :     $data .= $_;
436 :     }
437 :     push @dbs, $data if length $data;
438 :    
439 :     $entry{ DBSOURCE } = \@dbs if @dbs;
440 : golsen 1.3 }
441 : golsen 1.7
442 : golsen 1.3 elsif ( $tag eq 'ORGANISM' )
443 : gdpusch 1.1 {
444 : golsen 1.7 my $org = shift @value;
445 : golsen 1.3 $entry{ ORGANISM } = $org;
446 : golsen 1.7 my $tax = @value ? join( ' ', @value ) : '';
447 :     $tax =~ s/\s*\.$//;
448 : golsen 1.3 $entry{ TAXONOMY } = [ split /; */, $tax ] if $tax;
449 : gdpusch 1.1 }
450 : golsen 1.7
451 : golsen 1.3 else
452 : gdpusch 1.1 {
453 : golsen 1.7 my $imax = @value - 2;
454 :     if ( $imax > 0 )
455 :     {
456 :     foreach ( @value[ 0 .. $imax ] ) { $_ .= ' ' if ! /-$/ }
457 :     }
458 :     my $value = join( '', @value );
459 :    
460 :     if ( $tag eq 'ACCESSION' )
461 :     {
462 :     $entry{ ACCESSION } = [ split / +/, $value ];
463 :     }
464 :     elsif ( $tag eq 'VERSION' )
465 :     {
466 :     if ( $value =~ / +GI:(\d+)/ )
467 :     {
468 :     $entry{ gi } = $1;
469 :     $value =~ s/ +GI:\d+//;
470 :     }
471 :     $entry{ VERSION } = [ split / +/, $value ];
472 :     }
473 :     elsif ( $tag eq 'KEYWORDS' )
474 :     {
475 :     $value =~ s/\s*\.\s*$//;
476 :     if ( $value )
477 :     {
478 :     my @keys = split /; */, $value;
479 :     $entry{ KEYWORDS } = \@keys;
480 :     $entry{ key_index } = { map { $_ => 1 } @keys };
481 :     }
482 :     }
483 :    
484 :     # Generic case:
485 :    
486 :     else
487 :     {
488 :     $entry{ $tag } = $value;
489 :     }
490 : golsen 1.3 }
491 : gdpusch 1.1
492 : golsen 1.3 # To know that we are at end of continuations, we must have
493 :     # read another line.
494 : gdpusch 1.1
495 : golsen 1.3 defined() or $state = -1;
496 :     }
497 : gdpusch 1.1
498 : golsen 1.3 else # This is really a format error, but let's skip it.
499 :     {
500 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
501 :     }
502 :     }
503 : gdpusch 1.1
504 : golsen 1.8 #
505 : golsen 1.7 # Reading FEATURES requires merging continuations, then dealing
506 : golsen 1.3 # with the data:
507 : golsen 1.8 #
508 : golsen 1.7 while ( $state == 2 && ( /^ (\S+)\s+(\S+)/ ) )
509 : golsen 1.3 {
510 : golsen 1.7 my ( $type, $loc ) = ( $1, $2 );
511 :    
512 :     # Collect the rest of the location:
513 :    
514 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
515 :     while ( $state >= 0 && /^ {15}\s*([^\/\s]\S*)\s*$/ )
516 : golsen 1.3 {
517 : golsen 1.7 $loc .= $1;
518 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
519 : golsen 1.3 }
520 : gdpusch 1.1
521 : golsen 1.7 # Collect qualiiers:
522 :    
523 :     my ( $qualif, $value, %qualifs, @qualif_order );
524 :     while ( $state == 2 && ( $_ =~ /^\s*\/\w+/ ) )
525 : golsen 1.3 {
526 : golsen 1.7 # Qualifiers without = get an undef value (intentionally)
527 :    
528 :     ( $qualif, undef, $value ) = /^\s*\/(\w+)(=(.*))?/;
529 :     # Quoted strings can have value lines that start with /, so
530 :     # we must track quotation marks.
531 :    
532 :     my $nquote = $value ? $value =~ tr/"// : 0;
533 : gdpusch 1.1
534 : golsen 1.3 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
535 : golsen 1.7 while ( $state >= 0 && /^ {15}/ && ( $nquote % 2 || ( ! /^\s*\/\w/ ) ) )
536 : gdpusch 1.1 {
537 : golsen 1.3 s/^ +//;
538 : golsen 1.7 $nquote += tr/"//;
539 : golsen 1.8 $value .= ( $value =~ /\S-$/ ? '' : ' ' ) . $_;
540 : gdpusch 1.1 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
541 :     }
542 :    
543 : golsen 1.7 if ( $nquote % 2 )
544 :     {
545 :     print STDERR "Feature quotation nesting error: $type, $loc, $qualif=$value\n";
546 :     exit;
547 :     }
548 : gdpusch 1.1
549 : golsen 1.7 if ( $qualif )
550 : gdpusch 1.1 {
551 : golsen 1.8 if ( $valueless_qual{ $qualif } )
552 :     {
553 :     $value = 1;
554 :     }
555 :     elsif ( ! defined $value || ! length $value )
556 :     {
557 :     next;
558 :     }
559 :     else
560 :     {
561 :     $value =~ s/""/"/g if $value =~ s/^"(.*)"$/$1/;
562 :     $value =~ s/ +//g if $qualif eq 'translation';
563 :     }
564 : golsen 1.3
565 : golsen 1.7 push @qualif_order, $qualif if ! $qualifs{ $qualif };
566 :     push @{ $qualifs{ $qualif } }, $value;
567 :     }
568 :     }
569 :    
570 :     $qualifs{ '_qualif_order' } = \@qualif_order if @qualif_order;
571 :     push @{ $entry{ FEATURES }->{ $type } }, [ $loc, \%qualifs ] if ( $type && $loc );
572 :     push @{ $entry{ ftr_list } }, [ $type, $loc, \%qualifs ] if ( $type && $loc );
573 :    
574 :     defined() or $state = -1;
575 :     }
576 :    
577 :     $state = 3 if $state > 0;
578 : gdpusch 1.1
579 : golsen 1.7 # Only a few fields are allowed in the standard, but:
580 : gdpusch 1.1
581 : golsen 1.7 while ( $state == 3 )
582 :     {
583 :     if ( /^ORIGIN / )
584 :     {
585 :     $entry{ ORIGIN } = $1 if /^ORIGIN \s+(\S.*\S)\s*$/;
586 :     $state = 4;
587 :     }
588 : gdpusch 1.1
589 : golsen 1.8 # 2011/01/23 -- GJO
590 :     #
591 :     # NCBI has stopped including BASE COUNT. Their example is a disaster
592 :     # (http://www.ncbi.nlm.nih.gov/genbank/genomesubmit-Examples.html):
593 :     #
594 :     # BASE COUNT 1165552 a 648314 c 647106 g1169556 t
595 :     # ORIGIN
596 :     #
597 : golsen 1.7 elsif ( s/^BASE COUNT\s+// )
598 :     {
599 : golsen 1.8 # $entry{ BASE_COUNT } = { reverse split };
600 : golsen 1.7 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
601 :     }
602 : gdpusch 1.1
603 : golsen 1.7 elsif ( /^(.{10}) (.*)$/ ) # Any other keyword
604 :     {
605 :     my ( $tag, @value ) = map { s/^ +| +$//g; $_ } ( uc( $1 ), $2 );
606 : gdpusch 1.1
607 : golsen 1.7 # Merge continuations:
608 : gdpusch 1.1
609 : golsen 1.7 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
610 : gdpusch 1.1
611 : golsen 1.7 while ( $state >= 0 && s/^ {12}\s*// )
612 :     {
613 :     s/\s+$//;
614 :     push @value, $_;
615 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
616 : gdpusch 1.1 }
617 :    
618 : golsen 1.7 # Special case formats
619 : gdpusch 1.1
620 : golsen 1.7 if ( $tag eq 'COMMENT' )
621 :     {
622 :     push @{ $entry{ COMMENT } }, @value;
623 :     }
624 : gdpusch 1.1
625 : golsen 1.7 else
626 :     {
627 :     my $imax = @value - 2;
628 :     if ( $imax > 0 )
629 :     {
630 :     foreach ( @value[ 0 .. $imax ] ) { $_ .= ' ' if ! /-$/ }
631 :     }
632 :     $entry{ $tag } = join( '', @value );
633 :     }
634 : golsen 1.3 }
635 : gdpusch 1.1
636 : golsen 1.7 else # This is really a format error, but let's skip it.
637 : gdpusch 1.1 {
638 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
639 : golsen 1.3 }
640 :     }
641 : gdpusch 1.1
642 : golsen 1.3 # Read the sequence:
643 : gdpusch 1.1
644 : golsen 1.7 while ( $state == 4 )
645 : golsen 1.3 {
646 : golsen 1.7 my @sequence;
647 :    
648 : golsen 1.3 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
649 : golsen 1.7 while ( $state >= 0 && s/^ *\d+ +// )
650 : golsen 1.3 {
651 :     s/[^A-Za-z]+//g;
652 :     push @sequence, $_;
653 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
654 : gdpusch 1.1 }
655 : golsen 1.7 $entry{ SEQUENCE } = join '', @sequence;
656 : golsen 1.3
657 :     $state = ( $_ eq '//' ) ? 0 : -1 if $state >= 0;
658 : gdpusch 1.1 }
659 :    
660 : golsen 1.3 $state >= 0 ? \%entry : undef;
661 : gdpusch 1.1 }
662 :    
663 :    
664 :     #-------------------------------------------------------------------------------
665 :     # Parse a reference.
666 :     #-------------------------------------------------------------------------------
667 :     sub read_ref
668 :     {
669 :     my ( $line, $fh ) = @_;
670 :     my $state = 1;
671 : golsen 1.7 my %ref = ();
672 :    
673 :     if ( $line =~ /\s*(\d+)/ )
674 :     {
675 :     $ref{ ref_num } = $1;
676 :     $ref{ ref_note } = $1 if $line =~ /\s*\d+\s+\((.*)\)\s*$/;
677 :     }
678 :    
679 : gdpusch 1.1 if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
680 :    
681 :     my ( $tag, $value );
682 :     while ( ( $state >= 0 ) && /^ / )
683 :     {
684 :     if ( substr( $_, 0, 10 ) =~ /\S/ )
685 :     {
686 :     ( $tag, $value ) = $_ =~ /\s*(\w+)\s+(.*)$/;
687 :     }
688 :     elsif ( /\S/ )
689 :     {
690 :     s/^ +//;
691 :     $value .= " $_" if $_;
692 :     }
693 :    
694 :     if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
695 :    
696 :     if ( ( $state < 0 ) || ( /^ {0,9}\S/ ) )
697 :     {
698 :     if ( $tag && $value )
699 :     {
700 :     $ref{ $tag } = $value;
701 :     }
702 :     $tag = $value = undef;
703 :     }
704 :     }
705 :    
706 : golsen 1.3 ( ( keys %ref ? \%ref : undef ), $state, $_ )
707 : gdpusch 1.1 }
708 :    
709 :    
710 :    
711 :     #===============================================================================
712 :     # Access methods for some features and feature data
713 :     #===============================================================================
714 : golsen 1.4 #
715 :     # @types = feature_types( $entry );
716 :     # \@types = feature_types( $entry );
717 :     #
718 :     #-------------------------------------------------------------------------------
719 :     sub feature_types
720 :     {
721 :     my $entry = shift;
722 :     return wantarray ? () : [] if ! ( $entry && ref $entry eq 'HASH' );
723 :    
724 :     my $ftrs = $entry->{ FEATURES };
725 :     return wantarray ? () : [] if ! ( $ftrs && ref $ftrs eq 'HASH' );
726 :    
727 : golsen 1.7 my @types = sort { lc $a cmp lc $b } grep { ! /^_/ } keys %$ftrs;
728 : golsen 1.4 wantarray ? @types : \@types;
729 :     }
730 :    
731 :    
732 : golsen 1.7 #-------------------------------------------------------------------------------
733 :     # This is really lame for more than one type because the individual features
734 :     # are not marked with their types. This should be more effecient than the
735 :     # following method for finding all the features of one type.
736 :     #
737 :     # @ftrs = features_of_type( $entry, @types );
738 :     # \@ftrs = features_of_type( $entry, @types );
739 :     # @ftrs = features_of_type( $entry, \@types );
740 :     # \@ftrs = features_of_type( $entry, \@types );
741 :     #
742 :     #-------------------------------------------------------------------------------
743 : golsen 1.4 sub features_of_type
744 :     {
745 :     my $entry = shift;
746 :     return wantarray ? () : [] if ! ( $entry && ref $entry eq 'HASH' );
747 :    
748 :     my $ftrs = $entry->{ FEATURES };
749 :     return wantarray ? () : [] if ! ( $ftrs && ref $ftrs eq 'HASH' );
750 :    
751 :     my @types = ( ! $_[0] ) ? sort { lc $a cmp lc $b } keys %$ftrs
752 :     : ( ref $_[0] eq 'ARRAY' ) ? @{ $_[0] }
753 :     : @_;
754 :    
755 :     my @ftrs = map { @{ $ftrs->{ $_ } || [] } } @types;
756 :     wantarray ? @ftrs : \@ftrs;
757 :     }
758 :    
759 :    
760 : golsen 1.7 #-------------------------------------------------------------------------------
761 :     # Get a list of a features as triples [ $type, $loc, \%quals ] in genome
762 :     # order:
763 :     #
764 :     # @features = feature_list( $entry )
765 :     # \@features = feature_list( $entry )
766 :     # @features = feature_list( $entry, $type )
767 :     # \@features = feature_list( $entry, $type )
768 :     #
769 :     # Uses $entry->{ ftr_list } to cache the results (before filtering by type)
770 :     #-------------------------------------------------------------------------------
771 :     sub feature_list
772 :     {
773 :     my $entry = shift;
774 :    
775 :     if ( ! $entry->{ ftr_list } || ref $entry->{ ftr_list } ne 'ARRAY' )
776 :     {
777 :     my $features = $entry->{ FEATURES };
778 :    
779 :     # Is it a list of [ $type, $location, \@qualifiers ]?
780 :    
781 :     if ( ref $features eq 'ARRAY' )
782 :     {
783 :     $entry->{ ftr_list } = $features;
784 :     }
785 :    
786 :     # Is it a hash of ( $type => [ [ $location, \@qualifiers ], ... ] )?
787 :     # Build a list of features:
788 :    
789 :     elsif ( ref $features eq 'HASH' )
790 :     {
791 :     my @features;
792 :     foreach my $type ( keys %$features )
793 :     {
794 :     push @features, map { [ $type, @$_ ] } @{ $features->{ $type } };
795 :     }
796 :     $entry->{ ftr_list } = sort_feature_list( \@features );
797 :     }
798 :     }
799 :    
800 :     if ( $_[0] )
801 :     {
802 :     my @features = grep { $_->[0] eq $_[0] } @{ $entry->{ ftr_list } };
803 :     return wantarray ? @features : \@features;
804 :     }
805 :    
806 :     wantarray ? @{ $entry->{ ftr_list } } : $entry->{ ftr_list };
807 :     }
808 :    
809 :    
810 :     #-------------------------------------------------------------------------------
811 :     # Order features by location.
812 :     #
813 :     # @feautres = sort_feature_list( @features )
814 :     # \@feautres = sort_feature_list( \@features )
815 :     #
816 :     # Handles [ $location, \%quals] and [ $type, $location, \%quals ]
817 :     #-------------------------------------------------------------------------------
818 :     sub sort_feature_list
819 :     {
820 :     my $by_ref = $_[0] && ( ref $_[0] eq 'ARRAY' )
821 :     && $_[0]->[0] && ( ref $_[0]->[0] eq 'ARRAY' );
822 :    
823 :     my @features = map { $_->[0] }
824 :     sort { $a->[1] <=> $b->[1] || $b->[2] <=> $a->[2] }
825 :     map { [ $_, end_coordinates( $_->[-2] ) ] }
826 :     $by_ref ? @{ $_[0] } : @_;
827 :    
828 :     wantarray ? @features : \@features;
829 :     }
830 :    
831 :    
832 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
833 :     # For a feature location, find its left and right end coordinates.
834 :     #
835 :     # ( $left, $right ) = end_coordinates( $location )
836 :     #
837 :     # Rather than parsing, extract a list of coordinates, and take the extremes.
838 :     #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
839 :     sub end_coordinates
840 :     {
841 :     local $_ = shift;
842 :     s/([,(])[^,(]+\:[^,)]+([,)])/$1$2/g; # Delete locations in other contigs
843 :     ( sort { $a <=> $b } m/(\d+)/g )[0,-1]; # Return the lowest and highest values
844 :     }
845 :    
846 :    
847 : golsen 1.4 #-------------------------------------------------------------------------------
848 : golsen 1.3 # Sequence of a feature. In list context, include information on partial ends.
849 :     # Can get extract the data from a DNA string, reference to a string, or from
850 :     # the SEQUENCE in an entry.
851 :     #
852 :     # $seq = ftr_dna( $dna, $ftr )
853 :     # $seq = ftr_dna( \$dna, $ftr )
854 :     # $seq = ftr_dna( $entry, $ftr )
855 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( $dna, $ftr )
856 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna, $ftr )
857 :     # ( $seq, $partial_5, $partial_3 ) = ftr_dna( $entry, $ftr )
858 : gdpusch 1.1 #
859 : golsen 1.7 # Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
860 : gdpusch 1.1 #-------------------------------------------------------------------------------
861 :     sub ftr_dna
862 :     {
863 : golsen 1.3 my ( $dna, $ftr ) = @_;
864 :     return undef if ! ( $dna && $ftr );
865 :    
866 :     my $dnaR = ref $dna eq 'SCALAR' ? $dna
867 :     : ref $dna eq 'HASH' && $dna->{ SEQUENCE } ? \$dna->{ SEQUENCE }
868 :     : ! ref $dna ? \$dna
869 :     : undef;
870 :     return undef if ! $dnaR;
871 : gdpusch 1.1
872 : golsen 1.6 my $have_lib = 0;
873 :     eval { require gjoseqlib; $have_lib = 1; };
874 :     return undef if ! $have_lib;
875 :    
876 : gdpusch 1.1 my $loc = &location( $ftr );
877 :     $loc or return undef;
878 :    
879 :     my $loc0 = $loc;
880 :     my $complement = ( $loc =~ s/^complement\((.*)\)$/$1/ );
881 :     $loc =~ s/^join\((.*)\)$/$1/;
882 :     my @spans = split /,/, $loc;
883 :     if ( grep { ! /^<?\d+\.\.>?\d+$/ } @spans )
884 :     {
885 :     print STDERR "*** Feature location parse error: $loc0\n";
886 :     return undef;
887 :     }
888 :    
889 :     my $partial_5 = $spans[ 0] =~ s/^<//;
890 :     my $partial_3 = $spans[-1] =~ s/\.\.>/../;
891 :     ( $partial_5, $partial_3 ) = ( $partial_3, $partial_5 ) if $complement;
892 :    
893 :     my $seq = join( '', map { extract_span( $dnaR, $_ ) } @spans );
894 :     $seq = gjoseqlib::complement_DNA_seq( $seq ) if $complement;
895 :    
896 :     # Sequences that run off the end can start at other than the first
897 :     # nucleotide of a codon.
898 :    
899 :     my $qual = &qualifiers( $ftr );
900 :     my $codon_start = $qual->{ codon_start } ? $qual->{ codon_start }->[0] : 1;
901 :     $seq = substr( $seq, $codon_start-1 ) if $codon_start > 1;
902 :    
903 :     wantarray ? ( $seq, $partial_5, $partial_3 ) : $seq;
904 :     }
905 :    
906 : golsen 1.7
907 : gdpusch 1.1 sub extract_span
908 :     {
909 :     my ( $dnaR, $span ) = @_;
910 :     my ( $beg, $end ) = $span =~ /^<?(\d+)\.\.>?(\d+)$/;
911 :     ( $beg > 0 ) && ( $beg <= $end ) && ( $end <= length( $$dnaR ) ) or return '';
912 :    
913 :     substr( $$dnaR, $beg-1, $end-$beg+1 );
914 :     }
915 :    
916 :    
917 :     #-------------------------------------------------------------------------------
918 : golsen 1.7 # Get the location of a feature
919 :     #
920 :     # $ftr_location = location( $ftr ) # Returns empty string on failure.
921 :     #
922 :     # Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
923 :     #-------------------------------------------------------------------------------
924 :     sub location
925 :     {
926 :     my ( $ftr ) = @_;
927 :    
928 :     ( defined( $ftr )
929 :     && ( ref( $ftr ) eq 'ARRAY' )
930 :     && ( @$ftr > 1 ) )
931 :     ? $ftr->[-2]
932 :     : '';
933 :     }
934 :    
935 :    
936 :     #-------------------------------------------------------------------------------
937 : golsen 1.3 # Identify features with partial 5' or 3' ends.
938 :     #
939 :     # $partial_5_prime = partial_5_prime( $ftr )
940 :     # $partial_3_prime = partial_3_prime( $ftr )
941 :     #
942 : golsen 1.7 # Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
943 : golsen 1.3 #-------------------------------------------------------------------------------
944 :     sub partial_5_prime
945 :     {
946 :     my $ftr = shift or return undef;
947 :     my $loc = &location( $ftr ) or return undef;
948 :     my $complement = ( $loc =~ s/^complement\((.*)\)$/$1/ );
949 :     $loc =~ s/^join\((.*)\)$/$1/;
950 :     my @spans = split /,/, $loc;
951 :    
952 :     $complement ? $spans[-1] =~ /\.\.>/ : $spans[0] =~ /^</;
953 :     }
954 :    
955 :    
956 :     sub partial_3_prime
957 :     {
958 :     my $ftr = shift or return undef;
959 :     my $loc = &location( $ftr ) or return undef;
960 :     my $complement = ( $loc =~ s/^complement\((.*)\)$/$1/ );
961 :     $loc =~ s/^join\((.*)\)$/$1/;
962 :     my @spans = split /,/, $loc;
963 :    
964 :     $complement ? $spans[0] =~ /^</ : $spans[-1] =~ /\.\.>/;
965 :     }
966 :    
967 :    
968 :     #-------------------------------------------------------------------------------
969 : golsen 1.7 # Get the qualifier hash for a feature:
970 : gdpusch 1.1 #
971 : golsen 1.7 # \%ftr_qualifiers = qualifiers( $ftr )
972 : gdpusch 1.1 #
973 : golsen 1.7 # Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
974 :     # Returns empty hash reference on failure.
975 :     # Note that this list may include _qualif_order => \@keys
976 : gdpusch 1.1 #-------------------------------------------------------------------------------
977 :     sub qualifiers
978 :     {
979 :     my ( $ftr ) = @_;
980 :     my $qual;
981 :     ( defined( $ftr )
982 :     && ( ref( $ftr ) eq 'ARRAY' )
983 :     && ( @$ftr > 1 )
984 : golsen 1.7 && defined( $qual = $ftr->[-1] )
985 : gdpusch 1.1 && ( ref( $qual ) eq 'HASH' ) )
986 :     ? $qual
987 :     : {};
988 :     }
989 :    
990 :    
991 :     #-------------------------------------------------------------------------------
992 :     #
993 :     # $gene = gene( $ftr )
994 :     # @gene_and_synonyms = gene( $ftr )
995 :     #
996 :     #-------------------------------------------------------------------------------
997 :     sub gene
998 :     {
999 :     my $qual = &qualifiers( @_ );
1000 :     my %seen;
1001 :     my @gene = grep { ! $seen{ $_ }++ }
1002 :     ( ( $qual->{ gene } ? @{ $qual->{ gene } } : () ),
1003 :     ( $qual->{ gene_synonym } ? @{ $qual->{ gene_synonym } } : () )
1004 :     );
1005 :    
1006 :     wantarray ? @gene : $gene[0];
1007 :     }
1008 :    
1009 :    
1010 :     #-------------------------------------------------------------------------------
1011 :     # Prefer protein_id as id:
1012 :     #
1013 :     # $id = CDS_id( $ftr )
1014 :     #
1015 :     #-------------------------------------------------------------------------------
1016 :     sub CDS_id
1017 :     {
1018 :     my $qual = &qualifiers( @_ );
1019 :     my $id;
1020 :    
1021 :     ( $id ) = @{ $qual->{ protein_id } } if $qual->{ protein_id };
1022 :     ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if ! $id && $qual->{ db_xref };
1023 :     ( $id ) = @{ $qual->{ locus_tag } } if ! $id && $qual->{ locus_tag };
1024 :    
1025 :     $id;
1026 :     }
1027 :    
1028 :    
1029 :     #-------------------------------------------------------------------------------
1030 :     # Prefer gi number as id:
1031 :     #
1032 :     # $id = CDS_gi_or_id( $ftr )
1033 :     #
1034 :     #-------------------------------------------------------------------------------
1035 :     sub CDS_gi_or_id
1036 :     {
1037 :     my $qual = &qualifiers( @_ );
1038 :     my $id;
1039 :    
1040 :     ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if $qual->{ db_xref };
1041 :     ( $id ) = @{ $qual->{ protein_id } } if ! $id && $qual->{ protein_id };
1042 :     ( $id ) = @{ $qual->{ locus_tag } } if ! $id && $qual->{ locus_tag };
1043 :    
1044 :     $id;
1045 :     }
1046 :    
1047 :    
1048 :     #-------------------------------------------------------------------------------
1049 :     # gi number or nothing:
1050 :     #
1051 :     # $gi = CDS_gi( $ftr )
1052 :     #
1053 :     #-------------------------------------------------------------------------------
1054 :     sub CDS_gi
1055 :     {
1056 :     my $qual = &qualifiers( @_ );
1057 :    
1058 :     my ( $id ) = map { m/^GI:(.+)$/i ? $1 : () } @{ $qual->{ db_xref } } if $qual->{ db_xref };
1059 :    
1060 :     $id;
1061 :     }
1062 :    
1063 :    
1064 :     #-------------------------------------------------------------------------------
1065 :     #
1066 :     # $product = product( $ftr )
1067 :     #
1068 :     #-------------------------------------------------------------------------------
1069 :     sub product
1070 :     {
1071 :     my $qual = &qualifiers( @_ );
1072 :     my $prod;
1073 :    
1074 :     ( $prod ) = @{ $qual->{ product } } if $qual->{ product };
1075 :     ( $prod ) = @{ $qual->{ function } } if ! $prod && $qual->{ function };
1076 :     ( $prod ) = @{ $qual->{ note } } if ! $prod && $qual->{ note };
1077 :    
1078 :     $prod;
1079 :     }
1080 :    
1081 :    
1082 :     #-------------------------------------------------------------------------------
1083 :     #
1084 :     # @EC_number = EC_number( $ftr )
1085 :     # \@EC_number = EC_number( $ftr )
1086 :     #
1087 :     #-------------------------------------------------------------------------------
1088 :     sub EC_number
1089 :     {
1090 :     my $qual = &qualifiers( @_ );
1091 :     my @EC = $qual->{ EC_number } ? @{ $qual->{ EC_number } } : ();
1092 :    
1093 :     wantarray ? @EC : \@EC;
1094 :     }
1095 :    
1096 :    
1097 :     #-------------------------------------------------------------------------------
1098 :     # This is the in situ translation. Will extract from the DNA sequence if
1099 :     # supplied.
1100 :     #
1101 :     # $translation = CDS_translation( $ftr )
1102 :     # $translation = CDS_translation( $ftr, $dna )
1103 :     # $translation = CDS_translation( $ftr, \$dna )
1104 : golsen 1.3 # $translation = CDS_translation( $ftr, $entry )
1105 :     #
1106 : gdpusch 1.1 #
1107 :     #-------------------------------------------------------------------------------
1108 :     sub CDS_translation
1109 :     {
1110 : golsen 1.3 my ( $ftr, $dna ) = @_;
1111 : gdpusch 1.1 my $qual = &qualifiers( $ftr );
1112 :    
1113 :     return $qual->{ translation }->[0] if $qual->{ translation };
1114 :    
1115 : golsen 1.3 return undef if ! $dna;
1116 : gdpusch 1.1
1117 :     my $have_lib = 0;
1118 :     eval { require gjoseqlib; $have_lib = 1; };
1119 : golsen 1.3 return undef if ! $have_lib;
1120 :    
1121 :     my $CDS_dna = ftr_dna( $dna, $ftr ) or return undef;
1122 :     my $pep = gjoseqlib::translate_seq( $CDS_dna, ! partial_5_prime( $ftr ) );
1123 :     $pep =~ s/\*$// if $pep;
1124 : gdpusch 1.1
1125 : golsen 1.3 $pep;
1126 : gdpusch 1.1 }
1127 :    
1128 :    
1129 :     #===============================================================================
1130 :     # Utilities for locations and location strings.
1131 :     #===============================================================================
1132 :     # Convert GenBank location to a SEED location.
1133 :     #
1134 :     # $loc = genbank_loc_2_seed( $acc, $loc )
1135 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_seed( $acc, $loc )
1136 :     #
1137 :     #-------------------------------------------------------------------------------
1138 :     sub genbank_loc_2_seed
1139 :     {
1140 :     my ( $acc, $loc ) = @_;
1141 :     $acc && $loc or return undef;
1142 :     genbank_loc_2_string( $acc, $loc, 'seed' );
1143 :     }
1144 :    
1145 :    
1146 :     #-------------------------------------------------------------------------------
1147 :     # Convert GenBank location to a Sapling location.
1148 :     #
1149 :     # $loc = genbank_loc_2_sapling( $acc, $loc )
1150 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_sapling( $acc, $loc )
1151 :     #
1152 :     #-------------------------------------------------------------------------------
1153 :     sub genbank_loc_2_sapling
1154 :     {
1155 :     my ( $acc, $loc ) = @_;
1156 :     $acc && $loc or return undef;
1157 :     genbank_loc_2_string( $acc, $loc, 'sapling' );
1158 :     }
1159 :    
1160 :    
1161 :     #-------------------------------------------------------------------------------
1162 :     # Convert GenBank location to another location format.
1163 :     # At present, only 'sapling' (D) and 'seed' are supported.
1164 :     #
1165 :     # $loc = genbank_loc_2_string( $acc, $loc, $format )
1166 :     # ( $loc, $partial_5, $partial_3 ) = genbank_loc_2_string( $acc, $loc, $format )
1167 :     #
1168 :     #-------------------------------------------------------------------------------
1169 :     sub genbank_loc_2_string
1170 :     {
1171 :     my ( $acc, $loc, $format ) = @_;
1172 :     $acc && $loc or return undef;
1173 :    
1174 :     my ( $cbdl, $partial_5, $partial_3 ) = genbank_loc_2_cbdl( $loc, $acc );
1175 :     my $str = cbdl_2_string( $cbdl, $format || 'sapling' );
1176 :    
1177 :     wantarray ? ( $str, $partial_5, $partial_3 ) : $str;
1178 :     }
1179 :    
1180 :    
1181 :     #-------------------------------------------------------------------------------
1182 :     # Convert GenBank location to a list of [contig, begin, dir, len ] locations.
1183 :     # order() is treated as join(). Nesting is allowed (unlike the standard).
1184 :     #
1185 :     # \@cbdl = genbank_loc_2_cbdl( $loc, $accession )
1186 :     # ( \@cbdl, $partial_5, $partial_3 ) = genbank_loc_2_cbdl( $loc, $accession )
1187 :     #
1188 :     # Elements are:
1189 :     #
1190 :     # (accession:)?<?\d+..>?\d+ # range of sites
1191 :     # (accession:)?<?\d+^>?\d+ # site between residues
1192 :     # (accession:)?\d+ # single residue
1193 :     # (accession:)?complement\(element\)
1194 :     # (accession:)?join\(element,element,...\)
1195 :     # (accession:)?order\(element,element,...\)
1196 :     #
1197 :     #-------------------------------------------------------------------------------
1198 :     # Paterns used in the parsing. They are in each subroutine due to a very
1199 :     # strange initialization issue.
1200 :     #
1201 :     # Because $paranthetical is self-referential, it must be declared before it is
1202 :     # defined.
1203 :     #
1204 :     # my $paranthetical;
1205 :     # $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
1206 :     # my $contigid = qr/[^\s:(),]+/;
1207 :     # my $complement = qr/(?:$contigid:)?complement$paranthetical/;
1208 :     # my $complement_parts = qr/(?:($contigid):)?complement($paranthetical)/;
1209 :     # my $join = qr/(?:$contigid:)?join$paranthetical/;
1210 :     # my $join_parts = qr/(?:($contigid):)?join($paranthetical)/;
1211 :     # my $order = qr/(?:$contigid:)?order$paranthetical/;
1212 :     # my $order_parts = qr/(?:($contigid):)?order($paranthetical)/;
1213 :     # my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
1214 :     # my $range_parts = qr/(?:($contigid):)?(<?)(\d+)\.\.(>?)(\d+)/;
1215 :     # my $site = qr/(?:$contigid:)?<?\d+^>?\d+/;
1216 :     # my $site_parts = qr/(?:($contigid):)?(<?)(\d+)^(>?)(\d+)/;
1217 :     # my $position = qr/(?:$contigid:)?\d+/;
1218 :     # my $position_parts = qr/(?:($contigid):)?(\d+)/;
1219 :     # my $element = qr/$range|$position|$complement|$join|$order/;
1220 :     # my $elementlist = qr/$element(?:,$element)*/;
1221 :     #
1222 :     #-------------------------------------------------------------------------------
1223 :    
1224 :     sub genbank_loc_2_cbdl
1225 :     {
1226 :     my ( $loc, $acc ) = @_;
1227 :    
1228 :     my $contigid = qr/[^\s:(),]+/;
1229 :    
1230 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
1231 :     return gb_loc_range( $loc, $acc ) if $loc =~ /^$range$/;
1232 :    
1233 :     # This cannot by part of any other format except complement
1234 :     my $site = qr/(?:$contigid:)?<?\d+\^>?\d+/;
1235 :     return gb_loc_site( $loc, $acc ) if $loc =~ /^$site$/;
1236 :    
1237 :     my $position = qr/(?:$contigid:)?\d+/;
1238 :     return gb_loc_position( $loc, $acc ) if $loc =~ /^$position$/;
1239 :    
1240 :     my $paranthetical;
1241 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
1242 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
1243 :     return gb_loc_complement( $loc, $acc ) if $loc =~ /^$complement$/;
1244 :    
1245 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
1246 :     return gb_loc_join( $loc, $acc ) if $loc =~ /^$join$/;
1247 :    
1248 :     # Treated as a join
1249 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
1250 :     return gb_loc_order( $loc, $acc ) if $loc =~ /^$order$/;
1251 :    
1252 :     return ();
1253 :     }
1254 :    
1255 :    
1256 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1257 :     # A range of positions, with optional accession number prefix, optional less
1258 :     # than first position (5' partial), begin, .., optional greater than end
1259 :     # position (3' partial), and end position.
1260 :     #
1261 :     # (\S+:)?<?\d+\.\.>?\d+
1262 :     #
1263 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_range( $loc, $acc )
1264 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1265 :     sub gb_loc_range
1266 :     {
1267 :     my ( $loc, $acc ) = @_;
1268 :    
1269 :     my $contigid = qr/[^\s:(),]+/;
1270 :     my $range_parts = qr/(?:($contigid):)?(<?)(\d+)\.\.(>?)(\d+)/;
1271 :     my ( $acc2, $p5, $beg, $p3, $end ) = $loc =~ /^$range_parts$/;
1272 :     $beg && $end or return ();
1273 :     $acc2 ||= $acc;
1274 :    
1275 :     # GenBank standard is always $beg <= $end. We will relax that.
1276 :    
1277 :     ( [ [ $acc2, $beg, (($end>=$beg)?'+':'-'), abs($end-$beg)+1 ] ], $p5, $p3 );
1278 :     }
1279 :    
1280 :    
1281 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1282 :     # A range of positions, with optional accession number prefix, optional less
1283 :     # than first position (5' partial), begin, ^, optional greater than end
1284 :     # position (3' partial), and end position.
1285 :     #
1286 :     # (\S+:)?<?\d+^>?\d+
1287 :     #
1288 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_site( $loc, $acc )
1289 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1290 :     sub gb_loc_site
1291 :     {
1292 :     my ( $loc, $acc ) = @_;
1293 :    
1294 :     my $contigid = qr/[^\s:(),]+/;
1295 :     my $site_parts = qr/(?:($contigid):)?(<?)(\d+)\^(>?)(\d+)/;
1296 :     my ( $acc2, $p5, $beg, $p3, $end ) = $loc =~ /^$site_parts$/;
1297 :     $beg && $end or return ();
1298 :     $acc2 ||= $acc;
1299 :    
1300 :     # GenBank standard is always $beg <= $end. We will relax that.
1301 :    
1302 :     ( [ [ $acc2, $beg, (($end>=$beg)?'+':'-'), abs($end-$beg)+1 ] ], $p5, $p3 );
1303 :     }
1304 :    
1305 :    
1306 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1307 :     # A singe position, with optional accession number prefix.
1308 :     #
1309 :     # (\S+:)?\d+
1310 :     #
1311 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_position( $loc, $acc )
1312 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1313 :     sub gb_loc_position
1314 :     {
1315 :     my ( $loc, $acc ) = @_;
1316 :    
1317 :     my $contigid = qr/[^\s:(),]+/;
1318 :     my $position_parts = qr/(?:($contigid):)?(\d+)/;
1319 :     my ( $acc2, $beg ) = $loc =~ /^$position_parts$/;
1320 :     $beg or return ();
1321 :    
1322 :     ( [ [ $acc2 || $acc, $beg, '+', 1 ] ], '', '' );
1323 :     }
1324 :    
1325 :    
1326 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1327 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_complement( $loc, $acc )
1328 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1329 :     sub gb_loc_complement
1330 :     {
1331 :     my ( $loc, $acc ) = @_;
1332 :    
1333 :     my $paranthetical;
1334 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
1335 :     my $contigid = qr/[^\s:(),]+/;
1336 :     my $complement_parts = qr/(?:($contigid):)?complement($paranthetical)/;
1337 :    
1338 :     my ( $acc2, $loc2 ) = $loc =~ /^$complement_parts$/;
1339 :     $loc2 && $loc2 =~ s/^\(// && $loc2 =~ s/\)$// or return ();
1340 :     my ( $locs, $p5, $p3 ) = genbank_loc_2_cbdl( $loc2, $acc2 || $acc );
1341 :     $locs && ref( $locs ) eq 'ARRAY' && @$locs or return ();
1342 :    
1343 :     ( [ map { complement_cbdl( @$_ ) } reverse @$locs ], $p3, $p5 );
1344 :     }
1345 :    
1346 :    
1347 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1348 :     #
1349 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_join( $loc, $acc )
1350 :     #
1351 :     # There is no warning about partial sequences internal to list.
1352 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1353 :     sub gb_loc_join
1354 :     {
1355 :     my ( $loc, $acc ) = @_;
1356 :    
1357 :     my $paranthetical;
1358 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
1359 :     my $contigid = qr/[^\s:(),]+/;
1360 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
1361 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
1362 :     my $join_parts = qr/(?:($contigid):)?join($paranthetical)/;
1363 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
1364 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
1365 :     my $position = qr/(?:$contigid:)?\d+/;
1366 :     my $element = qr/$range|$position|$complement|$join|$order/;
1367 :     my $elementlist = qr/$element(?:,$element)*/;
1368 :    
1369 :     my ( $acc2, $locs ) = $loc =~ /^$join_parts$/;
1370 :     $locs && $locs =~ s/^\(// && $locs =~ s/\)$//
1371 :     && $locs =~ /^$elementlist$/
1372 :     or return ();
1373 :     $acc2 ||= $acc;
1374 :    
1375 :     my @elements = map { [ genbank_loc_2_cbdl( $_, $acc2 ) ] }
1376 :     $locs =~ m/($element)/g;
1377 :     @elements or return ();
1378 :    
1379 :     ( [ map { @{ $_->[0] } } @elements ], $elements[0]->[1], $elements[-1]->[2] );
1380 :     }
1381 :    
1382 :    
1383 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1384 :     # Ordered list is treated as a join:
1385 :     #
1386 :     # ( \@cbdl_list, $partial_5, $partial_3 ) = gb_loc_order( $loc, $acc )
1387 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1388 :     sub gb_loc_order
1389 :     {
1390 :     my ( $loc, $acc ) = @_;
1391 :    
1392 :     my $paranthetical;
1393 :     $paranthetical = qr/\([^()]*(?:(??{$paranthetical})[^()]*)*\)/;
1394 :     my $contigid = qr/[^\s:(),]+/;
1395 :     my $complement = qr/(?:$contigid:)?complement$paranthetical/;
1396 :     my $join = qr/(?:$contigid:)?join$paranthetical/;
1397 :     my $order = qr/(?:$contigid:)?order$paranthetical/;
1398 :     my $order_parts = qr/(?:($contigid):)?order($paranthetical)/;
1399 :     my $range = qr/(?:$contigid:)?<?\d+\.\.>?\d+/;
1400 :     my $position = qr/(?:$contigid:)?\d+/;
1401 :     my $element = qr/$range|$position|$complement|$join|$order/;
1402 :     my $elementlist = qr/$element(?:,$element)*/;
1403 :    
1404 :     my ( $acc2, $locs ) = $loc =~ /^$order_parts$/;
1405 :     $locs && $locs =~ s/^\(// && $locs =~ s/\)$//
1406 :     && $locs =~ /^$elementlist$/
1407 :     or return ();
1408 :    
1409 :     gb_loc_join( "join($locs)", $acc2 || $acc );
1410 :     }
1411 :    
1412 :    
1413 :     #-------------------------------------------------------------------------------
1414 :     # $cbdl = complement_cbdl( $contig, $beg, $dir, $len )
1415 :     # $cbdl = complement_cbdl( [ $contig, $beg, $dir, $len ] )
1416 :     #-------------------------------------------------------------------------------
1417 :     sub complement_cbdl
1418 :     {
1419 :     defined $_[0] or return ();
1420 :     my ( $contig, $beg, $dir, $len ) = ref( $_[0] ) ? @{$_[0]} : @_;
1421 :    
1422 :     ( $dir =~ /^-/ ) ? [ $contig, $beg -= $len - 1, '+', $len ]
1423 :     : [ $contig, $beg += $len - 1, '-', $len ];
1424 :     }
1425 :    
1426 :    
1427 :     #-------------------------------------------------------------------------------
1428 :     # $loc = cbdl_2_string( \@cbdl, $format )
1429 :     #-------------------------------------------------------------------------------
1430 :     sub cbdl_2_string
1431 :     {
1432 :     my ( $cbdl, $format ) = @_;
1433 :     $cbdl && ( ref( $cbdl ) eq 'ARRAY' ) or return undef;
1434 :     $format = 'sapling' if ! defined $format;
1435 : gdpusch 1.2 return cbdl_2_genbank( $cbdl ) if $format =~ m/genbank/i;
1436 : gdpusch 1.1 join( ',', map { cbdl_part_2_string( $_, $format ) } @$cbdl );
1437 :     }
1438 :    
1439 :    
1440 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1441 :     # Support function for formatting one contiguous part of location.
1442 :     #
1443 :     # $loc_part = cbdl_part_2_string( $cbdl_part, $format )
1444 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1445 :     sub cbdl_part_2_string
1446 :     {
1447 :     my ( $part, $format ) = @_;
1448 :     $part && ( ref( $part ) eq 'ARRAY' ) && ( @$part == 4 ) or return ();
1449 :     my ( $contig, $beg, $dir, $len ) = @$part;
1450 :     $dir = $dir =~ /^-/ ? '-' : '+';
1451 :    
1452 :     if ( $format =~ m/seed/i )
1453 :     {
1454 :     my $n2 = ( $dir eq '+' ) ? $beg + $len - 1 : $beg - $len + 1;
1455 :     return join( '_', $contig, $beg, $n2 );
1456 :     }
1457 :    
1458 :     # Default is sapling:
1459 :    
1460 :     return $contig . '_' . $beg . $dir . $len;
1461 :     }
1462 :    
1463 :     #-------------------------------------------------------------------------------
1464 :     # Convert a [ [ contig, begin, dir, length ], ... ] location to GenBank.
1465 :     #
1466 :     # $gb_location = cbdl_2_genbank( \@cbdl )
1467 :     # ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )
1468 :     #-------------------------------------------------------------------------------
1469 :     sub cbdl_2_genbank
1470 :     {
1471 :     my ( $cbdl, $contig ) = @_;
1472 :     $cbdl && ref( $cbdl ) eq 'ARRAY' && @$cbdl or return '';
1473 :     my @cbdl = ref( $cbdl->[0] ) ? @$cbdl : ( $cbdl );
1474 :     @cbdl or return '';
1475 :    
1476 :     my $dir = $cbdl[0]->[2];
1477 :     @cbdl = map { complement_cbdl( $_ ) } reverse @cbdl if $dir =~ /^-/;
1478 :    
1479 : gdpusch 1.2 $contig = $cbdl[0]->[0];
1480 : gdpusch 1.1 my @gb = map { cbdl_part_2_genbank( $_, $contig ) } @cbdl;
1481 :    
1482 :     my $gb = ( @gb > 1 ) ? 'join(' . join( ',', @gb ) . ')' : $gb[0];
1483 :    
1484 :     $gb = "complement($gb)" if $dir =~ /^-/;
1485 :    
1486 :     return wantarray ? ( $contig, $gb ) : $gb;
1487 :     }
1488 :    
1489 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1490 :     # Support function for formatting one contiguous part of location.
1491 :     #
1492 :     # $loc_part = cbdl_part_2_genbank( $cbdl_part, $contig )
1493 :     # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1494 :     sub cbdl_part_2_genbank
1495 :     {
1496 :     my ( $part, $contig0 ) = @_;
1497 :     $part && ( ref( $part ) eq 'ARRAY' ) && ( @$part == 4 ) or return ();
1498 :     my ( $contig, $beg, $dir, $len ) = @$part;
1499 :    
1500 :     my $gb;
1501 :     if ( $dir =~ /^-/ )
1502 :     {
1503 :     my $end = $beg - ( $len-1);
1504 :     $gb = "complement($end..$beg)";
1505 :     }
1506 :     else
1507 :     {
1508 :     my $end = $beg + ($len-1);
1509 :     $gb = "$beg..$end";
1510 :     }
1511 :    
1512 :     $gb = "$contig:$gb" if $contig0 && $contig ne $contig0;
1513 :    
1514 :     return $gb;
1515 :     }
1516 :    
1517 : golsen 1.7
1518 :    
1519 :     #===============================================================================
1520 :     # Write a GenBank entry:
1521 :     #
1522 :     # write_genbank( \%entry );
1523 :     # write_genbank( $filename, \%entry );
1524 :     # write_genbank( \*FH, \%entry );
1525 :     #
1526 :     #===============================================================================
1527 :     # Recognized record types are:
1528 :     #
1529 :     # + ACCESSION => [ Accession, ... ]
1530 :     # COMMENT => [ Comment, ... ]
1531 :     # DBLINK => [ Field, ... ]
1532 :     # DBSOURCE => [ Field, ... ]
1533 :     # DEFINITION => Definition
1534 :     # KEYWORDS => { Key_phrase => 1, Key_phrase => 1, ... }
1535 :     # + LOCUS => Locus_id # locus_id is preferred
1536 :     # * ORGANISM => Organism_name
1537 :     # ORIGIN => Description_of_sequence_origin
1538 :     # REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]
1539 :     # * SEQUENCE => Sequence
1540 :     # SOURCE => Source_string
1541 :     # TAXONOMY => [ Taxon, Taxon, ... ]
1542 :     # + VERSION => [ Version, Other_information ]
1543 :     #
1544 :     # Data that are not in record types, but supply data to build records:
1545 :     #
1546 :     # date => Date
1547 :     # geometry => linear | circular
1548 :     # gi => Entry_gi_number
1549 :     # is_protein => Bool
1550 :     # locus_id => Locus_id
1551 :     # mol_type => Mol_type
1552 :     #
1553 :     # The record types marked with a * are required.
1554 :     # At least of of the types marked with a + must also be present.
1555 :     #
1556 :     # Feature records are merged by type. Slash is removed from qualifier name.
1557 :     # Surrounding quotation marks are stripped from qualifier values.
1558 :     #
1559 :     # FEATURES => { Type => [ [ Location, { Qualifier => \@values } ],
1560 :     # [ Location, { Qualifier => \@values } ],
1561 :     # ...
1562 :     # ],
1563 :     # Type => ...
1564 :     # }
1565 :     #
1566 :     # 1 2 3 4 5 6 7 8
1567 :     # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
1568 :     # LOCUS @<<<<<<<<<<<<<<< @########## @< DNA @<<<<<<< @<< @<<<<<<<<<<
1569 :     # LOCUS CYC_HUMAN 105 aa linear PRI 15-MAR-2004
1570 :     # LOCUS NC_000909 1664970 bp DNA circular BCT 03-DEC-2005
1571 :     #
1572 :     #===============================================================================
1573 :    
1574 :    
1575 :     sub write_genbank
1576 :     {
1577 :     my ( $fh, $close );
1578 :     if ( $_[0] && ( ( ! ref( $_[0] ) ) || ( ref( $_[0] ) eq 'GLOB' ) ) )
1579 :     {
1580 :     ( $fh, $close ) = output_filehandle( shift );
1581 :     }
1582 :    
1583 :     my $entry = shift;
1584 :     $entry && ref $entry eq 'HASH' or return 0;
1585 :     $entry->{ SEQUENCE } or return 0;
1586 :     $entry->{ ORGANISM } or return 0;
1587 :    
1588 :     # Prepare the data for the forms:
1589 :    
1590 :     my $key;
1591 :    
1592 :     if ( ! $entry->{ locus_id } )
1593 :     {
1594 :     ( $key ) = grep { m/locus_id/i } keys %$entry;
1595 :     ( $key ) = grep { m/locus/i } keys %$entry if ! $key;
1596 :     ( $key ) = grep { m/contig/i } keys %$entry if ! $key;
1597 :     ( $key ) = grep { m/accession/i } keys %$entry if ! $key;
1598 :     ( $key ) = grep { m/version/i } keys %$entry if ! $key;
1599 :     my $locus_id = $key ? $entry->{ $key } : 'undefined';
1600 :     $locus_id = $locus_id->[0] if ref $locus_id;
1601 :     $locus_id =~ s/\s.*$//;
1602 :     $entry->{ locus_id } = $locus_id;
1603 :     }
1604 :    
1605 :     if ( ! $entry->{ mol_type } )
1606 :     {
1607 :     my $mol_type;
1608 :     my $is_prot = $entry->{ is_protein };
1609 :     $mol_type = ' ' if $is_prot;
1610 :     if ( ! $mol_type )
1611 :     {
1612 :     $mol_type = mol_type( \$entry->{ SEQUENCE } );
1613 :     $mol_type = ' ' if $mol_type =~ m/^prot/i;
1614 :     $entry->{ is_protein } = ( $mol_type =~ m/NA$/i ) ? 0 : 1;
1615 :     }
1616 :     $entry->{ mol_type } = $mol_type;
1617 :     }
1618 :    
1619 :     if ( ! $entry->{ DEFINITION } )
1620 :     {
1621 :     ( $key ) = grep { m/definition/i } keys %$entry;
1622 :     ( $key ) = grep { m/^def/i } keys %$entry if ! $key;
1623 :     my $def = $key ? $entry->{ $key } : $entry->{ locus_id };
1624 :     $entry->{ DEFINITION } = $def;
1625 :     }
1626 :    
1627 :     if ( ! $entry->{ ACCESSION } || ! ref $entry->{ ACCESSION } )
1628 :     {
1629 :     ( $key ) = grep { m/accession/i } keys %$entry;
1630 :     ( $key ) = grep { m/^acc/i } keys %$entry if ! $key;
1631 :     ( $key ) = grep { m/version/i } keys %$entry if ! $key;
1632 :     ( $key ) = grep { m/^ver/i } keys %$entry if ! $key;
1633 :     $key = 'locus_id' if ! $key;
1634 :     my $acc = $entry->{ $key };
1635 :     $entry->{ ACCESSION } = ref $acc ? $acc : [ $acc ];
1636 :     }
1637 :    
1638 :     if ( ! $entry->{ VERSION } || ! ref $entry->{ VERSION } )
1639 :     {
1640 :     ( $key ) = grep { m/version/i } keys %$entry;
1641 :     ( $key ) = grep { m/^ver/i } keys %$entry if ! $key;
1642 :     $key = 'ACCESSION' if ! $key;
1643 :     my $ver = $entry->{ $key };
1644 :     $entry->{ VERSION } = ref $ver ? $ver : [ $ver ];
1645 :     }
1646 :    
1647 :     if ( ! grep { m/^gi:/ } @{ $entry->{ VERSION } } )
1648 :     {
1649 :     ( $key ) = grep { m/^gi$/i } keys %$entry;
1650 :     push @{ $entry->{ VERSION } }, "GI:$entry->{$key}" if $key;
1651 :     }
1652 :    
1653 :    
1654 :     # Fill out the forms
1655 :    
1656 :     my $old_fh;
1657 :     $old_fh = select( $fh ) if ref $fh eq 'GLOB';
1658 :    
1659 :     write_locus( $entry );
1660 :     write_definition( $entry );
1661 :     write_accession( $entry );
1662 :     write_version( $entry );
1663 :     write_dblink( $entry ) if $entry->{ DBLINK };
1664 :     write_dbsource( $entry ) if $entry->{ DBSOURCE };
1665 :     write_keywords( $entry );
1666 :     write_source( $entry );
1667 :     write_organism( $entry );
1668 :     write_taxonomy( $entry );
1669 :     write_references( $entry );
1670 :     write_comment( $entry ) if $entry->{ COMMENT };
1671 :     write_features( $entry );
1672 : golsen 1.8 # write_base_count( $entry ); # Just not useful, and error prone
1673 : golsen 1.7 write_origin( $entry );
1674 :     write_sequence( $entry );
1675 :     write_end_of_entry();
1676 :    
1677 :     select( $old_fh ) if $old_fh;
1678 :    
1679 :     close( $fh ) if $fh && $close;
1680 :    
1681 :     return 1;
1682 :     }
1683 :    
1684 :    
1685 :     sub mol_type
1686 :     {
1687 :     return undef if ! defined $_[0];
1688 :     local $_ = ref $_[0] ? shift : \$_[0];
1689 :     my $n_nt = $$_ =~ tr/ACGNTUacgntu//;
1690 :     my $n = length( $$_ );
1691 :     return ( $n_nt < 0.8 * $n ) ? 'protein'
1692 :     : ( $$_ =~ tr/Tt// > $$_ =~ tr/Uu// ) ? 'DNA' : 'RNA';
1693 :     }
1694 :    
1695 :    
1696 :    
1697 :     sub write_locus
1698 :     {
1699 :     my $entry = shift;
1700 :     my ( $locus, $mol_type, $geometry, $div, $date )
1701 :     = map { $entry->{ $_ } || ' ' }
1702 :     qw( locus_id mol_type geometry division date );
1703 :    
1704 :     my $length = $entry->{ length } ||= length( $entry->{ SEQUENCE } );
1705 :    
1706 :     my $unit = ( $mol_type =~ m/NA$/i ) ? 'bp' : 'aa';
1707 :    
1708 :     my $stranded = ' ';
1709 :     ( $stranded, $mol_type ) = ( $1, $2 ) if $mol_type =~ /^(\S+-)(\S+)$/;
1710 :    
1711 :     my $key;
1712 :     ( $key ) = grep { m/^geometry/i } keys %$entry;
1713 :     ( $key ) = grep { m/^geom/i } keys %$entry if ! $key;
1714 :     $geometry = $entry->{ geometry } ||= $key ? $entry->{ $key }
1715 :     : ' ';
1716 :    
1717 :     ( $key ) = grep { m/^div/i } keys %$entry;
1718 :     $div = $entry->{ division } ||= $key ? $entry->{ $key } : 'UNK';
1719 :    
1720 :     ( $key ) = grep { m/^date/i } keys %$entry;
1721 :     ( $key ) = grep { m/date/i } keys %$entry if ! $key;
1722 :     $date = $entry->{ date } ||= $key ? $entry->{ $key }
1723 :     : ( map { chomp; uc } `date '+%d-%b-%Y'` )[0];
1724 :    
1725 :     $^A = '';
1726 :     formline( <<'End_locus', $locus, $length, $unit, $stranded, $mol_type, $geometry, $div, $date );
1727 :     LOCUS @<<<<<<<<<<<<<<<<< @######## @< @>>@<<<<<< @<<<<<<< @<< @<<<<<<<<<<
1728 :     End_locus
1729 :     print $^A;
1730 :     $^A = '';
1731 :     }
1732 :    
1733 :    
1734 :     sub write_definition
1735 :     {
1736 :     local $_ = $_[0]->{ DEFINITION } || '';
1737 :     $_ .= '.' unless /\.$/;
1738 :     output_record( 'DEFINITION', $_ );
1739 :     }
1740 :    
1741 :    
1742 :     sub output_record
1743 :     {
1744 :     my ( $type, $text ) = @_;
1745 :    
1746 :     $^A = '';
1747 :    
1748 :     formline( <<'End_of_line_1', $type, $text );
1749 :     @<<<<<<<<< ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1750 :     End_of_line_1
1751 :    
1752 :     formline( <<'End_of_continuation', $text ) if defined $text && length $text;
1753 :     ~~ ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1754 :     End_of_continuation
1755 :    
1756 :     print $^A;
1757 :     $^A = '';
1758 :     }
1759 :    
1760 :    
1761 :     sub write_continuation
1762 :     {
1763 :     local $_ = shift;
1764 :    
1765 :     $^A = '';
1766 :     formline( <<'End_of_continuation', $_ ) if defined $_ && length $_;
1767 :     ~~ ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1768 :     End_of_continuation
1769 :     print $^A;
1770 :     $^A = '';
1771 :     }
1772 :    
1773 :    
1774 :     sub write_accession
1775 :     {
1776 :     local $_ = join( ' ', @{ $_[0]->{ ACCESSION } || [] } );
1777 :     output_record( 'ACCESSION', $_ ) if length $_;
1778 :     }
1779 :    
1780 :    
1781 :     sub write_version
1782 :     {
1783 :     local $_ = join( ' ', @{ $_[0]->{ VERSION } || [] } );
1784 :     $_ =~ s/ GI:/ GI:/g;
1785 :     output_record( 'VERSION', $_ ) if length $_;
1786 :     }
1787 :    
1788 :    
1789 :     sub write_dblink
1790 :     {
1791 :     local $_ = join( "\r", @{ $_[0]->{ DBLINK } || [] } );
1792 :     output_record( 'DBLINK', $_ ) if length $_;
1793 :     }
1794 :    
1795 :    
1796 :     sub write_dbsource
1797 :     {
1798 :     local $_ = join( "\r", @{ $_[0]->{ DBSOURCE } || [] } );
1799 :     output_record( 'DBSOURCE', $_ ) if length $_;
1800 :     }
1801 :    
1802 :    
1803 :     sub write_keywords
1804 :     {
1805 :     local $_ = join( '; ', @{ $_[0]->{ KEYWORDS } || [] } );
1806 :     $_ .= '.' if ! /\.$/;
1807 :     output_record( 'KEYWORDS', $_ );
1808 :     }
1809 :    
1810 :    
1811 :     sub write_source
1812 :     {
1813 :     local $_ = $_[0]->{ SOURCE } || 'unknown';
1814 :     output_record( 'SOURCE', $_ );
1815 :     }
1816 :    
1817 :    
1818 :     sub write_organism
1819 :     {
1820 :     local $_ = $_[0]->{ ORGANISM } || 'unclasified';
1821 :     output_record( ' ORGANISM', $_ );
1822 :     }
1823 :    
1824 :    
1825 :     sub write_taxonomy
1826 :     {
1827 :     local $_ = join( '; ', @{ $_[0]->{ TAXONOMY } || [] } );
1828 :     $_ .= '.' if ! /\.$/;
1829 :     output_record( ' ', $_ );
1830 :     }
1831 :    
1832 :    
1833 :     sub write_references
1834 :     {
1835 :     my $n = 0;
1836 :     foreach ( @{ $_[0]->{ REFERENCES } || [] } )
1837 :     {
1838 :     write_reference( $_, ++$n );
1839 :     write_authors( $_ );
1840 :     write_consortium( $_ );
1841 :     write_title( $_ );
1842 :     write_journal( $_ );
1843 :     write_pubmed( $_ );
1844 :     write_remark( $_ );
1845 :     }
1846 :     }
1847 :    
1848 :    
1849 :     sub write_reference
1850 :     {
1851 :     local $_ = $_[0]->{ ref_num } || $_[1] || '1';
1852 :     $_ .= " ($_[0]->{ref_note})" if $_[0]->{ ref_note };
1853 :     output_record( 'REFERENCE', $_ );
1854 :     }
1855 :    
1856 :    
1857 :     sub write_authors
1858 :     {
1859 :     local $_ = $_[0]->{ AUTHORS };
1860 :     output_record( ' AUTHORS', $_ ) if defined $_ && length $_;
1861 :     }
1862 :    
1863 :    
1864 :     sub write_consortium
1865 :     {
1866 :     local $_ = $_[0]->{ CONSRTM };
1867 :     output_record( ' CONSRTM', $_ ) if defined $_ && length $_;
1868 :     }
1869 :    
1870 :    
1871 :     sub write_title
1872 :     {
1873 :     local $_ = $_[0]->{ TITLE };
1874 :     output_record( ' TITLE', $_ ) if defined $_ && length $_;
1875 :     }
1876 :    
1877 :    
1878 :     sub write_journal
1879 :     {
1880 :     local $_ = $_[0]->{ JOURNAL };
1881 :     output_record( ' JOURNAL', $_ ) if defined $_ && length $_;
1882 :     }
1883 :    
1884 :    
1885 :     sub write_pubmed
1886 :     {
1887 :     local $_ = $_[0]->{ PUBMED };
1888 :     output_record( ' PUBMED', $_ ) if defined $_ && length $_;
1889 :     }
1890 :    
1891 :    
1892 :     sub write_remark
1893 :     {
1894 :     local $_ = $_[0]->{ REMARK };
1895 :     output_record( ' REMARK', $_ ) if defined $_ && length $_;
1896 :     }
1897 :    
1898 :    
1899 :     sub write_comment
1900 :     {
1901 :     local $_ = join( "\r", @{ $_[0]->{ COMMENT } || [] } );
1902 :     output_record( 'COMMENT', $_ ) if length $_;
1903 :     }
1904 :    
1905 :    
1906 :     sub write_features
1907 :     {
1908 :     my $ftr_list = feature_list( $_[0] );
1909 :     if ( $ftr_list && @$ftr_list )
1910 :     {
1911 :     write_feature_header();
1912 :     foreach ( @$ftr_list ) { write_feature( $_ ) }
1913 :     }
1914 :     }
1915 :    
1916 :    
1917 :     sub write_feature_header
1918 :     {
1919 :     print "FEATURES Location/Qualifiers\n";
1920 :     }
1921 :    
1922 :    
1923 :     sub write_feature
1924 :     {
1925 :     my ( $type, $location, $qualifiers ) = @{ $_[0] };
1926 :    
1927 :     $^A = '';
1928 :    
1929 :     my ( $loc, $rest ) = wrap_location( $location );
1930 :    
1931 :     formline( <<'End_of_line_1', $type, $loc );
1932 :     @<<<<<<<<<<<<<< ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1933 :     End_of_line_1
1934 :    
1935 :    
1936 :     while ( $rest )
1937 :     {
1938 :     ( $loc, $rest ) = wrap_location( $rest );
1939 :     formline( <<'End_of_continuation', $loc ) if defined $loc && length( $loc );
1940 :     ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1941 :     End_of_continuation
1942 :     }
1943 :    
1944 :     if ( defined $qualifiers && ref $qualifiers eq 'ARRAY' )
1945 :     {
1946 :     foreach ( @$qualifiers )
1947 :     {
1948 :     add_feature_qualifier( @$_ ) if $_->[0] !~ /^_/;
1949 :     }
1950 :     }
1951 :     elsif ( defined $qualifiers && ref $qualifiers eq 'HASH' )
1952 :     {
1953 :     foreach my $key ( ordered_qualifiers( $qualifiers ) )
1954 :     {
1955 :     foreach ( @{ $qualifiers->{ $key } } ) { add_feature_qualifier( $key, $_ ) }
1956 :     }
1957 :     }
1958 :    
1959 :     print $^A;
1960 :     $^A = '';
1961 :     }
1962 :    
1963 :    
1964 :     sub wrap_location
1965 :     {
1966 :     local $_ = shift;
1967 :     return ( $_, '' ) if length( $_ ) <= 58;
1968 :    
1969 :     return ( $1, $2 ) if m/^(.{1,57},)(.*)$/;
1970 :     return ( $1, $2 ) if m/^(.{1,56}\.\.)(.*)$/;
1971 :     return ( $1, $2 ) if m/^(.{1,57}\()(.*)$/;
1972 :     return $_;
1973 :     }
1974 :    
1975 :    
1976 :     sub ordered_qualifiers
1977 :     {
1978 :     my $qualifiers = shift;
1979 :    
1980 :     if ( ! $qualifiers->{ '_qualif_order' } )
1981 :     {
1982 :     @{ $qualifiers->{ '_qualif_order' } } = sort { qualifier_priority( $a ) <=> qualifier_priority( $b )
1983 :     || lc $a cmp lc $b
1984 :     }
1985 :     grep { ! /^_/ }
1986 :     keys %$qualifiers;
1987 :     }
1988 :    
1989 :     @{ $qualifiers->{ '_qualif_order' } };
1990 :     }
1991 :    
1992 :    
1993 :     sub qualifier_priority { $qualifier_order{ lc $_[0] } || 999999 }
1994 :    
1995 :    
1996 :     sub add_feature_qualifier
1997 :     {
1998 :     my ( $key, $value ) = @_;
1999 :    
2000 :     my $qualif;
2001 :     if ( ! defined $value || $value eq '' )
2002 :     {
2003 :     $qualif = "/$key";
2004 :     }
2005 :     else
2006 :     {
2007 :     if ( $value !~ /^\d+$/ ) { $value =~ s/"/""/g; $value = qq("$value"); }
2008 :     $qualif = "/$key=$value";
2009 :     }
2010 :    
2011 :     formline( <<'End_of_continuation', $qualif ) if defined $qualif && length $qualif;
2012 :     ~~ ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2013 :     End_of_continuation
2014 :     }
2015 :    
2016 :    
2017 : golsen 1.8 # 2011/01/23 -- GJO
2018 :     #
2019 :     # Not required, and introduces more problems than it solves.
2020 :     #
2021 : golsen 1.7 sub write_base_count
2022 :     {
2023 : golsen 1.8 return;
2024 :     #
2025 :     # my $entry = shift;
2026 :     # return if $entry->{ is_protein };
2027 :     # return if $entry->{ mol_type } =~ /\S/ && $entry->{ mol_type } !~ m/NA$/i;
2028 :     #
2029 :     # my ( $a, $c, $g, $t, $o );
2030 :     # my $counts = $entry->{ BASE_COUNT };
2031 :     # if ( $counts && ref $counts eq 'HASH' && keys %$counts)
2032 :     # {
2033 :     # ( $a, $c, $g, $t, $o ) = map { $counts->{ $_ } || 0 }
2034 :     # qw( a c g t other );
2035 :     # $a || $c || $g || $t || $o or return;
2036 :     # }
2037 :     # else
2038 :     # {
2039 :     # $entry->{ SEQUENCE } or return;
2040 :     # my $s = \$entry->{ SEQUENCE };
2041 :     # $a = $$s =~ tr/Aa//;
2042 :     # $c = $$s =~ tr/Cc//;
2043 :     # $g = $$s =~ tr/Gg//;
2044 :     # $t = $$s =~ tr/Tt//;
2045 :     # $o = length( $$s ) - ( $a + $c + $g + $t );
2046 :     # $a || $c || $g || $t || $o or return;
2047 :     #
2048 :     # $entry->{ BASE_COUNT } = { a => $a,
2049 :     # c => $c,
2050 :     # g => $g,
2051 :     # t => $t
2052 :     # };
2053 :     # $entry->{ BASE_COUNT }->{ other } = $o if $o;
2054 :     # }
2055 :     #
2056 :     # printf "BASE COUNT %6d a %6d c %6d g %6d t%s\n",
2057 :     # $a, $c, $g, $t, $o ? sprintf( ' %6d other', $o ) : '';
2058 : golsen 1.7 }
2059 :    
2060 :    
2061 :     sub write_origin
2062 :     {
2063 :     my $entry = shift;
2064 :     my $info = $entry && defined $entry->{ ORIGIN } ? $entry->{ ORIGIN } : '';
2065 :     print "ORIGIN $info\n"; # Use explicit print to get blank padding
2066 :     }
2067 :    
2068 :    
2069 :     sub write_sequence
2070 :     {
2071 :     my $entry = shift;
2072 :     $entry->{ SEQUENCE } or return;
2073 :    
2074 :     my $start = 1;
2075 :     foreach ( $entry->{ SEQUENCE } =~ m/(.{1,60})/g )
2076 :     {
2077 :     print join( ' ', sprintf( '%9d', $start ), m/(.{1,10})/g ), "\n";
2078 :     $start += 60;
2079 :     }
2080 :     }
2081 :    
2082 :    
2083 :     sub write_end_of_entry
2084 :     {
2085 :     print "//\n";
2086 :     }
2087 :    
2088 :    
2089 : gdpusch 1.1 #===============================================================================
2090 :     # Helper function for defining an input filehandle:
2091 :     #
2092 :     # filehandle is passed through
2093 : golsen 1.7 # string is taken as file name to be opened
2094 :     # undef or "" defaults to STDIN
2095 : gdpusch 1.1 #
2096 :     # \*FH = input_filehandle( $file );
2097 :     # ( \*FH, $close ) = input_filehandle( $file );
2098 :     #
2099 :     #===============================================================================
2100 :     sub input_filehandle
2101 :     {
2102 :     my $file = shift;
2103 :    
2104 :     # Null string or undef
2105 :    
2106 :     if ( ! defined( $file ) || ( $file eq '' ) )
2107 :     {
2108 :     return wantarray ? ( \*STDIN, 0 ) : \*STDIN;
2109 :     }
2110 :    
2111 :     # FILEHANDLE?
2112 :    
2113 :     if ( ref( $file ) eq "GLOB" )
2114 :     {
2115 :     return wantarray ? ( $file, 0 ) : $file;
2116 :     }
2117 :    
2118 :     # File name
2119 :    
2120 :     if ( ! ref( $file ) )
2121 :     {
2122 : golsen 1.7 -f $file or die "Could not find input file '$file'.\n";
2123 : gdpusch 1.1 my $fh;
2124 : golsen 1.7 open( $fh, "<$file" ) || die "Could not open '$file' for input.\n";
2125 : gdpusch 1.1 return wantarray ? ( $fh, 1 ) : $fh;
2126 :     }
2127 :    
2128 :     return wantarray ? ( \*STDIN, undef ) : \*STDIN;
2129 :     }
2130 :    
2131 :    
2132 : golsen 1.7 #===============================================================================
2133 :     # Helper function for defining an output filehandle:
2134 :     #
2135 :     # filehandle is passed through
2136 :     # string is taken as file name to be opened
2137 :     # undef or "" defaults to STDOUT
2138 :     #
2139 :     # \*FH = output_filehandle( $file );
2140 :     # ( \*FH, $close ) = output_filehandle( $file );
2141 :     #
2142 :     #===============================================================================
2143 :     sub output_filehandle
2144 :     {
2145 :     my $file = shift;
2146 :    
2147 :     # Null string or undef
2148 :    
2149 :     if ( ! defined( $file ) || ( $file eq '' ) )
2150 :     {
2151 :     return wantarray ? ( \*STDOUT, 0 ) : \*STDOUT;
2152 :     }
2153 :    
2154 :     # FILEHANDLE?
2155 :    
2156 :     if ( ref( $file ) eq "GLOB" )
2157 :     {
2158 :     return wantarray ? ( $file, 0 ) : $file;
2159 :     }
2160 :    
2161 :     # File name
2162 :    
2163 :     if ( ! ref( $file ) )
2164 :     {
2165 :     my $fh;
2166 :     open( $fh, ">$file" ) || die "Could not open '$file' for output.\n";
2167 :     return wantarray ? ( $fh, 1 ) : $fh;
2168 :     }
2169 :    
2170 :     return wantarray ? ( \*STDOUT, undef ) : \*STDOUT;
2171 :     }
2172 :    
2173 :    
2174 : gdpusch 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3