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

Annotation of /FigKernelPackages/gjogenbank.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3