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

Diff of /FigKernelPackages/gjogenbank.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.6, Mon Dec 6 20:03:56 2010 UTC revision 1.7, Sun Jan 2 18:54:32 2011 UTC
# Line 25  Line 25 
25  #  Each entry is a hash with key value pairs, some of which have special  #  Each entry is a hash with key value pairs, some of which have special
26  #  processing.  #  processing.
27  #  #
 #     Key => Value  
 #  
28  #     ACCESSION  => [ Accession, ... ]  #     ACCESSION  => [ Accession, ... ]
29  #     DATE       =>   Date  #     COMMENT    => [ Comment_line, ... ]   # Original wrap is preserved
30    #     DBLINK     => [ Field, ... ]
31    #     DBSOURCE   => [ Field, ... ]
32  #     DEFINITION =>   Definition  #     DEFINITION =>   Definition
33  #     GEOMETRY   =>   Geometry  #     KEYWORDS   => [ Keyword, ... ]
34  #     GI         =>   GI_number  #     LOCUS      =>   Locus_id
 #     KEYWORDS   => { Key_phrase => 1, Key_phrase => 1, ... }  
 #     LOCUS      =>   Locus  
35  #     ORGANISM   =>   Organism_name  #     ORGANISM   =>   Organism_name
36    #     ORIGIN     =>   Description_of_sequence_origin
37  #     REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]  #     REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]
38  #     SEQUENCE   =>   Sequence  #     SEQUENCE   =>   Sequence
39  #     SOURCE     =>   Source_string  #     SOURCE     =>   Source_string
40  #     TAXONOMY   => [ Taxon, Taxon, ... ]  #     TAXONOMY   => [ Taxon, Taxon, ... ]
41  #     VERSION    => [ Version, Other_information ]  #     VERSION    => [ Version, Other_information ]
42  #  #
43    #  Data that are more derived or parsed:
44    #
45    #     date       =>   Date
46    #     division   =>   Division
47    #     ftr_list   => [ [ type, loc, quals ], ... ]
48    #     geometry   =>   linear | circular
49    #     gi         =>   Entry_gi_number
50    #     is_protein =>   Bool
51    #     key_index  => { Keyword => 1, ... }
52    #     locus_id   =>   Locus_id
53    #     mol_type   =>   Mol_type
54    #
55  #  Feature records are merged by type.  Slash is removed from qualifier name.  #  Feature records are merged by type.  Slash is removed from qualifier name.
56  #  Surrounding quotation marks are stripped from qualifier values.  #  Surrounding quotation marks are stripped from qualifier values.
57  #  #
# Line 113  Line 124 
124  #   ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )  #   ( $contig, $gb_location ) = cbdl_2_genbank( \@cbdl )
125  #  #
126  #===============================================================================  #===============================================================================
127    #  Write a GenBank entry:
128    #
129    #     write_genbank(             \%entry );
130    #     write_genbank( $file_name, \%entry );
131    #     write_genbank( \*FH,       \%entry );
132    #
133    #===============================================================================
134    
135  use strict;  use strict;
136    
137    use Data::Dumper;
138    
139  require Exporter;  require Exporter;
140  our @ISA    = qw( Exporter );  our @ISA    = qw( Exporter );
141  our @EXPORT = qw( parse_genbank  our @EXPORT = qw( parse_genbank
142                    parse_next_genbank                    parse_next_genbank
143                    feature_types                    feature_types
144                    features_of_type                    features_of_type
145                      feature_list
146                      ftr_dna
147    
148                      genbank_loc_2_seed
149                      genbank_loc_2_sapling
150                      genbank_loc_2_cbdl
151                      cbdl_2_genbank
152    
153                      write_genbank
154                  );                  );
155    
156  use Data::Dumper;  
157    #  An approximate ordering of the common qualifiers in GenBank feature table
158    #  entries:
159    
160    my %qualifier_order;
161    {
162        my $tmp_n = 1;
163        %qualifier_order = map { $_ => $tmp_n++ }
164                           qw( organism
165                               sub_species
166                               sex
167                               mating_type
168                               chromosome
169                               macronuclear
170                               plasmid
171                               organelle
172                               strain
173                               sub_strain
174                               cultivar
175                               serotype
176                               serovar
177                               variety
178                               isolate
179                               tissue_type
180                               cell_type
181                               dev_stage
182                               cell_line
183                               rearranged
184                               clone_lib
185                               clone
186                               sub_clone
187                               host
188                               lab_host
189                               mol_type
190                               direction
191    
192                               rpt_type
193                               rpt_family
194    
195                               operon
196                               gene
197                               gene_synonym
198                               allele
199                               locus_tag
200                               codon_start
201                               transl_table
202    
203                               anticodon
204                               product
205                               function
206                               EC_number
207                               protein_id
208                               locus_peptide
209                               pseudo
210    
211                               note
212                               db_xref
213    
214                               translation
215                            );
216    }
217    
218    
219  #===============================================================================  #===============================================================================
220  #  #
# Line 185  Line 275 
275      return $entry;      return $entry;
276  }  }
277    
   
