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

Annotation of /FigKernelPackages/gjogenbank.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3