278  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
279  #  If it should be necessary to close a stream openned by parse_next_genbank()  #  If it should be necessary to close a stream openned by parse_next_genbank()
280  #  before it reaches the end-of-file, this will do it.  #  before it reaches the end-of-file, this will do it.
# Line 215  Line 304 
304  sub parse_one_genbank_entry  sub parse_one_genbank_entry
305  {  {
306      my $fh = shift;      my $fh = shift;
307        local $_;
308    
309      my $state = 0;      my $state = 0;
310      if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }      if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
# Line 226  Line 316 
316      # -1 = Error      # -1 = Error
317    
318      my %entry = ();      my %entry = ();
     my @sequence;  
319    
320  #          1         2         3         4         5         6         7         8  #          1         2         3         4         5         6         7         8
321  # 12345678901234567890123456789012345678901234567890123456789012345678901234567890  # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
# Line 238  Line 327 
327          if ( s/^LOCUS\s+// )          if ( s/^LOCUS\s+// )
328          {          {
329              my @parts = split;              my @parts = split;
330              $entry{ LOCUS }    = shift @parts;              $entry{ locus_id } = $entry{ LOCUS } = shift @parts;
331              $entry{ DATE }     = pop @parts if $parts[-1] =~ m/^\d+-[A-Z][A-Z][A-Z]-\d+$/i;              $entry{ date }     = pop @parts if $parts[-1] =~ m/^\d+-[A-Z][A-Z][A-Z]-\d+$/i;
332              $entry{ DIVISION } = pop @parts;              $entry{ division } = pop @parts if $parts[-1] =~ m/^\S\S\S$/;
333              $entry{ GEOMETRY } = pop @parts if $parts[-1] =~ m/^(lin|circ)/i;              $entry{ geometry } = pop @parts if $parts[-1] =~ m/^(lin|circ)/i;
334              $entry{ MOL_TYPE } = pop @parts if $parts[-1] =~ m/na$/i;              $entry{ mol_type } = pop @parts if $parts[-1] =~ m/na$/i;
335              $state = 1;              $state = 1;
336          }          }
337          if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }          if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
# Line 262  Line 351 
351          elsif ( /^ORIGIN / )          elsif ( /^ORIGIN / )
352          {          {
353              $state = 3;              $state = 3;
354                last;
355          }          }
356    
357          elsif ( /^REFERENCE / )          elsif ( /^REFERENCE / )
# Line 272  Line 362 
362              defined() or $state = -1;              defined() or $state = -1;
363          }          }
364    
365          elsif ( /^(..........)  (.*\S)\s*$/ )  # Any other keyword          elsif ( /^(.{10})  +(.*\S)\s*$/ )  # Any other keyword
366          {          {
367              my ( $tag, $value ) = ( $1, $2 );              my ( $tag, @value ) = ( uc( $1 ), $2 );
368              $tag =~ s/^ +//;              $tag =~ s/^ +| +$//g;
             $tag =~ s/ +$//;  
369    
370              # Merge continuations:              # Merge continuations:
371    
             my $sep = ( $tag eq 'ORGANISM' ) ? "\t" :  
                       ( $tag eq 'COMMENT'  ) ? "\n" : ' ';  
   
372              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
373              while ( $state >= 0 && s/^ {12}/$sep/ )  
374                while ( $state >= 0 && s/^ {12}\s*// )
375              {              {
376                  $value .= $_ if length() > 1;                  s/\s+$//;
377                  $sep = ' ' if $sep eq "\t";                  push @value, $_;
378                  if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }                  if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
379              }              }
380    
381              #  Special case formats              #  Special case formats
382    
383                if ( $tag eq 'COMMENT' )
384                {
385                    push @{ $entry{ COMMENT } }, @value;
386                }
387    
388                elsif ( $tag eq 'DBLINK' )
389                {
390                    my $data = '';
391                    my @dbs;
392                    foreach ( @value )
393                    {
394                        if ( $data && /: / )
395                        {
396                            push @dbs, $data;
397                            $data = '';
398                        }
399                        $data .= $_;
400                    }
401                    push @dbs, $data  if length $data;
402    
403                    $entry{ DBLINK } = \@dbs  if @dbs;
404                }
405    
406                elsif ( $tag eq 'DBSOURCE' )
407                {
408                    my $data = '';
409                    my @dbs;
410                    foreach ( @value )
411                    {
412                        if ( $data && /: / )
413                        {
414                            push @dbs, $data;
415                            $data = '';
416                        }
417                        $data .= $_;
418                    }
419                    push @dbs, $data  if length $data;
420    
421                    $entry{ DBSOURCE } = \@dbs  if @dbs;
422                }
423    
424                elsif ( $tag eq 'ORGANISM' )
425                {
426                    my $org = shift @value;
427                    $entry{ ORGANISM } = $org;
428                    my $tax = @value ? join( ' ', @value ) : '';
429                    $tax =~ s/\s*\.$//;
430                    $entry{ TAXONOMY } = [ split /; */, $tax ] if $tax;
431                }
432    
433                else
434                {
435                    my $imax = @value - 2;
436                    if ( $imax > 0 )
437                    {
438                        foreach ( @value[ 0 .. $imax ] ) { $_ .= ' ' if ! /-$/ }
439                    }
440                    my $value = join( '', @value );
441    
442              if ( $tag eq 'ACCESSION' )              if ( $tag eq 'ACCESSION' )
443              {              {
444                  $entry{ ACCESSION } = [ split / +/, $value ];                  $entry{ ACCESSION } = [ split / +/, $value ];
445              }              }
446              elsif ( $tag eq 'VERSION' )              elsif ( $tag eq 'VERSION' )
447              {              {
448                  my ( $gi ) = $value =~ s/ +GI:(\d+)//;                      if ( $value =~ / +GI:(\d+)/ )
449                  $entry{ GI } = $gi if $gi;                      {
450                            $entry{ gi } = $1;
451                            $value =~ s/ +GI:\d+//;
452                        }
453                  $entry{ VERSION } = [ split / +/, $value ];                  $entry{ VERSION } = [ split / +/, $value ];
454              }              }
455              elsif ( $tag eq 'KEYWORDS' )              elsif ( $tag eq 'KEYWORDS' )
456              {              {
457                  $value =~ s/\s*\.$//;                      $value =~ s/\s*\.\s*$//;
458                  $entry{ KEYWORDS } = { map { $_ => 1 } split /; */, $value } if $value;                      if ( $value )
             }  
             elsif ( $tag eq 'ORGANISM' )  
459              {              {
460                  $value =~ s/\.$//;                          my @keys = split /; */, $value;
461                  my ( $org, $tax ) = split /\t/, $value;                          $entry{ KEYWORDS  } = \@keys;
462                  $entry{ ORGANISM } = $org;                          $entry{ key_index } = { map { $_ => 1 } @keys };
                 $entry{ TAXONOMY } = [ split /; */, $tax ] if $tax;  
463              }              }
464                    }
465    
466                    #  Generic case:
467    
468              else              else
469              {              {
470                  $entry{ $tag } = $value;                  $entry{ $tag } = $value;
471              }              }
472                }
473    
474              # To know that we are at end of continuations, we must have              # To know that we are at end of continuations, we must have
475              # read another line.              # read another line.
# Line 332  Line 483 
483          }          }
484      }      }
485    
486      #  Reading the features requires merging continuations, then dealing      #  Reading FEATURES requires merging continuations, then dealing
487      #  with the data:      #  with the data:
488    
489      while ( $state == 2 )      while ( $state == 2 && ( /^     (\S+)\s+(\S+)/ ) )
     {  
         if ( /^ORIGIN/ || /^BASE COUNT/ )  
         {  
             $state = 3;  
         }  
   
         elsif ( /^     (\S+)\s+(\S+)/ )  
490          {          {
491              my ( $type, $loc ) = ( $1, $2 );              my ( $type, $loc ) = ( $1, $2 );
             my ( $qualif, $value, %qualifs );  
492    
493              #  Collect the rest of the location:              #  Collect the rest of the location:
494    
495              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
496              while ( $state >= 0 && /^ {15}/ && ( $_ !~ /^\s*\/\w/ ) )          while ( $state >= 0 && /^ {15}\s*([^\/\s]\S*)\s*$/ )
497              {              {
498                  s/^ +//;              $loc .= $1;
                 $loc .= $_;  
499                  if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }                  if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
500              }              }
501    
502              #  Collect qualiiers:              #  Collect qualiiers:
503    
504            my ( $qualif, $value, %qualifs, @qualif_order );
505              while ( $state == 2 && ( $_ =~ /^\s*\/\w+/ ) )              while ( $state == 2 && ( $_ =~ /^\s*\/\w+/ ) )
506              {              {
507                  #  Qualifiers without = get an undef value (intentionally)                  #  Qualifiers without = get an undef value (intentionally)
# Line 387  Line 530 
530    
531                  if ( $qualif )                  if ( $qualif )
532                  {                  {
533                      if ( $value && $value =~ /^".*"$/ )                  if ( $value && $value =~ s/^"(.*)"$/$1/ ) { $value =~ s/""/"/g }
                     {  
                         $value =~ s/^"//;  
                         $value =~ s/"$//;  
                         $value =~ s/""/"/g;  
                     }  
   
534                      if ( $qualif eq 'translation' ) { $value =~ s/ +//g }                      if ( $qualif eq 'translation' ) { $value =~ s/ +//g }
535    
536                    push @qualif_order, $qualif  if ! $qualifs{ $qualif };
537                      push @{ $qualifs{ $qualif } }, $value;                      push @{ $qualifs{ $qualif } }, $value;
538                  }                  }
539              }              }
540    
541            $qualifs{ '_qualif_order' } = \@qualif_order if @qualif_order;
542              push @{ $entry{ FEATURES }->{ $type } }, [ $loc, \%qualifs ] if ( $type && $loc );              push @{ $entry{ FEATURES }->{ $type } }, [ $loc, \%qualifs ] if ( $type && $loc );
543            push @{ $entry{ ftr_list } },     [ $type, $loc, \%qualifs ] if ( $type && $loc );
544    
545              defined() or $state = -1;              defined() or $state = -1;
546          }          }
547    
548          elsif ( /^\s{0,4}\S/ )  # Not feature and not origin      $state = 3 if $state > 0;
549    
550        #  Only a few fields are allowed in the standard, but:
551    
552        while ( $state == 3 )
553        {
554            if ( /^ORIGIN / )
555          {          {
556              $state = -1;              $entry{ ORIGIN } = $1 if /^ORIGIN \s+(\S.*\S)\s*$/;
557                $state = 4;
558          }          }
559    
560          else          elsif ( s/^BASE COUNT\s+// )
561          {          {
562                $entry{ BASE_COUNT } = { reverse split };
563              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
564          }          }
565    
566            elsif ( /^(.{10})  (.*)$/ )  # Any other keyword
567            {
568                my ( $tag, @value ) = map { s/^ +| +$//g; $_ } ( uc( $1 ), $2 );
569    
570                # Merge continuations:
571    
572                if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
573    
574                while ( $state >= 0 && s/^ {12}\s*// )
575                {
576                    s/\s+$//;
577                    push @value, $_;
578                    if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
579      }      }
580    
581      if ( $state == 3 && s/^BASE COUNT\s+// )              #  Special case formats
582    
583                if ( $tag eq 'COMMENT' )
584      {      {
585          $entry{ BASE_COUNT } = { reverse split };                  push @{ $entry{ COMMENT } }, @value;
586          $state = ( defined( $_ = <$fh> ) && /^ORIGIN/ ) ? 3 : -1;              }
587    
588                else
589                {
590                    my $imax = @value - 2;
591                    if ( $imax > 0 )
592                    {
593                        foreach ( @value[ 0 .. $imax ] ) { $_ .= ' ' if ! /-$/ }
594                    }
595                    $entry{ $tag } = join( '', @value );
596                }
597            }
598    
599            else  # This is really a format error, but let's skip it.
600            {
601                if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
602            }
603      }      }
604    
605      #  Read the sequence:      #  Read the sequence:
606    
607      while ( $state == 3 )      while ( $state == 4 )
608      {      {
609            my @sequence;
610    
611          if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }          if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
612          while ( $state >= 0 && /^[ 0-9]{10}/ )          while ( $state >= 0 && s/^ *\d+ +// )
613          {          {
614              s/[^A-Za-z]+//g;              s/[^A-Za-z]+//g;
615              push @sequence, $_;              push @sequence, $_;
616              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }              if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
617          }          }
618          $entry{ SEQUENCE } = join "", @sequence;          $entry{ SEQUENCE } = join '', @sequence;
619    
620          $state = ( $_ eq '//' ) ? 0 : -1 if $state >= 0;          $state = ( $_ eq '//' ) ? 0 : -1 if $state >= 0;
621      }      }
# Line 449  Line 631 
631  {  {
632      my ( $line, $fh ) = @_;      my ( $line, $fh ) = @_;
633      my $state = 1;      my $state = 1;
634        my %ref = ();
635    
636        if ( $line =~ /\s*(\d+)/ )
637        {
638            $ref{ ref_num  } = $1;
639            $ref{ ref_note } = $1 if $line =~ /\s*\d+\s+\((.*)\)\s*$/;
640        }
641    
642      if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }      if ( defined( $_ = <$fh> ) ) { chomp } else { $state = -1 }
643    
644      my ( $tag, $value );      my ( $tag, $value );
     my %ref = ();  
645      while ( ( $state >= 0 ) && /^  / )      while ( ( $state >= 0 ) && /^  / )
646      {      {
647          if ( substr( $_, 0, 10 ) =~ /\S/ )          if ( substr( $_, 0, 10 ) =~ /\S/ )
# Line 489  Line 678 
678  #     @types = feature_types( $entry );  #     @types = feature_types( $entry );
679  #    \@types = feature_types( $entry );  #    \@types = feature_types( $entry );
680  #  #
 #     @ftrs = features_of_type( $entry,  @types );  
 #    \@ftrs = features_of_type( $entry,  @types );  
 #     @ftrs = features_of_type( $entry, \@types );  
 #    \@ftrs = features_of_type( $entry, \@types );  
 #  
681  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
   
682  sub feature_types  sub feature_types
683  {  {
684      my $entry = shift;      my $entry = shift;
# Line 504  Line 687 
687      my $ftrs = $entry->{ FEATURES };      my $ftrs = $entry->{ FEATURES };
688      return wantarray ? () : [] if ! ( $ftrs && ref $ftrs eq 'HASH' );      return wantarray ? () : [] if ! ( $ftrs && ref $ftrs eq 'HASH' );
689    
690      my @types = sort { lc $a cmp lc $b } keys %$ftrs;      my @types = sort { lc $a cmp lc $b } grep { ! /^_/ } keys %$ftrs;
691      wantarray ? @types : \@types;      wantarray ? @types : \@types;
692  }  }
693    
694    
695    #-------------------------------------------------------------------------------
696    #  This is really lame for more than one type because the individual features
697    #  are not marked with their types.  This should be more effecient than the
698    #  following method for finding all the features of one type.
699    #
700    #     @ftrs = features_of_type( $entry,  @types );
701    #    \@ftrs = features_of_type( $entry,  @types );
702    #     @ftrs = features_of_type( $entry, \@types );
703    #    \@ftrs = features_of_type( $entry, \@types );
704    #
705    #-------------------------------------------------------------------------------
706  sub features_of_type  sub features_of_type
707  {  {
708      my $entry = shift;      my $entry = shift;
# Line 527  Line 721 
721    
722    
723  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
724    #  Get a list of a features as triples [ $type, $loc, \%quals ] in genome
725    #  order:
726    #
727    #      @features = feature_list( $entry )
728    #     \@features = feature_list( $entry )
729    #      @features = feature_list( $entry, $type )
730    #     \@features = feature_list( $entry, $type )
731    #
732    #  Uses $entry->{ ftr_list } to cache the results (before filtering by type)
733    #-------------------------------------------------------------------------------
734    sub feature_list
735    {
736        my $entry = shift;
737    
738        if ( ! $entry->{ ftr_list } || ref $entry->{ ftr_list } ne 'ARRAY' )
739        {
740            my $features = $entry->{ FEATURES };
741    
742            #  Is it a list of [ $type, $location, \@qualifiers ]?
743    
744            if ( ref $features eq 'ARRAY' )
745            {
746                $entry->{ ftr_list } = $features;
747            }
748    
749            #  Is it a hash of ( $type => [ [ $location, \@qualifiers ], ... ] )?
750            #  Build a list of features:
751    
752            elsif ( ref $features eq 'HASH' )
753            {
754                my @features;
755                foreach my $type ( keys %$features )
756                {
757                    push @features, map { [ $type, @$_ ] } @{ $features->{ $type } };
758                }
759                $entry->{ ftr_list } = sort_feature_list( \@features );
760            }
761        }
762    
763        if ( $_[0] )
764        {
765            my @features = grep { $_->[0] eq $_[0] } @{ $entry->{ ftr_list } };
766            return wantarray ? @features : \@features;
767        }
768    
769        wantarray ? @{ $entry->{ ftr_list } } : $entry->{ ftr_list };
770    }
771    
772    
773    #-------------------------------------------------------------------------------
774    #  Order features by location.
775    #
776    #      @feautres = sort_feature_list(  @features )
777    #     \@feautres = sort_feature_list( \@features )
778    #
779    #  Handles [ $location, \%quals] and [ $type, $location, \%quals ]
780    #-------------------------------------------------------------------------------
781    sub sort_feature_list
782    {
783        my $by_ref = $_[0]      && ( ref $_[0]      eq 'ARRAY' )
784                  && $_[0]->[0] && ( ref $_[0]->[0] eq 'ARRAY' );
785    
786        my @features = map  { $_->[0] }
787                       sort { $a->[1] <=> $b->[1] || $b->[2] <=> $a->[2] }
788                       map  { [ $_, end_coordinates( $_->[-2] ) ] }
789                       $by_ref ? @{ $_[0] } : @_;
790    
791        wantarray ? @features : \@features;
792    }
793    
794    
795    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
796    #  For a feature location, find its left and right end coordinates.
797    #
798    #      ( $left, $right ) = end_coordinates( $location )
799    #
800    #  Rather than parsing, extract a list of coordinates, and take the extremes.
801    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
802    sub end_coordinates
803    {
804        local $_ = shift;
805        s/([,(])[^,(]+\:[^,)]+([,)])/$1$2/g;     #  Delete locations in other contigs
806        ( sort { $a <=> $b } m/(\d+)/g )[0,-1];  #  Return the lowest and highest values
807    }
808    
809    
810    #-------------------------------------------------------------------------------
811  #  Sequence of a feature.  In list context, include information on partial ends.  #  Sequence of a feature.  In list context, include information on partial ends.
812  #  Can get extract the data from a DNA string, reference to a string, or from  #  Can get extract the data from a DNA string, reference to a string, or from
813  #  the SEQUENCE in an entry.  #  the SEQUENCE in an entry.
# Line 538  Line 819 
819  #   ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna,   $ftr )  #   ( $seq, $partial_5, $partial_3 ) = ftr_dna( \$dna,   $ftr )
820  #   ( $seq, $partial_5, $partial_3 ) = ftr_dna(  $entry, $ftr )  #   ( $seq, $partial_5, $partial_3 ) = ftr_dna(  $entry, $ftr )
821  #  #
822    #  Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
823  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
824  sub ftr_dna  sub ftr_dna
825  {  {
# Line 584  Line 866 
866      wantarray ? ( $seq, $partial_5, $partial_3 ) : $seq;      wantarray ? ( $seq, $partial_5, $partial_3 ) : $seq;
867  }  }
868    
869    
870  sub extract_span  sub extract_span
871  {  {
872      my ( $dnaR, $span ) = @_;      my ( $dnaR, $span ) = @_;
# Line 595  Line 878 
878    
879    
880  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
881    #  Get the location of a feature
882    #
883    #    $ftr_location = location( $ftr )   #  Returns empty string on failure.
884    #
885    #  Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
886    #-------------------------------------------------------------------------------
887    sub location
888    {
889        my ( $ftr ) = @_;
890    
891        ( defined( $ftr )
892             && ( ref( $ftr ) eq 'ARRAY' )
893             && ( @$ftr > 1 ) )
894             ? $ftr->[-2]
895             : '';
896    }
897    
898    
899    #-------------------------------------------------------------------------------
900  #  Identify features with partial 5' or 3' ends.  #  Identify features with partial 5' or 3' ends.
901  #  #
902  #     $partial_5_prime = partial_5_prime( $ftr )  #     $partial_5_prime = partial_5_prime( $ftr )
903  #     $partial_3_prime = partial_3_prime( $ftr )  #     $partial_3_prime = partial_3_prime( $ftr )
904  #  #
905    #  Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
906  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
907  sub partial_5_prime  sub partial_5_prime
908  {  {
# Line 626  Line 929 
929    
930    
931  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
932    #  Get the qualifier hash for a feature:
933  #  #
934  #    $ftr_location = location( $ftr )   #  Returns empty string on failure.  #    \%ftr_qualifiers = qualifiers( $ftr )
 #  
 #-------------------------------------------------------------------------------  
 sub location  
 {  
     my ( $ftr ) = @_;  
   
     ( defined( $ftr )  
          && ( ref( $ftr ) eq 'ARRAY' )  
          && ( @$ftr > 1 ) )  
          ? $ftr->[0]  
          : '';  
 }  
   
   
 #-------------------------------------------------------------------------------  
 #  
 #  \%ftr_qualifiers = qualifiers( $ftr )   #  Returns empty hash reference on failure.  
935  #  #
936    #  Handles both [ $location, \%quals ] and [ $type, $location, \%quals ]
937    #  Returns empty hash reference on failure.
938    #  Note that this list may include  _qualif_order => \@keys
939  #-------------------------------------------------------------------------------  #-------------------------------------------------------------------------------
940  sub qualifiers  sub qualifiers
941  {  {
# Line 654  Line 944 
944      ( defined( $ftr )      ( defined( $ftr )
945           && ( ref( $ftr ) eq 'ARRAY' )           && ( ref( $ftr ) eq 'ARRAY' )
946           && ( @$ftr > 1 )           && ( @$ftr > 1 )
947           && defined( $qual = $ftr->[1] )           && defined( $qual = $ftr->[-1] )
948           && ( ref( $qual ) eq 'HASH' ) )           && ( ref( $qual ) eq 'HASH' ) )
949           ? $qual           ? $qual
950           : {};           : {};
# Line 1187  Line 1477 
1477      return $gb;      return $gb;
1478  }  }
1479    
1480    
1481    
1482    #===============================================================================
1483    #  Write a GenBank entry:
1484    #
1485    #     write_genbank(            \%entry );
1486    #     write_genbank( $filename, \%entry );
1487    #     write_genbank( \*FH,      \%entry );
1488    #
1489    #===============================================================================
1490    #  Recognized record types are:
1491    #
1492    #   + ACCESSION  => [ Accession, ... ]
1493    #     COMMENT    => [ Comment, ... ]
1494    #     DBLINK     => [ Field, ... ]
1495    #     DBSOURCE   => [ Field, ... ]
1496    #     DEFINITION =>   Definition
1497    #     KEYWORDS   => { Key_phrase => 1, Key_phrase => 1, ... }
1498    #   + LOCUS      =>   Locus_id            # locus_id is preferred
1499    #   * ORGANISM   =>   Organism_name
1500    #     ORIGIN     =>   Description_of_sequence_origin
1501    #     REFERENCES => [ { Field => Value, Field => Value, ... }, ... ]
1502    #   * SEQUENCE   =>   Sequence
1503    #     SOURCE     =>   Source_string
1504    #     TAXONOMY   => [ Taxon, Taxon, ... ]
1505    #   + VERSION    => [ Version, Other_information ]
1506    #
1507    #  Data that are not in record types, but supply data to build records:
1508    #
1509    #     date       =>   Date
1510    #     geometry   =>   linear | circular
1511    #     gi         =>   Entry_gi_number
1512    #     is_protein =>   Bool
1513    #     locus_id   =>   Locus_id
1514    #     mol_type   =>   Mol_type
1515    #
1516    #  The record types marked with a * are required.
1517    #  At least of of the types marked with a + must also be present.
1518    #
1519    #  Feature records are merged by type.  Slash is removed from qualifier name.
1520    #  Surrounding quotation marks are stripped from qualifier values.
1521    #
1522    #     FEATURES   => { Type => [ [ Location, { Qualifier => \@values } ],
1523    #                               [ Location, { Qualifier => \@values } ],
1524    #                               ...
1525    #                             ],
1526    #                     Type => ...
1527    #                   }
1528    #
1529    #          1         2         3         4         5         6         7         8
1530    # 12345678901234567890123456789012345678901234567890123456789012345678901234567890
1531    # LOCUS       @<<<<<<<<<<<<<<< @########## @< DNA        @<<<<<<< @<< @<<<<<<<<<<
1532    # LOCUS       CYC_HUMAN                105 aa            linear   PRI 15-MAR-2004
1533    # LOCUS       NC_000909            1664970 bp    DNA     circular BCT 03-DEC-2005
1534    #
1535    #===============================================================================
1536    
1537    
1538    sub write_genbank
1539    {
1540        my ( $fh, $close );
1541        if ( $_[0] && ( ( ! ref( $_[0] ) ) || ( ref( $_[0] ) eq 'GLOB' ) ) )
1542        {
1543            ( $fh, $close ) = output_filehandle( shift );
1544        }
1545    
1546        my $entry = shift;
1547        $entry && ref $entry eq 'HASH' or return 0;
1548        $entry->{ SEQUENCE } or return 0;
1549        $entry->{ ORGANISM } or return 0;
1550    
1551        #  Prepare the data for the forms:
1552    
1553        my $key;
1554    
1555        if ( ! $entry->{ locus_id } )
1556        {
1557            ( $key ) = grep { m/locus_id/i  } keys %$entry;
1558            ( $key ) = grep { m/locus/i     } keys %$entry if ! $key;
1559            ( $key ) = grep { m/contig/i    } keys %$entry if ! $key;
1560            ( $key ) = grep { m/accession/i } keys %$entry if ! $key;
1561            ( $key ) = grep { m/version/i   } keys %$entry if ! $key;
1562            my $locus_id = $key ? $entry->{ $key } : 'undefined';
1563            $locus_id = $locus_id->[0] if ref $locus_id;
1564            $locus_id =~ s/\s.*$//;
1565            $entry->{ locus_id } = $locus_id;
1566        }
1567    
1568        if ( ! $entry->{ mol_type } )
1569        {
1570            my $mol_type;
1571            my $is_prot  = $entry->{ is_protein };
1572            $mol_type = ' ' if $is_prot;
1573            if ( ! $mol_type )
1574            {
1575                $mol_type = mol_type( \$entry->{ SEQUENCE } );
1576                $mol_type = ' ' if $mol_type =~ m/^prot/i;
1577                $entry->{ is_protein } = ( $mol_type =~ m/NA$/i ) ? 0 : 1;
1578            }
1579            $entry->{ mol_type } = $mol_type;
1580        }
1581    
1582        if ( ! $entry->{ DEFINITION } )
1583        {
1584            ( $key ) = grep { m/definition/i } keys %$entry;
1585            ( $key ) = grep { m/^def/i       } keys %$entry if ! $key;
1586            my $def = $key ? $entry->{ $key } : $entry->{ locus_id };
1587            $entry->{ DEFINITION } = $def;
1588        }
1589    
1590        if ( ! $entry->{ ACCESSION } || ! ref $entry->{ ACCESSION } )
1591        {
1592            ( $key ) = grep { m/accession/i } keys %$entry;
1593            ( $key ) = grep { m/^acc/i      } keys %$entry if ! $key;
1594            ( $key ) = grep { m/version/i   } keys %$entry if ! $key;
1595            ( $key ) = grep { m/^ver/i      } keys %$entry if ! $key;
1596              $key   = 'locus_id'                          if ! $key;
1597            my $acc = $entry->{ $key };
1598            $entry->{ ACCESSION } = ref $acc ? $acc : [ $acc ];
1599        }
1600    
1601        if ( ! $entry->{ VERSION } || ! ref $entry->{ VERSION } )
1602        {
1603            ( $key ) = grep { m/version/i } keys %$entry;
1604            ( $key ) = grep { m/^ver/i    } keys %$entry if ! $key;
1605              $key   = 'ACCESSION'                       if ! $key;
1606            my $ver = $entry->{ $key };
1607            $entry->{ VERSION } = ref $ver ? $ver : [ $ver ];
1608        }
1609    
1610        if ( ! grep { m/^gi:/ } @{ $entry->{ VERSION } } )
1611        {
1612            ( $key ) = grep { m/^gi$/i } keys %$entry;
1613            push @{ $entry->{ VERSION } }, "GI:$entry->{$key}" if $key;
1614        }
1615    
1616    
1617        #  Fill out the forms
1618    
1619        my $old_fh;
1620        $old_fh = select( $fh ) if ref $fh eq 'GLOB';
1621    
1622        write_locus( $entry );
1623        write_definition( $entry );
1624        write_accession( $entry );
1625        write_version( $entry );
1626        write_dblink( $entry )       if $entry->{ DBLINK };
1627        write_dbsource( $entry )     if $entry->{ DBSOURCE };
1628        write_keywords( $entry );
1629        write_source( $entry );
1630        write_organism( $entry );
1631        write_taxonomy( $entry );
1632        write_references( $entry );
1633        write_comment( $entry )      if $entry->{ COMMENT };
1634        write_features( $entry );
1635        write_base_count( $entry );
1636        write_origin( $entry );
1637        write_sequence( $entry );
1638        write_end_of_entry();
1639    
1640        select( $old_fh ) if $old_fh;
1641    
1642        close( $fh ) if $fh && $close;
1643    
1644        return 1;
1645    }
1646    
1647    
1648    sub mol_type
1649    {
1650        return undef if ! defined $_[0];
1651        local $_ = ref $_[0] ? shift : \$_[0];
1652        my $n_nt = $$_ =~ tr/ACGNTUacgntu//;
1653        my $n    = length( $$_ );
1654        return ( $n_nt < 0.8 * $n )                ? 'protein'
1655             : ( $$_ =~ tr/Tt// > $$_ =~ tr/Uu// ) ? 'DNA' : 'RNA';
1656    }
1657    
1658    
1659    
1660    sub write_locus
1661    {
1662        my $entry = shift;
1663        my ( $locus, $mol_type, $geometry, $div, $date )
1664            = map { $entry->{ $_ } || ' ' }
1665              qw( locus_id mol_type geometry division date );
1666    
1667        my $length = $entry->{ length } ||= length( $entry->{ SEQUENCE } );
1668    
1669        my $unit = ( $mol_type =~ m/NA$/i ) ? 'bp' : 'aa';
1670    
1671        my $stranded = ' ';
1672        ( $stranded, $mol_type ) = ( $1, $2 ) if $mol_type =~ /^(\S+-)(\S+)$/;
1673    
1674        my $key;
1675        ( $key ) = grep { m/^geometry/i } keys %$entry;
1676        ( $key ) = grep { m/^geom/i     } keys %$entry if ! $key;
1677        $geometry = $entry->{ geometry } ||= $key ? $entry->{ $key }
1678                                              : ' ';
1679    
1680        ( $key ) = grep { m/^div/i } keys %$entry;
1681        $div = $entry->{ division } ||= $key ? $entry->{ $key } : 'UNK';
1682    
1683        ( $key ) = grep { m/^date/i } keys %$entry;
1684        ( $key ) = grep { m/date/i  } keys %$entry if ! $key;
1685        $date = $entry->{ date } ||= $key ? $entry->{ $key }
1686                                             : ( map { chomp; uc } `date '+%d-%b-%Y'` )[0];
1687    
1688        $^A = '';
1689        formline( <<'End_locus', $locus, $length, $unit, $stranded, $mol_type, $geometry, $div, $date );
1690    LOCUS       @<<<<<<<<<<<<<<<<< @######## @< @>>@<<<<<< @<<<<<<< @<< @<<<<<<<<<<
1691    End_locus
1692        print $^A;
1693        $^A = '';
1694    }
1695    
1696    
1697    sub write_definition
1698    {
1699        local $_ = $_[0]->{ DEFINITION } || '';
1700        $_ .= '.' unless /\.$/;
1701        output_record( 'DEFINITION', $_ );
1702    }
1703    
1704    
1705    sub output_record
1706    {
1707        my ( $type, $text ) = @_;
1708    
1709        $^A = '';
1710    
1711        formline( <<'End_of_line_1', $type, $text );
1712    @<<<<<<<<<  ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1713    End_of_line_1
1714    
1715        formline( <<'End_of_continuation', $text ) if defined $text && length $text;
1716    ~~          ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1717    End_of_continuation
1718    
1719        print $^A;
1720        $^A = '';
1721    }
1722    
1723    
1724    sub write_continuation
1725    {
1726        local $_ = shift;
1727    
1728        $^A = '';
1729        formline( <<'End_of_continuation', $_ ) if defined $_ && length $_;
1730    ~~          ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1731    End_of_continuation
1732        print $^A;
1733        $^A = '';
1734    }
1735    
1736    
1737    sub write_accession
1738    {
1739        local $_ = join( ' ', @{ $_[0]->{ ACCESSION } || [] } );
1740        output_record( 'ACCESSION', $_ ) if length $_;
1741    }
1742    
1743    
1744    sub write_version
1745    {
1746        local $_ = join( ' ', @{ $_[0]->{ VERSION } || [] } );
1747        $_ =~ s/ GI:/  GI:/g;
1748        output_record( 'VERSION', $_ ) if length $_;
1749    }
1750    
1751    
1752    sub write_dblink
1753    {
1754        local $_ = join( "\r", @{ $_[0]->{ DBLINK } || [] } );
1755        output_record( 'DBLINK', $_ ) if length $_;
1756    }
1757    
1758    
1759    sub write_dbsource
1760    {
1761        local $_ = join( "\r", @{ $_[0]->{ DBSOURCE } || [] } );
1762        output_record( 'DBSOURCE', $_ ) if length $_;
1763    }
1764    
1765    
1766    sub write_keywords
1767    {
1768        local $_ = join( '; ', @{ $_[0]->{ KEYWORDS } || [] } );
1769        $_ .= '.' if ! /\.$/;
1770        output_record( 'KEYWORDS', $_ );
1771    }
1772    
1773    
1774    sub write_source
1775    {
1776        local $_ = $_[0]->{ SOURCE } || 'unknown';
1777        output_record( 'SOURCE', $_ );
1778    }
1779    
1780    
1781    sub write_organism
1782    {
1783        local $_ = $_[0]->{ ORGANISM } || 'unclasified';
1784        output_record( '  ORGANISM', $_ );
1785    }
1786    
1787    
1788    sub write_taxonomy
1789    {
1790        local $_ = join( '; ', @{ $_[0]->{ TAXONOMY } || [] } );
1791        $_ .= '.' if ! /\.$/;
1792        output_record( ' ', $_ );
1793    }
1794    
1795    
1796    sub write_references
1797    {
1798        my $n = 0;
1799        foreach ( @{ $_[0]->{ REFERENCES } || [] } )
1800        {
1801            write_reference( $_, ++$n );
1802            write_authors( $_ );
1803            write_consortium( $_ );
1804            write_title( $_ );
1805            write_journal( $_ );
1806            write_pubmed( $_ );
1807            write_remark( $_ );
1808        }
1809    }
1810    
1811    
1812    sub write_reference
1813    {
1814        local $_ = $_[0]->{ ref_num } || $_[1] || '1';
1815        $_ .= "  ($_[0]->{ref_note})" if $_[0]->{ ref_note };
1816        output_record( 'REFERENCE', $_ );
1817    }
1818    
1819    
1820    sub write_authors
1821    {
1822        local $_ = $_[0]->{ AUTHORS };
1823        output_record( '  AUTHORS', $_ ) if defined $_ && length $_;
1824    }
1825    
1826    
1827    sub write_consortium
1828    {
1829        local $_ = $_[0]->{ CONSRTM };
1830        output_record( '  CONSRTM', $_ ) if defined $_ && length $_;
1831    }
1832    
1833    
1834    sub write_title
1835    {
1836        local $_ = $_[0]->{ TITLE };
1837        output_record( '  TITLE', $_ ) if defined $_ && length $_;
1838    }
1839    
1840    
1841    sub write_journal
1842    {
1843        local $_ = $_[0]->{ JOURNAL };
1844        output_record( '  JOURNAL', $_ ) if defined $_ && length $_;
1845    }
1846    
1847    
1848    sub write_pubmed
1849    {
1850        local $_ = $_[0]->{ PUBMED };
1851        output_record( '   PUBMED', $_ ) if defined $_ && length $_;
1852    }
1853    
1854    
1855    sub write_remark
1856    {
1857        local $_ = $_[0]->{ REMARK };
1858        output_record( '  REMARK', $_ ) if defined $_ && length $_;
1859    }
1860    
1861    
1862    sub write_comment
1863    {
1864        local $_ = join( "\r", @{ $_[0]->{ COMMENT } || [] } );
1865        output_record( 'COMMENT', $_ ) if length $_;
1866    }
1867    
1868    
1869    sub write_features
1870    {
1871        my $ftr_list = feature_list( $_[0] );
1872        if ( $ftr_list && @$ftr_list )
1873        {
1874            write_feature_header();
1875            foreach ( @$ftr_list ) { write_feature( $_ ) }
1876        }
1877    }
1878    
1879    
1880    sub write_feature_header
1881    {
1882        print "FEATURES             Location/Qualifiers\n";
1883    }
1884    
1885    
1886    sub write_feature
1887    {
1888        my ( $type, $location, $qualifiers ) = @{ $_[0] };
1889    
1890        $^A = '';
1891    
1892        my ( $loc, $rest ) = wrap_location( $location );
1893    
1894        formline( <<'End_of_line_1', $type, $loc );
1895         @<<<<<<<<<<<<<< ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1896    End_of_line_1
1897    
1898    
1899        while ( $rest )
1900        {
1901            ( $loc, $rest ) = wrap_location( $rest );
1902            formline( <<'End_of_continuation', $loc ) if defined $loc && length( $loc );
1903                         ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1904    End_of_continuation
1905        }
1906    
1907        if ( defined $qualifiers && ref $qualifiers eq 'ARRAY' )
1908        {
1909            foreach ( @$qualifiers )
1910            {
1911                add_feature_qualifier( @$_ ) if $_->[0] !~ /^_/;
1912            }
1913        }
1914        elsif ( defined $qualifiers && ref $qualifiers eq 'HASH' )
1915        {
1916            foreach my $key ( ordered_qualifiers( $qualifiers ) )
1917            {
1918                foreach ( @{ $qualifiers->{ $key } } ) { add_feature_qualifier( $key, $_ ) }
1919            }
1920        }
1921    
1922        print $^A;
1923        $^A = '';
1924    }
1925    
1926    
1927    sub wrap_location
1928    {
1929        local $_ = shift;
1930        return ( $_, '' ) if length( $_ ) <= 58;
1931    
1932        return ( $1, $2 ) if m/^(.{1,57},)(.*)$/;
1933        return ( $1, $2 ) if m/^(.{1,56}\.\.)(.*)$/;
1934        return ( $1, $2 ) if m/^(.{1,57}\()(.*)$/;
1935        return $_;
1936    }
1937    
1938    
1939    sub ordered_qualifiers
1940    {
1941        my $qualifiers = shift;
1942    
1943        if ( ! $qualifiers->{ '_qualif_order' } )
1944        {
1945            @{ $qualifiers->{ '_qualif_order' } } = sort { qualifier_priority( $a ) <=> qualifier_priority( $b )
1946                                                        || lc $a cmp lc $b
1947                                                         }
1948                                                    grep { ! /^_/ }
1949                                                    keys %$qualifiers;
1950        }
1951    
1952        @{ $qualifiers->{ '_qualif_order' } };
1953    }
1954    
1955    
1956    sub qualifier_priority { $qualifier_order{ lc $_[0] } || 999999 }
1957    
1958    
1959    sub add_feature_qualifier
1960    {
1961        my ( $key, $value ) = @_;
1962    
1963        my $qualif;
1964        if ( ! defined $value || $value eq '' )
1965        {
1966            $qualif = "/$key";
1967        }
1968        else
1969        {
1970            if ( $value !~ /^\d+$/ ) { $value =~ s/"/""/g; $value = qq("$value"); }
1971            $qualif = "/$key=$value";
1972        }
1973    
1974        formline( <<'End_of_continuation', $qualif ) if defined $qualif && length $qualif;
1975    ~~                   ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1976    End_of_continuation
1977    }
1978    
1979    
1980    sub write_base_count
1981    {
1982        my $entry = shift;
1983        return if $entry->{ is_protein };
1984        return if $entry->{ mol_type } =~ /\S/ && $entry->{ mol_type } !~ m/NA$/i;
1985    
1986        my ( $a, $c, $g, $t, $o );
1987        my $counts = $entry->{ BASE_COUNT };
1988        if ( $counts && ref $counts eq 'HASH' && keys %$counts)
1989        {
1990            ( $a, $c, $g, $t, $o ) = map { $counts->{ $_ } || 0 }
1991                                     qw( a c g t other );
1992            $a || $c || $g || $t || $o or return;
1993        }
1994        else
1995        {
1996            $entry->{ SEQUENCE } or return;
1997            my $s = \$entry->{ SEQUENCE };
1998            $a = $$s =~ tr/Aa//;
1999            $c = $$s =~ tr/Cc//;
2000            $g = $$s =~ tr/Gg//;
2001            $t = $$s =~ tr/Tt//;
2002            $o = length( $$s ) - ( $a + $c + $g + $t );
2003            $a || $c || $g || $t || $o or return;
2004    
2005            $entry->{ BASE_COUNT } = { a => $a,
2006                                       c => $c,
2007                                       g => $g,
2008                                       t => $t
2009                                     };
2010            $entry->{ BASE_COUNT }->{ other } = $o if $o;
2011        }
2012    
2013        printf "BASE COUNT  %6d a %6dc %6d g %6d t%s\n",
2014               $a, $c, $g, $t, $o ? sprintf( ' %6d', $o ) : '';
2015    }
2016    
2017    
2018    sub write_origin
2019    {
2020        my $entry = shift;
2021        my $info = $entry && defined $entry->{ ORIGIN } ? $entry->{ ORIGIN } : '';
2022        print "ORIGIN      $info\n";  # Use explicit print to get blank padding
2023    }
2024    
2025    
2026    sub write_sequence
2027    {
2028        my $entry = shift;
2029        $entry->{ SEQUENCE } or return;
2030    
2031        my $start = 1;
2032        foreach ( $entry->{ SEQUENCE } =~ m/(.{1,60})/g )
2033        {
2034            print join( ' ', sprintf( '%9d', $start ), m/(.{1,10})/g ), "\n";
2035            $start += 60;
2036        }
2037    }
2038    
2039    
2040    sub write_end_of_entry
2041    {
2042        print "//\n";
2043    }
2044    
2045    
2046  #===============================================================================  #===============================================================================
2047  #  Helper function for defining an input filehandle:  #  Helper function for defining an input filehandle:
2048  #  #
2049  #     filehandle is passed through  #     filehandle is passed through
2050  #     string is taken as file name to be openend  #     string is taken as file name to be opened
2051  #     undef or "" defaults to STDOUT  #     undef or "" defaults to STDIN
2052  #  #
2053  #      \*FH           = input_filehandle( $file );  #      \*FH           = input_filehandle( $file );
2054  #    ( \*FH, $close ) = input_filehandle( $file );  #    ( \*FH, $close ) = input_filehandle( $file );
# Line 1220  Line 2076 
2076    
2077      if ( ! ref( $file ) )      if ( ! ref( $file ) )
2078      {      {
2079          -f $file or die "Could not find input file \"$file\"\n";          -f $file or die "Could not find input file '$file'.\n";
2080          my $fh;          my $fh;
2081          open( $fh, "<$file" ) || die "Could not open \"$file\" for input\n";          open( $fh, "<$file" ) || die "Could not open '$file' for input.\n";
2082          return wantarray ? ( $fh, 1 ) : $fh;          return wantarray ? ( $fh, 1 ) : $fh;
2083      }      }
2084    
# Line 1230  Line 2086 
2086  }  }
2087    
2088    
2089    #===============================================================================
2090    #  Helper function for defining an output filehandle:
2091    #
2092    #     filehandle is passed through
2093    #     string is taken as file name to be opened
2094    #     undef or "" defaults to STDOUT
2095    #
2096    #      \*FH           = output_filehandle( $file );
2097    #    ( \*FH, $close ) = output_filehandle( $file );
2098    #
2099    #===============================================================================
2100    sub output_filehandle
2101    {
2102        my $file = shift;
2103    
2104        #  Null string or undef
2105    
2106        if ( ! defined( $file ) || ( $file eq '' ) )
2107        {
2108            return wantarray ? ( \*STDOUT, 0 ) : \*STDOUT;
2109        }
2110    
2111        #  FILEHANDLE?
2112    
2113        if ( ref( $file ) eq "GLOB" )
2114        {
2115            return wantarray ? ( $file, 0 ) : $file;
2116        }
2117    
2118        #  File name
2119    
2120        if ( ! ref( $file ) )
2121        {
2122            my $fh;
2123            open( $fh, ">$file" ) || die "Could not open '$file' for output.\n";
2124            return wantarray ? ( $fh, 1 ) : $fh;
2125        }
2126    
2127        return wantarray ? ( \*STDOUT, undef ) : \*STDOUT;
2128    }
2129    
2130    
2131  1;  1;

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3