← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/bioperl/oneiric

« back to all changes in this revision

Viewing changes to Bio/SeqIO/genbank.pm

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/SeqIO/genbank.pm
1
 
# $Id: genbank.pm,v 1.50.2.1 2002/03/13 16:48:28 jason Exp $
 
1
# $Id: genbank.pm,v 1.99 2003/12/22 18:33:15 heikki Exp $
2
2
#
3
3
# BioPerl module for Bio::SeqIO::GenBank
4
4
#
5
 
# Cared for by Elia Stupka <elia@ebi.ac.uk>
 
5
# Cared for by Elia Stupka <elia@tll.org.sg>
6
6
#
7
7
# Copyright Elia Stupka
8
8
#
25
25
        # do something with $seq
26
26
    }
27
27
 
 
28
 
28
29
=head1 DESCRIPTION
29
30
 
30
31
This object can transform Bio::Seq objects to and from GenBank flat
33
34
There is alot of flexibility here about how to dump things which I need
34
35
to document fully.
35
36
 
 
37
=head2 Mapping of record properties to object properties
 
38
 
 
39
This section is supposed to document which sections and properties of
 
40
a GenBank databank record end up where in the Bioperl object model. It
 
41
is far from complete and presently focuses only on those mappings
 
42
which may be non-obvious. $seq in the text refers to the
 
43
Bio::Seq::RichSeqI implementing object returned by the parser for each
 
44
record.
 
45
 
 
46
=over 4
 
47
 
 
48
=item GI number
 
49
 
 
50
$seq-E<gt>primary_id
 
51
 
 
52
=back
36
53
 
37
54
=head2 Optional functions
38
55
 
56
73
To generate the ID line. If it is not there, it generates a sensible ID
57
74
line using a number of tools.
58
75
 
 
76
 
 
77
If you want to output annotations in genbank format they need to be
 
78
stored in a Bio::Annotation::Collection object which is accessible
 
79
through the Bio::SeqI interface method L<annotation()|annotation>.  
 
80
 
 
81
The following are the names of the keys which are polled from a
 
82
L<Bio::Annotation::Collection> object.
 
83
 
 
84
reference       - Should contain Bio::Annotation::Reference objects
 
85
comment         - Should contain Bio::Annotation::Comment objects
 
86
 
 
87
segment         - Should contain a Bio::Annotation::SimpleValue object
 
88
origin          - Should contain a Bio::Annotation::SimpleValue object
 
89
 
59
90
=back
60
91
 
 
92
=head1 Where does the data go?
 
93
 
 
94
Data parsed in Bio::SeqIO::genbank is stored in a variety of data
 
95
fields in the sequence object that is returned.  More information in
 
96
the HOWTOs about exactly what each field means and where it goes.
 
97
Here is a partial list of fields.
 
98
 
 
99
Items listed as RichSeq or Seq or PrimarySeq and then NAME() tell you
 
100
the top level object which defines a function called NAME() which
 
101
stores this information.
 
102
 
 
103
Items listed as Annotation 'NAME' tell you the data is stored the
 
104
associated Bio::Annotation::Colection object which is associated with
 
105
Bio::Seq objects.  If it is explictly requested that no annotations
 
106
should be stored when parsing a record of course they won't be
 
107
available when you try and get them.  If you are having this problem
 
108
look at the type of SeqBuilder that is being used to contruct your
 
109
sequence object.
 
110
 
 
111
Comments             Annotation 'comment'
 
112
References           Annotation 'reference'
 
113
Segment              Annotation 'segment'
 
114
Origin               Annotation 'origin'
 
115
 
 
116
Accessions           PrimarySeq accession_number()
 
117
Secondary accessions RichSeq get_secondary_accessions()
 
118
Keywords             RichSeq keywords()
 
119
Dates                RichSeq get_dates()
 
120
Molecule             RichSeq molecule()
 
121
Seq Version          RichSeq seq_version()
 
122
PID                  RichSeq pid()
 
123
Division             RichSeq division()
 
124
Features             Seq get_SeqFeatures()
 
125
Alphabet             PrimarySeq alphabet()
 
126
Definition           PrimarySeq description() or desc()
 
127
Version              PrimarySeq version()
 
128
 
 
129
Sequence             PrimarySeq seq()
 
130
 
61
131
=head1 FEEDBACK
62
132
 
63
133
=head2 Mailing Lists
77
147
 Bug reports can be submitted via email or the web:
78
148
 
79
149
  bioperl-bugs@bio.perl.org
80
 
  http://bio.perl.org/bioperl-bugs/
 
150
  http://bugzilla.bioperl.org/
81
151
 
82
152
=head1 AUTHOR - Elia Stupka
83
153
 
84
 
Email elia@ebi.ac.uk
 
154
Email elia@tll.org.sg
85
155
 
86
156
=head1 CONTRIBUTORS
87
157
 
90
160
Chris Mungall cjm@fruitfly.bdgp.berkeley.edu
91
161
Lincoln Stein lstein@cshl.org
92
162
Heikki Lehvaslaiho, heikki@ebi.ac.uk
 
163
Hilmar Lapp, hlapp@gmx.net
 
164
Donald G. Jackson, donald.jackson@bms.com
93
165
 
94
166
=head1 APPENDIX
95
167
 
96
 
The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _
 
168
The rest of the documentation details each of the object
 
169
methods. Internal methods are usually preceded with a _
97
170
 
98
171
=cut
99
172
 
100
173
# Let the code begin...
101
174
 
102
175
package Bio::SeqIO::genbank;
103
 
use vars qw(@ISA);
 
176
use vars qw(@ISA %FTQUAL_NO_QUOTE);
104
177
use strict;
105
178
 
106
 
use Bio::Seq::RichSeq;
107
179
use Bio::SeqIO;
108
180
use Bio::SeqIO::FTHelper;
109
181
use Bio::SeqFeature::Generic;
110
182
use Bio::Species;
 
183
use Bio::Seq::SeqFactory;
 
184
use Bio::Annotation::Collection;
 
185
use Bio::Annotation::Comment;
 
186
use Bio::Annotation::Reference;
 
187
use Bio::Annotation::DBLink;
111
188
 
112
189
@ISA = qw(Bio::SeqIO);
113
 
 
 
190
 
 
191
%FTQUAL_NO_QUOTE=(
 
192
                  'anticodon'    => 1,
 
193
                  'citation'     => 1,
 
194
                  'codon'        => 1,
 
195
                  'codon_start'  => 1,
 
196
                  'cons_splice'  => 1,
 
197
                  'direction'    => 1,
 
198
                  'evidence'     => 1,
 
199
                  'label'        => 1,
 
200
                  'mod_base'     => 1,
 
201
                  'number'       => 1,
 
202
                  'rpt_type'     => 1,
 
203
                  'rpt_unit'     => 1,
 
204
                  'transl_except'=> 1,
 
205
                  'transl_table' => 1,
 
206
                  'usedin'       => 1,
 
207
                  );
 
208
 
114
209
sub _initialize {
115
210
    my($self,@args) = @_;
116
 
    $self->SUPER::_initialize(@args);
117
 
 
 
211
    
 
212
    $self->SUPER::_initialize(@args); 
118
213
    # hash for functions for decoding keys.
119
214
    $self->{'_func_ftunit_hash'} = {}; 
120
215
    $self->_show_dna(1); # sets this to one by default. People can change it
 
216
    if( ! defined $self->sequence_factory ) {
 
217
        $self->sequence_factory(new Bio::Seq::SeqFactory
 
218
                                (-verbose => $self->verbose(), 
 
219
                                 -type => 'Bio::Seq::RichSeq'));
 
220
    }
121
221
}
122
222
 
123
223
=head2 next_seq
132
232
 
133
233
sub next_seq {
134
234
    my ($self,@args) = @_;
135
 
    my ($pseq,$c,$name,@desc,$seqc,$mol,$div,$date);
136
 
    my $line;
137
 
    my @acc = ();
138
 
    my $seq = Bio::Seq::RichSeq->new('-verbose' =>$self->verbose());
139
 
    
140
 
    while(defined($line = $self->_readline())) {
141
 
        $line =~ /^LOCUS\s+\S+/ && last;
142
 
    }
143
 
    return undef if( !defined $line ); # end of file
144
 
    
145
 
    $line =~ /^LOCUS\s+\S+/ ||
146
 
        $self->throw("GenBank stream with no LOCUS. Not GenBank in my book. Got $line");
147
 
    
148
 
    $line =~ /^LOCUS\s+(\S+)\s+\S+\s+(bp|aa)\s+(\S+)\s+(\S+)\s+(\S+)?/ || do {
149
 
        $line =~ /^LOCUS\s+(\S+)/ ||
150
 
            $self->throw("GenBank stream with no LOCUS. Not GenBank in my book. Got $line");
151
 
        # fall back to at least an ID
152
 
        $name = $1;
153
 
        $self->warn("GenBank record with LOCUS line in unexpected format. ".
154
 
                    "Attributes from this line other than ID will be missing.");
155
 
    };
156
 
    $name = $1;
157
 
    # this is important to have the id for display in e.g. FTHelper, otherwise
158
 
    # you won't know which entry caused an error
159
 
    $seq->display_id($name);
160
 
    # the alphabet of the entry
161
 
    if($2 eq 'bp') {
162
 
        $seq->alphabet('dna');
163
 
    } else {
164
 
        # $2 eq 'aa'
165
 
        $seq->alphabet('protein');
166
 
    }
167
 
    # for aa there is usually no 'molecule' (mRNA etc)
168
 
    if (($2 eq 'bp') || defined($5)) {
169
 
        if ($4 eq 'circular') {
170
 
            $seq->molecule($3);
171
 
            $seq->is_circular($4);
172
 
            $seq->division($5);
173
 
        } else {
174
 
            $seq->molecule($3);
175
 
            if (CORE::length($4) == 3 ) { 
176
 
                $seq->division($4);
177
 
            } else {
178
 
                $seq->division($5);
179
 
            }
180
 
        }
181
 
    } else {
182
 
        $seq->molecule('PRT') if($2 eq 'aa');
183
 
        $seq->division($3);
184
 
        $date = $4;
185
 
    }
186
 
    ($date) = ($line =~ /.*(\d\d-\w\w\w-\d\d\d\d)/);
187
 
    if ($date) {
188
 
        $seq->add_date($date);
189
 
    }
190
 
    my $buffer = $self->_readline();
191
 
        
192
 
    until( !defined ($buffer) ) {
193
 
        $_ = $buffer;
194
 
        
195
 
        # Description line(s)
196
 
        if (/^DEFINITION\s+(\S.*\S)/) {
197
 
            push (@desc, $1);
198
 
            while ( defined($_ = $self->_readline) ) { 
199
 
                /^\s+(.*)/ && do { push (@desc, $1); next;};
200
 
                last;
201
 
            }
202
 
        }                
203
 
        # accession number (there can be multiple accessions)
204
 
        if( /^ACCESSION\s+(\S.*\S)/ ) {
205
 
            push(@acc, split(' ',$1));
206
 
        }
207
 
        
208
 
        # PID
209
 
        if( /^PID\s+(\S+)/ ) {
210
 
            $seq->pid($1);
211
 
        }
212
 
 
213
 
        #Version number
214
 
        if( /^VERSION\s+(\S+)\.(\d+)\s*(GI:\d+)?/ ) {
215
 
            $seq->seq_version($2);
216
 
            $seq->primary_id(substr($3, 3)) if($3);
217
 
        }
218
 
 
219
 
        #Keywords
220
 
        if( /^KEYWORDS\s+(.*)/ ) {
221
 
            my $keywords = $1;
222
 
            $keywords =~ s/\;//g;
223
 
            $seq->keywords($keywords);
224
 
        }
225
 
        
226
 
        # Organism name and phylogenetic information
227
 
        if (/^SOURCE/) {
228
 
            my $species = $self->_read_GenBank_Species(\$buffer);
229
 
            $seq->species( $species );
230
 
            next;
231
 
        }
232
 
        
233
 
        #References
234
 
        if (/^REFERENCE/) {
235
 
            my @refs = $self->_read_GenBank_References(\$buffer);
236
 
            foreach my $ref ( @refs ) {
237
 
                $seq->annotation->add_Annotation('reference',$ref);
238
 
            }
239
 
            next;
240
 
        }
241
 
        
242
 
        #Comments
243
 
        if (/^COMMENT\s+(.*)/) {
244
 
            my $comment = $1;
245
 
            while (defined($_ = $self->_readline)) {
246
 
                if (/^\S/) {
247
 
                    $comment =~ s/\n/ /g;
248
 
                    $comment =~ s/  +/ /g;
249
 
                    my $commobj = Bio::Annotation::Comment->new();
250
 
                    $commobj->text($comment);
251
 
                    $seq->annotation->add_Annotation('comment',$commobj);
252
 
                    last;
253
 
                }
254
 
                $comment .= $_; 
255
 
            }
256
 
            $buffer = $_;
257
 
            next;
258
 
        }
259
 
        # Exit at start of Feature table, or start of sequence
260
 
        last if( /^(FEATURES)|(ORIGIN)/ );
261
 
        # Get next line and loop again
262
 
        $buffer = $self->_readline;
263
 
    }
264
 
    return undef if(! defined($buffer));
265
 
    $seq->desc(join(' ', @desc)); # don't forget!
266
 
    $seq->accession_number(shift(@acc));
267
 
    $seq->add_secondary_accession(@acc);
268
 
    
269
 
    # some "minimal" formats may not necessarily have a feature table
270
 
    if(defined($_) && /^FEATURES/) {
271
 
        # need to read the first line of the feature table
272
 
        $buffer = $self->_readline;
273
 
 
274
 
        # DO NOT read lines in the while condition -- this is done as a side
275
 
        # effect in _read_FTHelper_GenBank!
276
 
        while( defined($buffer) ) {
277
 
            # check immediately -- not at the end of the loop
278
 
            # note: GenPept entries obviously do not have a BASE line
279
 
            last if(($buffer =~ /^BASE/) || ($buffer =~ /^ORIGIN/) ||
280
 
                    ($buffer =~ /^CONTIG/) );
281
 
 
282
 
            # slurp in one feature at a time -- at return, the start of
283
 
            # the next feature will have been read already, so we need
284
 
            # to pass a reference, and the called method must set this
285
 
            # to the last line read before returning 
286
 
 
287
 
            my $ftunit = $self->_read_FTHelper_GenBank(\$buffer);
288
 
            
289
 
            # fix suggested by James Diggans
290
 
 
291
 
            if( !defined $ftunit ) {
292
 
                # GRRRR. We have fallen over. Try to recover
293
 
                $self->warn("Unexpected error in feature table for ".$seq->id." Skipping feature, attempting to recover");
294
 
                unless( ($buffer =~ /^\s{5,5}\S+/) or ($buffer =~ /^\S+/)) {
295
 
                    $buffer = $self->_readline();
296
 
                }
297
 
                next; # back to reading FTHelpers
298
 
            }
 
235
    my $builder = $self->sequence_builder();
 
236
    my $seq;
 
237
    my %params;
 
238
 
 
239
  RECORDSTART: while (1) {
 
240
      my $buffer;
 
241
      my (@acc, @features);
 
242
      my ($display_id, $annotation);
 
243
      my $species;
 
244
 
 
245
      # initialize; we may come here because of starting over
 
246
      @features = ();
 
247
      $annotation = undef;
 
248
      @acc = ();
 
249
      $species = undef;
 
250
      %params = (-verbose => $self->verbose); # reset hash
 
251
      local($/) = "\n";
 
252
      while(defined($buffer = $self->_readline())) {
 
253
          last if index($buffer,'LOCUS       ') == 0;
 
254
      }
 
255
      return undef if( !defined $buffer ); # end of file
 
256
      $buffer =~ /^LOCUS\s+(\S.*)$/o ||
 
257
          $self->throw("GenBank stream with bad LOCUS line. Not GenBank in my book. Got '$buffer'");
 
258
      
 
259
      my @tokens = split(' ', $1);
 
260
 
 
261
      # this is important to have the id for display in e.g. FTHelper,
 
262
      # otherwise you won't know which entry caused an error
 
263
      $display_id = shift(@tokens);
 
264
      $params{'-display_id'} = $display_id;
 
265
      # may still be useful if we don't want the seq
 
266
      $params{'-length'} = shift(@tokens);
 
267
      # the alphabet of the entry
 
268
      $params{'-alphabet'} = (lc(shift @tokens) eq 'bp') ? 'dna' : 'protein';
 
269
      # for aa there is usually no 'molecule' (mRNA etc)
 
270
      if (($params{'-alphabet'} eq 'dna') || (@tokens > 2)) {   
 
271
          $params{'-molecule'} = shift(@tokens);
 
272
          my $circ = shift(@tokens);
 
273
          if ($circ eq 'circular') {
 
274
              $params{'-is_circular'} = 1;
 
275
              $params{'-division'} = shift(@tokens);
 
276
          } else {
 
277
              # 'linear' or 'circular' may actually be omitted altogether
 
278
              $params{'-division'} =
 
279
                  (CORE::length($circ) == 3 ) ? $circ : shift(@tokens);
 
280
          }
 
281
      } else {
 
282
          $params{'-molecule'} = 'PRT' if($params{'-alphabet'} eq 'aa');
 
283
          $params{'-division'} = shift(@tokens);
 
284
      }
 
285
      my $date = join(' ', @tokens); # we lump together the rest
 
286
 
 
287
      # this is per request bug #1513
 
288
      # we can handle
 
289
      # 9-10-2003
 
290
      # 9-10-03
 
291
      #09-10-2003
 
292
      #09-10-03
 
293
      if($date =~ s/\s*((\d{1,2})-(\w{3})-(\d{2,4})).*/$1/) {
 
294
          if( length($date) < 11 ) { # improperly formatted date
 
295
                                     # But we'll be nice and fix it for them
 
296
              my ($d,$m,$y) = ($2,$3,$4);
 
297
              if( length($d) == 1 ) {
 
298
                  $d = "0$d";
 
299
              }
 
300
              # guess the century here
 
301
              if( length($y) == 2 ) {
 
302
                  if( $y > 60 ) {  # arbitrarily guess that '60' means 1960
 
303
                      $y = "19$y";
 
304
                  } else { 
 
305
                      $y = "20$y";
 
306
                  }
 
307
                  $self->warn("Date was malformed, guessing the century for $date to be $y\n");
 
308
              }
 
309
              $params{'-dates'} = [join('-',$d,$m,$y)];
 
310
          } else { 
 
311
              $params{'-dates'} = [$date];
 
312
          }
 
313
      }
 
314
      # set them all at once
 
315
      $builder->add_slot_value(%params);
 
316
      %params = ();
 
317
 
 
318
      # parse the rest if desired, otherwise start over
 
319
      if(! $builder->want_object()) {
 
320
          $builder->make_object();
 
321
          next RECORDSTART;
 
322
      }
 
323
      
 
324
      # set up annotation depending on what the builder wants
 
325
      if($builder->want_slot('annotation')) {
 
326
          $annotation = new Bio::Annotation::Collection;
 
327
      }
 
328
      $buffer = $self->_readline();
 
329
      until( !defined ($buffer) ) {
 
330
          $_ = $buffer;
 
331
          
 
332
          # Description line(s)
 
333
          if (/^DEFINITION\s+(\S.*\S)/) {
 
334
              my @desc = ($1);
 
335
              while ( defined($_ = $self->_readline) ) { 
 
336
                  if( /^\s+(.*)/ ) { push (@desc, $1); next };            
 
337
                  last;
 
338
              }
 
339
              $builder->add_slot_value(-desc => join(' ', @desc));
 
340
              # we'll continue right here because DEFINITION always comes
 
341
              # at the top of the entry
 
342
          }
 
343
          # accession number (there can be multiple accessions)
 
344
          if( /^ACCESSION\s+(\S.*\S)/ ) {
 
345
              push(@acc, split(/\s+/,$1));
 
346
              while( defined($_ = $self->_readline) ) { 
 
347
                  /^\s+(.*)/ && do { push (@acc, split(/\s+/,$1)); next };
 
348
                  last;
 
349
              }
 
350
              $buffer = $_;
 
351
              next;
 
352
          }
 
353
          # PID
 
354
          elsif( /^PID\s+(\S+)/ ) {
 
355
              $params{'-pid'} = $1;
 
356
          }
 
357
          #Version number
 
358
          elsif( /^VERSION\s+(.+)$/ ) {
 
359
              my ($acc,$gi) = split(' ',$1);
 
360
              if($acc =~ /^\w+\.(\d+)/) {
 
361
                  $params{'-version'} = $1;
 
362
                  $params{'-seq_version'} = $1;
 
363
              }
 
364
              if($gi && (index($gi,"GI:") == 0)) {
 
365
                  $params{'-primary_id'} = substr($gi,3);
 
366
              }
 
367
          }
 
368
          #Keywords
 
369
          elsif( /^KEYWORDS\s+(.*)/ ) {
 
370
              my @kw = split(/\s*\;\s*/,$1);
 
371
              while( defined($_ = $self->_readline) ) { 
 
372
                  chomp;
 
373
                  /^\s+(.*)/ && do { push (@kw, split(/\s*\;\s*/,$1)); next };
 
374
                  last;
 
375
              }
 
376
              
 
377
              @kw && $kw[-1] =~ s/\.$//;
 
378
              $params{'-keywords'} = \@kw;
 
379
              $buffer = $_;
 
380
              next;
 
381
          }
 
382
          # Organism name and phylogenetic information
 
383
          elsif (/^SOURCE/) {
 
384
              if($builder->want_slot('species')) {
 
385
                  $species = $self->_read_GenBank_Species(\$buffer);
 
386
                  $builder->add_slot_value(-species => $species);
 
387
              } else {
 
388
                  while(defined($buffer = $self->_readline())) {
 
389
                      last if substr($buffer,0,1) ne ' ';
 
390
                  }
 
391
              }
 
392
              next;
 
393
          }
 
394
          #References
 
395
          elsif (/^REFERENCE/) {
 
396
              if($annotation) {
 
397
                  my @refs = $self->_read_GenBank_References(\$buffer);
 
398
                  foreach my $ref ( @refs ) {
 
399
                      $annotation->add_Annotation('reference',$ref);
 
400
                  }
 
401
              } else {
 
402
                  while(defined($buffer = $self->_readline())) {
 
403
                      last if substr($buffer,0,1) ne ' ';
 
404
                  }
 
405
              }
 
406
              next;
 
407
          }
 
408
          #Comments
 
409
          elsif (/^COMMENT\s+(.*)/) {
 
410
              if($annotation) {
 
411
                  my $comment = $1;
 
412
                  while (defined($_ = $self->_readline)) {
 
413
                      last if (/^\S/);
 
414
                      $comment .= $_; 
 
415
                  }
 
416
                  $comment =~ s/\n/ /g;
 
417
                  $comment =~ s/  +/ /g;
 
418
                  $annotation->add_Annotation(
 
419
                            'comment',
 
420
                            Bio::Annotation::Comment->new(-text => $comment));
 
421
                  $buffer = $_;
 
422
              } else {
 
423
                  while(defined($buffer = $self->_readline())) {
 
424
                      last if substr($buffer,0,1) ne ' ';
 
425
                  }
 
426
              }
 
427
              next;
 
428
          } elsif( /^SEGMENT\s+(.+)/ ) {
 
429
              if($annotation) {
 
430
                  my $segment = $1;
 
431
                  while (defined($_ = $self->_readline)) {
 
432
                      last if (/^\S/);
 
433
                      $segment .= $_; 
 
434
                  }
 
435
                  $segment =~ s/\n/ /og;
 
436
                  $segment =~ s/ {2,}/ /og;
 
437
                  $annotation->add_Annotation(
 
438
                                              'segment',
 
439
                                              Bio::Annotation::SimpleValue->new(-value => $segment));
 
440
                  $buffer = $_;
 
441
              } else {
 
442
                  while(defined($buffer = $self->_readline())) {
 
443
                      last if substr($buffer,0,1) ne ' ';
 
444
                  }
 
445
              }
 
446
              next;
 
447
          }
 
448
          # Exit at start of Feature table, or start of sequence
 
449
          last if( /^(FEATURES|ORIGIN)/ );
 
450
          # Get next line and loop again
 
451
          $buffer = $self->_readline;
 
452
      }
 
453
      return undef if(! defined($buffer));
 
454
 
 
455
      # add them all at once for efficiency
 
456
      $builder->add_slot_value(-accession_number => shift(@acc),
 
457
                               -secondary_accessions => \@acc,
 
458
                               %params);
 
459
      $builder->add_slot_value(-annotation => $annotation) if $annotation;
 
460
      %params = (); # reset before possible re-use to avoid setting twice
 
461
 
 
462
      # start over if we don't want to continue with this entry
 
463
      if(! $builder->want_object()) {
 
464
          $builder->make_object();
 
465
          next RECORDSTART;
 
466
      }      
 
467
      
 
468
      # some "minimal" formats may not necessarily have a feature table
 
469
      if($builder->want_slot('features') && defined($_) && /^FEATURES/o) {
 
470
          # need to read the first line of the feature table
 
471
          $buffer = $self->_readline;
 
472
          
 
473
          # DO NOT read lines in the while condition -- this is done as a side
 
474
          # effect in _read_FTHelper_GenBank!
 
475
          while( defined($buffer) ) {
 
476
              # check immediately -- not at the end of the loop
 
477
              # note: GenPept entries obviously do not have a BASE line
 
478
              last if(($buffer =~ /^BASE/o) || ($buffer =~ /^ORIGIN/o) ||
 
479
                      ($buffer =~ /^CONTIG/o) );
 
480
 
 
481
              # slurp in one feature at a time -- at return, the start of
 
482
              # the next feature will have been read already, so we need
 
483
              # to pass a reference, and the called method must set this
 
484
              # to the last line read before returning 
 
485
              
 
486
              my $ftunit = $self->_read_FTHelper_GenBank(\$buffer);
 
487
              
 
488
              # fix suggested by James Diggans
 
489
 
 
490
              if( !defined $ftunit ) {
 
491
                  # GRRRR. We have fallen over. Try to recover
 
492
                  $self->warn("Unexpected error in feature table for ".$params{'-display_id'}." Skipping feature, attempting to recover");
 
493
                  unless( ($buffer =~ /^\s{5,5}\S+/o) or 
 
494
                          ($buffer =~ /^\S+/o)) {
 
495
                      $buffer = $self->_readline();
 
496
                  }
 
497
                  next; # back to reading FTHelpers
 
498
              }
299
499
                
300
 
            # process ftunit
301
 
            $ftunit->_generic_seqfeature($seq);
302
 
        }
303
 
        $_ = $buffer;
304
 
    }
305
 
    if( defined ($_) ) {
306
 
        if( /^CONTIG/ ) {
307
 
            $b = "     $_"; # need 5 spaces to treat it like a feature
308
 
            my $ftunit = $self->_read_FTHelper_GenBank(\$b);
309
 
            if( ! defined $ftunit ) {
310
 
                $self->warn("unable to parse the CONTIG feature\n");
311
 
            } else { 
312
 
                $ftunit->_generic_seqfeature($seq);
313
 
            }   
314
 
        } elsif(! /^ORIGIN/) {     # advance to the section with the sequence
315
 
            $seqc = ""; 
316
 
            while (defined( $_ = $self->_readline) ) {
317
 
                last if /^ORIGIN/;
318
 
            }
319
 
        }
320
 
    }
321
 
    while( defined($_ = $self->_readline) ) {
322
 
        /^\/\// && last;
323
 
        $_ = uc($_);
324
 
        s/[^A-Za-z]//g;
325
 
        $seqc .= $_;
326
 
    }
327
 
    $seq->seq($seqc);
 
500
              # process ftunit
 
501
              my $feat =
 
502
                  $ftunit->_generic_seqfeature($self->location_factory(),
 
503
                                               $display_id);
 
504
              # add taxon_id from source if available
 
505
              if($species && ($feat->primary_tag eq 'source') &&
 
506
                 $feat->has_tag('db_xref') && (! $species->ncbi_taxid())) {
 
507
                  foreach my $tagval ($feat->get_tag_values('db_xref')) {
 
508
                      if(index($tagval,"taxon:") == 0) {
 
509
                          $species->ncbi_taxid(substr($tagval,6));
 
510
                      }
 
511
                  }
 
512
              }
 
513
              # add feature to list of features
 
514
              push(@features, $feat);
 
515
          }
 
516
          $builder->add_slot_value(-features => \@features);
 
517
          $_ = $buffer;
 
518
      }
 
519
      if( defined ($_) ) {
 
520
          if( /^CONTIG/o && $builder->want_slot('features')) {
 
521
              $b = "     $_"; # need 5 spaces to treat it like a feature
 
522
              my $ftunit = $self->_read_FTHelper_GenBank(\$b);
 
523
              if( ! defined $ftunit ) {
 
524
                  $self->warn("unable to parse the CONTIG feature\n");
 
525
              } else { 
 
526
                  push(@features,
 
527
                       $ftunit->_generic_seqfeature($self->location_factory(),
 
528
                                                    $display_id));
 
529
              } 
 
530
          } elsif(! /^(ORIGIN|\/\/)/ ) {    # advance to the sequence, if any
 
531
              while (defined( $_ = $self->_readline) ) {
 
532
                  last if /^(ORIGIN|\/\/)/;
 
533
              }
 
534
          }
 
535
      }
 
536
      if(! $builder->want_object()) {
 
537
          $builder->make_object(); # implicit end-of-object
 
538
          next RECORDSTART;
 
539
      }
 
540
      if($builder->want_slot('seq')) {
 
541
          # the fact that we want a sequence does not necessarily mean that
 
542
          # there also is a sequence ...
 
543
          if(defined($_) && s/^ORIGIN//) {
 
544
              chomp;
 
545
              if( $annotation && length($_) > 0 ) {
 
546
                  $annotation->add_Annotation('origin',
 
547
                                              Bio::Annotation::SimpleValue->new(-value => $_));
 
548
              }
 
549
              my $seqc = '';
 
550
              while( defined($_ = $self->_readline) ) {
 
551
                  /^\/\// && last;
 
552
                  $_ = uc($_);
 
553
                  s/[^A-Za-z]//g;
 
554
                  $seqc .= $_;
 
555
              }
 
556
              $self->debug("sequence length is ". length($seqc) ."\n");
 
557
              $builder->add_slot_value(-seq => $seqc);
 
558
          }
 
559
      } elsif ( defined($_) && (substr($_,0,2) ne '//')) {
 
560
          # advance to the end of the record
 
561
          while( defined($_ = $self->_readline) ) {
 
562
              last if substr($_,0,2) eq '//';
 
563
          }
 
564
      }
 
565
      # Unlikely, but maybe the sequence is so weird that we don't want it
 
566
      # anymore. We don't want to return undef if the stream's not exhausted
 
567
      # yet.
 
568
      $seq = $builder->make_object();
 
569
      next RECORDSTART unless $seq;
 
570
      last RECORDSTART;
 
571
  } # end while RECORDSTART
328
572
 
329
573
    return $seq;
330
574
}
335
579
 Usage   : $stream->write_seq($seq)
336
580
 Function: writes the $seq object (must be seq) to the stream
337
581
 Returns : 1 for success and 0 for error
338
 
 Args    : Bio::Seq
 
582
 Args    : array of 1 to n Bio::SeqI objects
339
583
 
340
584
 
341
585
=cut
342
586
 
343
587
sub write_seq {
344
 
    my ($self,$seq) = @_;
345
 
    
346
 
    if( !defined $seq ) {
347
 
        $self->throw("Attempting to write with no seq!");
348
 
    }
349
 
    
350
 
    if( ! ref $seq || ! $seq->isa('Bio::SeqI') ) {
351
 
        $self->warn(" $seq is not a SeqI compliant module. Attempting to dump, but may fail!");
352
 
    }
353
 
    
354
 
    my $i;
355
 
    my $str = $seq->seq;
356
 
    
357
 
    my ($div, $mol);
358
 
    my $len = $seq->length();
359
 
    
360
 
 
361
 
    if ( $seq->can('division') ) {
362
 
        $div=$seq->division;
363
 
    } 
364
 
    if( !defined $div || ! $div ) { $div = 'UNK'; }
365
 
 
366
 
    if( !$seq->can('molecule') || ! defined ($mol = $seq->molecule()) ) {
367
 
        $mol = $seq->alphabet || 'DNA';
368
 
    }
369
 
    
370
 
    my $circular = 'linear  ';
371
 
    $circular = 'circular' if $seq->is_circular;
372
 
 
373
 
    local($^W) = 0;   # supressing warnings about uninitialized fields.
374
 
    
375
 
    my $temp_line;
376
 
    if( $self->_id_generation_func ) {
377
 
        $temp_line = &{$self->_id_generation_func}($seq);
378
 
    } else {
379
 
        my $date = '';
380
 
        if( $seq->can('get_dates') ) {      
381
 
            ($date) = $seq->get_dates();
382
 
        }
383
 
        $temp_line = sprintf ("%-12s%-15s%13s %s%4s%-8s%-8s %3s %-s", 
384
 
                              'LOCUS', $seq->id(),$len,
385
 
                              (lc($mol) eq 'protein') ? ('aa','', '') : 
386
 
                              ('bp', '',$mol),$circular,
387
 
                              $div,$date);
388
 
    } 
389
 
    
390
 
    $self->_print("$temp_line\n");
391
 
    $self->_write_line_GenBank_regex("DEFINITION  ", "            ",
392
 
                                     $seq->desc(),"\\s\+\|\$",80);
393
 
    
394
 
    # if there, write the accession line
395
 
 
396
 
    if( $self->_ac_generation_func ) {
397
 
        $temp_line = &{$self->_ac_generation_func}($seq);
398
 
        $self->_print("ACCESSION   $temp_line\n");   
399
 
    } else {
400
 
        my @acc = ();
401
 
        push(@acc, $seq->accession_number());
402
 
        if( $seq->isa('Bio::Seq::RichSeqI') ) {
403
 
            push(@acc, $seq->get_secondary_accessions());
404
 
        }
405
 
        $self->_print("ACCESSION   ", join(" ", @acc), "\n");
406
 
        # otherwise - cannot print <sigh>
407
 
    } 
408
 
 
409
 
    # if PID defined, print it
410
 
    if($seq->isa('Bio::Seq::RichSeqI') && $seq->pid()) {
411
 
        $self->_print("PID         ", $seq->pid(), "\n");
412
 
    }
413
 
 
414
 
    # if there, write the version line
415
 
    
416
 
    if( defined $self->_sv_generation_func() ) {
417
 
        $temp_line = &{$self->_sv_generation_func}($seq);
418
 
        if( $temp_line ) {
419
 
            $self->_print("VERSION     $temp_line\n");   
420
 
        }
421
 
    } else {
422
 
        if($seq->isa('Bio::Seq::RichSeqI') && defined($seq->seq_version)) {
423
 
            my $id = $seq->primary_id(); # this may be a GI number
424
 
            $self->_print("VERSION     ",
425
 
                          $seq->accession_number(), ".", $seq->seq_version,
426
 
                          ($id && ($id =~ /^\d+$/) ? "  GI:".$id : ""),
427
 
                          "\n");
428
 
       }
429
 
    } 
430
 
    
431
 
    # if there, write the keywords line
432
 
    
433
 
    if( defined $self->_kw_generation_func() ) {
434
 
        $temp_line = &{$self->_kw_generation_func}($seq);
435
 
        $self->_print("KEYWORDS    $temp_line\n");   
436
 
    } else {
437
 
        if( $seq->can('keywords') ) {
438
 
            $self->_print("KEYWORDS    ",$seq->keywords,"\n");
439
 
        }
440
 
    } 
441
 
    
442
 
    
443
 
    # Organism lines
444
 
    if (my $spec = $seq->species) {
445
 
        my ($species, $genus, @class) = $spec->classification();
446
 
        my $OS;
447
 
        if( $spec->common_name ) {
448
 
            $OS = $spec->common_name;
449
 
        } else { 
450
 
            $OS = "$genus $species";
451
 
        }
452
 
        if (my $ssp = $spec->sub_species) {
453
 
            $OS .= " $ssp";
454
 
        }
455
 
        $self->_print("SOURCE      $OS.\n");
456
 
        $self->_print("  ORGANISM  ",
457
 
                      ($spec->organelle() ? $spec->organelle()." " : ""),
458
 
                      "$genus $species", "\n");
459
 
        my $OC = join('; ', (reverse(@class), $genus)) .'.';
460
 
        $self->_write_line_GenBank_regex(' 'x12,' 'x12,
461
 
                                         $OC,"\\s\+\|\$",80);
462
 
    }
463
 
    
464
 
    # Reference lines
465
 
    my $count = 1;
466
 
    foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
467
 
        $temp_line = sprintf ("REFERENCE   $count  (%s %d to %d)",
468
 
                              ($seq->alphabet() eq "protein" ?
469
 
                               "residues" : "bases"),
470
 
                              $ref->start,$ref->end);
 
588
    my ($self,@seqs) = @_;
 
589
 
 
590
    foreach my $seq ( @seqs ) {
 
591
        $self->throw("Attempting to write with no seq!") unless defined $seq;
 
592
 
 
593
        if( ! ref $seq || ! $seq->isa('Bio::SeqI') ) {
 
594
            $self->warn(" $seq is not a SeqI compliant module. Attempting to dump, but may fail!");
 
595
        }
 
596
 
 
597
        my $str = $seq->seq;
 
598
 
 
599
        my ($div, $mol);
 
600
        my $len = $seq->length();
 
601
 
 
602
        if ( $seq->can('division') ) {
 
603
            $div=$seq->division;
 
604
        } 
 
605
        if( !defined $div || ! $div ) { $div = 'UNK'; }
 
606
        my $alpha = $seq->alphabet;
 
607
        if( !$seq->can('molecule') || ! defined ($mol = $seq->molecule()) ) {
 
608
            $mol =  $alpha || 'DNA';
 
609
        }
 
610
 
 
611
        my $circular = 'linear  ';
 
612
        $circular = 'circular' if $seq->is_circular;
 
613
 
 
614
        local($^W) = 0; # supressing warnings about uninitialized fields.
 
615
 
 
616
        my $temp_line;
 
617
        if( $self->_id_generation_func ) {
 
618
            $temp_line = &{$self->_id_generation_func}($seq);
 
619
        } else {
 
620
            my $date = '';
 
621
            if( $seq->can('get_dates') ) {          
 
622
                ($date) = $seq->get_dates();
 
623
            }
 
624
 
 
625
            $self->warn("No whitespace allowed in GenBank display id [". $seq->display_id. "]")
 
626
                if $seq->display_id =~ /\s/;
 
627
 
 
628
            $temp_line = sprintf ("%-12s%-15s%13s %s%4s%-8s%-8s %3s %-s", 
 
629
                                  'LOCUS', $seq->id(),$len,
 
630
                                  (lc($alpha) eq 'protein') ? ('aa','', '') : 
 
631
                                  ('bp', '',$mol),$circular,
 
632
                                  $div,$date);
 
633
        } 
 
634
 
471
635
        $self->_print("$temp_line\n");
472
 
        $self->_write_line_GenBank_regex("  AUTHORS   ",' 'x12,
473
 
                                         $ref->authors,"\\s\+\|\$",80);
474
 
        $self->_write_line_GenBank_regex("  TITLE     "," "x12,
475
 
                                         $ref->title,"\\s\+\|\$",80);
476
 
        $self->_write_line_GenBank_regex("  JOURNAL   "," "x12,
477
 
                                         $ref->location,"\\s\+\|\$",80);
478
 
        if ($ref->comment) {
479
 
            $self->_write_line_GenBank_regex("  REMARK    "," "x12,
480
 
                                             $ref->comment,"\\s\+\|\$",80);
481
 
        }
482
 
        if( $ref->medline) {
483
 
            $self->_write_line_GenBank_regex("  MEDLINE   "," "x12,
484
 
                                             $ref->medline, "\\s\+\|\$",80);
485
 
             # I am assuming that pubmed entries only exist when there
486
 
             # are also MEDLINE entries due to the indentation
487
 
             # This could be a wrong assumption
488
 
             if( $ref->pubmed ) {
489
 
                 $self->_write_line_GenBank_regex("   PUBMED   "," "x12,
490
 
                                                  $ref->pubmed, "\\s\+\|\$",
491
 
                                                  80);
492
 
             }
493
 
        }
494
 
        $count++;
495
 
    }
496
 
    # Comment lines
497
 
    
498
 
    foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
499
 
        $self->_write_line_GenBank_regex("COMMENT     "," "x12,
500
 
                                         $comment->text,"\\s\+\|\$",80);
501
 
    }
502
 
    $self->_print("FEATURES             Location/Qualifiers\n");
503
 
    
504
 
    my $contig;
505
 
    if( defined $self->_post_sort ) {
506
 
        # we need to read things into an array. Process. Sort them. Print 'em
507
 
 
508
 
        my $post_sort_func = $self->_post_sort();
509
 
        my @fth;
510
 
 
511
 
        foreach my $sf ( $seq->top_SeqFeatures ) {
512
 
            push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
513
 
        }
514
 
 
515
 
        @fth = sort { &$post_sort_func($a,$b) } @fth;
516
 
 
517
 
        foreach my $fth ( @fth ) {
518
 
            $self->_print_GenBank_FTHelper($fth);
519
 
        }
520
 
    } else {
521
 
        # not post sorted. And so we can print as we get them.
522
 
        # lower memory load...
523
 
        
524
 
        foreach my $sf ( $seq->top_SeqFeatures ) {
525
 
            my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
 
636
        $self->_write_line_GenBank_regex("DEFINITION  ", "            ",
 
637
                                         $seq->desc(),"\\s\+\|\$",80);
 
638
 
 
639
        # if there, write the accession line
 
640
 
 
641
        if( $self->_ac_generation_func ) {
 
642
            $temp_line = &{$self->_ac_generation_func}($seq);
 
643
            $self->_print("ACCESSION   $temp_line\n");   
 
644
        } else {
 
645
            my @acc = ();
 
646
            push(@acc, $seq->accession_number());
 
647
            if( $seq->isa('Bio::Seq::RichSeqI') ) {
 
648
                push(@acc, $seq->get_secondary_accessions());
 
649
            }
 
650
            $self->_print("ACCESSION   ", join(" ", @acc), "\n");
 
651
            # otherwise - cannot print <sigh>
 
652
        } 
 
653
 
 
654
        # if PID defined, print it
 
655
        if($seq->isa('Bio::Seq::RichSeqI') && $seq->pid()) {
 
656
            $self->_print("PID         ", $seq->pid(), "\n");
 
657
        }
 
658
 
 
659
        # if there, write the version line
 
660
 
 
661
        if( defined $self->_sv_generation_func() ) {
 
662
            $temp_line = &{$self->_sv_generation_func}($seq);
 
663
            if( $temp_line ) {
 
664
                $self->_print("VERSION     $temp_line\n");   
 
665
            }
 
666
        } else {
 
667
            if($seq->isa('Bio::Seq::RichSeqI') && defined($seq->seq_version)) {
 
668
                my $id = $seq->primary_id(); # this may be a GI number
 
669
                $self->_print("VERSION     ",
 
670
                              $seq->accession_number(), ".", $seq->seq_version,
 
671
                              ($id && ($id =~ /^\d+$/) ? "  GI:".$id : ""),
 
672
                              "\n");
 
673
            }
 
674
        } 
 
675
 
 
676
        # if there, write the keywords line
 
677
 
 
678
        if( defined $self->_kw_generation_func() ) {
 
679
            $temp_line = &{$self->_kw_generation_func}($seq);
 
680
            $self->_print("KEYWORDS    $temp_line\n");   
 
681
        } else {
 
682
            if( $seq->can('keywords') ) {
 
683
                my $kw = $seq->keywords;
 
684
                $kw .= '.' if( $kw !~ /\.$/ );
 
685
                $self->_print("KEYWORDS    $kw\n");
 
686
            }
 
687
        } 
 
688
 
 
689
        # SEGMENT if it exists
 
690
        foreach my $ref ( $seq->annotation->get_Annotations('segment') ) {
 
691
            $self->_print(sprintf ("%-11s %s\n",'SEGMENT',
 
692
                                  $ref->value));
 
693
        }
 
694
 
 
695
        # Organism lines
 
696
        if (my $spec = $seq->species) {
 
697
            my ($species, $genus, @class) = $spec->classification();
 
698
            my $OS;
 
699
            if( $spec->common_name ) {
 
700
                $OS = $spec->common_name;
 
701
            } else { 
 
702
                $OS = "$genus $species";
 
703
            }
 
704
            if (my $ssp = $spec->sub_species) {
 
705
                $OS .= " $ssp";
 
706
            }
 
707
            $self->_print("SOURCE      $OS\n");
 
708
            $self->_print("  ORGANISM  ",
 
709
                          ($spec->organelle() ? $spec->organelle()." " : ""),
 
710
                          "$genus $species", "\n");
 
711
            my $OC = join('; ', (reverse(@class), $genus)) .'.';
 
712
            $self->_write_line_GenBank_regex(' 'x12,' 'x12,
 
713
                                             $OC,"\\s\+\|\$",80);
 
714
        }
 
715
 
 
716
        # Reference lines
 
717
        my $count = 1;
 
718
        foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
 
719
            $temp_line = "REFERENCE   $count";
 
720
            $temp_line .= sprintf ("  (%s %d to %d)",
 
721
                                  ($seq->alphabet() eq "protein" ?
 
722
                                   "residues" : "bases"),
 
723
                                  $ref->start,$ref->end)
 
724
                if $ref->start;
 
725
            $self->_print("$temp_line\n");
 
726
            $self->_write_line_GenBank_regex("  AUTHORS   ",' 'x12,
 
727
                                             $ref->authors,"\\s\+\|\$",80);
 
728
            $self->_write_line_GenBank_regex("  TITLE     "," "x12,
 
729
                                             $ref->title,"\\s\+\|\$",80);
 
730
            $self->_write_line_GenBank_regex("  JOURNAL   "," "x12,
 
731
                                             $ref->location,"\\s\+\|\$",80);
 
732
            if ($ref->comment) {
 
733
                $self->_write_line_GenBank_regex("  REMARK    "," "x12,
 
734
                                                 $ref->comment,"\\s\+\|\$",80);
 
735
            }
 
736
            if( $ref->medline) {
 
737
                $self->_write_line_GenBank_regex("  MEDLINE   "," "x12,
 
738
                                                 $ref->medline, "\\s\+\|\$",80);
 
739
                # I am assuming that pubmed entries only exist when there
 
740
                # are also MEDLINE entries due to the indentation
 
741
                # This could be a wrong assumption
 
742
                if( $ref->pubmed ) {
 
743
                    $self->_write_line_GenBank_regex("   PUBMED   "," "x12,
 
744
                                                     $ref->pubmed, "\\s\+\|\$",
 
745
                                                     80);
 
746
                }
 
747
            }
 
748
            $count++;
 
749
        }
 
750
        # Comment lines
 
751
 
 
752
        foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
 
753
            $self->_write_line_GenBank_regex("COMMENT     "," "x12,
 
754
                                             $comment->text,"\\s\+\|\$",80);
 
755
        }
 
756
        $self->_print("FEATURES             Location/Qualifiers\n");
 
757
 
 
758
        my $contig;
 
759
        if( defined $self->_post_sort ) {
 
760
            # we need to read things into an array. Process. Sort them. Print 'em
 
761
 
 
762
            my $post_sort_func = $self->_post_sort();
 
763
            my @fth;
 
764
 
 
765
            foreach my $sf ( $seq->top_SeqFeatures ) {
 
766
                push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
 
767
            }
 
768
 
 
769
            @fth = sort { &$post_sort_func($a,$b) } @fth;
 
770
 
526
771
            foreach my $fth ( @fth ) {
527
 
                if( ! $fth->isa('Bio::SeqIO::FTHelper') ) {
528
 
                    $sf->throw("Cannot process FTHelper... $fth");
529
 
                }               
530
772
                $self->_print_GenBank_FTHelper($fth);
531
773
            }
532
 
        }
533
 
    }
534
 
    if( $seq->length == 0 ) { $self->_show_dna(0) }
535
 
    
536
 
    if( $self->_show_dna() == 0 ) {
537
 
        $self->_print("\n//\n");
538
 
       return;
539
 
   }
540
 
    
 
774
        } else {
 
775
            # not post sorted. And so we can print as we get them.
 
776
            # lower memory load...
 
777
 
 
778
            foreach my $sf ( $seq->top_SeqFeatures ) {
 
779
                my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
 
780
                foreach my $fth ( @fth ) {
 
781
                    if( ! $fth->isa('Bio::SeqIO::FTHelper') ) {
 
782
                        $sf->throw("Cannot process FTHelper... $fth");
 
783
                    }           
 
784
                    $self->_print_GenBank_FTHelper($fth);
 
785
                }
 
786
            }
 
787
        }
 
788
        if( $seq->length == 0 ) { $self->_show_dna(0) }
 
789
 
 
790
        if( $self->_show_dna() == 0 ) {
 
791
            $self->_print("\n//\n");
 
792
            return;
 
793
        }
 
794
 
541
795
# finished printing features.
542
 
    
543
 
    $str =~ tr/A-Z/a-z/;
 
796
 
 
797
        $str =~ tr/A-Z/a-z/;
544
798
 
545
799
# Count each nucleotide
546
 
    unless(  $mol eq 'protein' ) {
547
 
        my $alen = $str =~ tr/a/a/;
548
 
        my $clen = $str =~ tr/c/c/;
549
 
        my $glen = $str =~ tr/g/g/;
550
 
        my $tlen = $str =~ tr/t/t/;
551
 
        
552
 
        my $olen = $len - ($alen + $tlen + $clen + $glen);
553
 
        if( $olen < 0 ) {
554
 
            $self->warn("Weird. More atgc than bases. Problem!");
555
 
        }
556
 
    
557
 
        my $base_count = sprintf("BASE COUNT %8s a %6s c %6s g %6s t%s\n",
558
 
                                 $alen,$clen,$glen,$tlen,
559
 
                                 ( $olen > 0 ) ? sprintf("%6s others",$olen) : '');
560
 
        $self->_print($base_count); 
561
 
    }
562
 
    $self->_print(sprintf("ORIGIN%6s\n",''));
563
 
    my $di;
564
 
    
565
 
    # push sequence blocks to be printed into an array
566
 
    # so we can comply with genbank format and not have 
567
 
    # block \n
568
 
    #     ^^ space after last block before newline is not allowed
569
 
    #        in the strictest sense of the format 
570
 
    my @seqline;
571
 
    for ($i = 0; $i < length($str); $i += 10) {
572
 
        
573
 
        $di=$i+11;
574
 
 
575
 
        #first line
576
 
        if ($i==0) {
577
 
            $self->_print(sprintf("%9d ",1));
578
 
        }
579
 
        #print sequence, spaced by 10
580
 
        push @seqline, substr($str,$i,10);
581
 
        #break line and print number at beginning of next line
582
 
        if(($i+10)%60 == 0) {       
583
 
            $self->_print(join(' ', @seqline), "\n");
584
 
            $self->_print(sprintf("%9d ",$di));
585
 
            @seqline = ();
586
 
        }
587
 
    }
588
 
    # fencepost problem, print whatever is left in the queue
589
 
    # again we are doing this to comply with  genbank in the 
590
 
    # strictest sense
591
 
    $self->_print(join(' ', @seqline), ' ') if( @seqline );
592
 
    
593
 
    
594
 
    $self->_print("\n//\n");
595
 
    return 1;
 
800
        unless(  $mol eq 'protein' ) {
 
801
            my $alen = $str =~ tr/a/a/;
 
802
            my $clen = $str =~ tr/c/c/;
 
803
            my $glen = $str =~ tr/g/g/;
 
804
            my $tlen = $str =~ tr/t/t/;
 
805
 
 
806
            my $olen = $len - ($alen + $tlen + $clen + $glen);
 
807
            if( $olen < 0 ) {
 
808
                $self->warn("Weird. More atgc than bases. Problem!");
 
809
            }
 
810
 
 
811
            my $base_count = sprintf("BASE COUNT %8s a %6s c %6s g %6s t%s\n",
 
812
                                     $alen,$clen,$glen,$tlen,
 
813
                                     ( $olen > 0 ) ? 
 
814
                                     sprintf("%6s others",$olen) : '');
 
815
            $self->_print($base_count); 
 
816
        }
 
817
 
 
818
        my ($o) = $seq->annotation->get_Annotations('origin');
 
819
        $self->_print(sprintf("%-6s%s\n",'ORIGIN',$o ? $o->value : ''));
 
820
        # print out the sequence
 
821
        my $nuc = 60;           # Number of nucleotides per line
 
822
        my $whole_pat = 'a10' x 6; # Pattern for unpacking a whole line
 
823
        my $out_pat   = 'A11' x 6; # Pattern for packing a line
 
824
        my $length = length($str);
 
825
 
 
826
        # Calculate the number of nucleotides which fit on whole lines
 
827
        my $whole = int($length / $nuc) * $nuc;
 
828
 
 
829
        # Print the whole lines
 
830
        my $i;
 
831
        for ($i = 0; $i < $whole; $i += $nuc) {
 
832
            my $blocks = pack $out_pat,
 
833
            unpack $whole_pat,
 
834
            substr($str, $i, $nuc);
 
835
            chop $blocks;
 
836
            $self->_print(sprintf("%9d $blocks\n", $i + $nuc - 59));
 
837
        }
 
838
 
 
839
        # Print the last line
 
840
        if (my $last = substr($str, $i)) {
 
841
            my $last_len = length($last);
 
842
            my $last_pat = 'a10' x int($last_len / 10) .'a'. $last_len % 10;
 
843
            my $blocks = pack $out_pat,
 
844
            unpack($last_pat, $last);
 
845
            $blocks =~ s/ +$//;
 
846
            $self->_print(sprintf("%9d $blocks\n", $length - $last_len + 1));
 
847
        }
 
848
 
 
849
        $self->_print("//\n");
 
850
 
 
851
        $self->flush if $self->_flush_on_write && defined $self->_fh;
 
852
        return 1;
 
853
    }
596
854
}
597
855
 
598
856
=head2 _print_GenBank_FTHelper
608
866
=cut
609
867
 
610
868
sub _print_GenBank_FTHelper {
611
 
   my ($self,$fth,$always_quote) = @_;
 
869
   my ($self,$fth) = @_;
612
870
   
613
871
   if( ! ref $fth || ! $fth->isa('Bio::SeqIO::FTHelper') ) {
614
872
       $fth->warn("$fth is not a FTHelper class. Attempting to print, but there could be tears!");   
623
881
                                        $fth->loc,"\,\|\$",80);
624
882
   }
625
883
 
626
 
   if( !defined $always_quote) { $always_quote = 0; }
627
 
   
628
884
   foreach my $tag ( keys %{$fth->field} ) {
629
885
       foreach my $value ( @{$fth->field->{$tag}} ) {
630
886
           $value =~ s/\"/\"\"/g;
633
889
                                                " "x21,
634
890
                                                "/$tag","\.\|\$",80);
635
891
           }
636
 
           elsif( $always_quote == 1 || $value !~ /^\d+$/ ) {
637
 
              my ($pat) = ($value =~ /\s/ ? '\s|$' : '.|$');          
 
892
           # there are almost 3x more quoted qualifier values and they
 
893
           # are more common too so we take quoted ones first
 
894
           elsif (!$FTQUAL_NO_QUOTE{$tag}) {
 
895
              my ($pat) = ($value =~ /\s/ ? '\s|$' : '.|$');
638
896
              $self->_write_line_GenBank_regex(" "x21,
639
897
                                               " "x21,
640
898
                                               "/$tag=\"$value\"",$pat,80);
 
899
 
641
900
           } else {
642
901
              $self->_write_line_GenBank_regex(" "x21,
643
902
                                               " "x21,
654
913
 Title   : _read_GenBank_References
655
914
 Usage   :
656
915
 Function: Reads references from GenBank format. Internal function really
657
 
 Example :
658
916
 Returns : 
659
917
 Args    :
660
918
 
676
934
 
677
935
   my (@title,@loc,@authors,@com,@medline,@pubmed);
678
936
 
679
 
   while( defined($_) || defined($_ = $self->_readline) ) {
680
 
       if (/^  AUTHORS\s+(.*)/) { 
 
937
   REFLOOP: while( defined($_) || defined($_ = $self->_readline) ) {
 
938
       if (/^\s{2}AUTHORS\s+(.*)/o) { 
681
939
           push (@authors, $1);   
682
940
           while ( defined($_ = $self->_readline) ) {
683
 
               /^\s{3,}(.*)/ && do { push (@authors, $1);next;};
 
941
               /^\s{3,}(.*)/o && do { push (@authors, $1);next;};
684
942
               last;
685
943
           }
686
944
           $ref->authors(join(' ', @authors));
687
945
       }
688
 
       if (/^  TITLE\s+(.*)/)  { 
 
946
       if (/^\s{2}TITLE\s+(.*)/o)  { 
689
947
           push (@title, $1);
690
948
           while ( defined($_ = $self->_readline) ) {
691
 
               /^\s{3,}(.*)/ && do { push (@title, $1);
 
949
               /^\s{3,}(.*)/o && do { push (@title, $1);
692
950
                                     next;
693
951
                                 };
694
952
               last;
695
953
           }
696
954
           $ref->title(join(' ', @title));
697
955
       }
698
 
       if (/^  JOURNAL\s+(.*)/) { 
 
956
       if (/^\s{2}JOURNAL\s+(.*)/o) { 
699
957
           push(@loc, $1);
700
958
           while ( defined($_ = $self->_readline) ) {
701
 
               /^\s{3,}(.*)/ && do { push(@loc, $1);
 
959
               /^\s{3,}(.*)/o && do { push(@loc, $1);
702
960
                                     next;
703
961
                                 };
704
962
               last;
705
963
           }
706
964
           $ref->location(join(' ', @loc));
 
965
           redo REFLOOP;
707
966
       }
708
 
       if (/^  REMARK\s+(.*)/) { 
 
967
       if (/^\s{2}REMARK\s+(.*)/o) { 
709
968
           push (@com, $1);
710
969
           while ( defined($_ = $self->_readline) ) {          
711
 
               /^\s{3,}(.*)/ && do { push(@com, $1);
 
970
               /^\s{3,}(.*)/o && do { push(@com, $1);
712
971
                                     next;
713
972
                                 };
714
973
               last;
715
974
           }
716
975
           $ref->comment(join(' ', @com));
 
976
           redo REFLOOP;
717
977
       }
718
 
       if( /^  MEDLINE\s+(.*)/ ) {
 
978
       if( /^\s{2}MEDLINE\s+(.*)/ ) {
719
979
           push(@medline,$1);
720
980
           while ( defined($_ = $self->_readline) ) {          
721
981
               /^\s{4,}(.*)/ && do { push(@medline, $1);
724
984
               last;
725
985
           }
726
986
           $ref->medline(join(' ', @medline));
727
 
       }       
728
 
       if( /^   PUBMED\s+(.*)/ ) {
 
987
           redo REFLOOP;
 
988
       }
 
989
       if( /^\s{3}PUBMED\s+(.*)/ ) {
729
990
           push(@pubmed,$1);
730
991
           while ( defined($_ = $self->_readline) ) {          
731
992
               /^\s{5,}(.*)/ && do { push(@pubmed, $1);
734
995
               last;
735
996
           }
736
997
           $ref->pubmed(join(' ', @pubmed));
 
998
           redo REFLOOP;
737
999
       }
738
1000
       
739
 
       /^REFERENCE/ && do {
 
1001
       /^REFERENCE/o && do {
740
1002
           # store current reference
741
1003
           $self->_add_ref_to_array(\@refs,$ref) if $ref;
742
1004
           # reset
755
1017
           }
756
1018
       };
757
1019
 
758
 
       /^(FEATURES)|(COMMENT)/ && last;
 
1020
       /^(FEATURES)|(COMMENT)/o && last;
759
1021
 
760
1022
       $_ = undef; # Empty $_ to trigger read of next line
761
1023
   }
804
1066
           lines.
805
1067
 Example :
806
1068
 Returns : A Bio::Species object
807
 
 Args    :
 
1069
 Args    : a reference to the current line buffer
808
1070
 
809
1071
=cut
810
1072
 
811
1073
sub _read_GenBank_Species {
812
1074
    my( $self,$buffer) = @_;
813
1075
    my @organell_names = ("chloroplast", "mitochondr"); 
814
 
             # only those carrying DNA, apart from the nucleus
 
1076
    # only those carrying DNA, apart from the nucleus
815
1077
 
816
1078
    $_ = $$buffer;
817
1079
    
818
 
    my( $sub_species, $species, $genus, $common, $organelle, @class );
 
1080
    my( $sub_species, $species, $genus, $common, $organelle, @class, $ns_name );
819
1081
    # upon first entering the loop, we must not read a new line -- the SOURCE
820
1082
    # line is already in the buffer (HL 05/10/2000)
821
1083
    while (defined($_) || defined($_ = $self->_readline())) {
822
1084
        # de-HTMLify (links that may be encountered here don't contain
823
1085
        # escaped '>', so a simple-minded approach suffices)
824
1086
        s/<[^>]+>//g;
825
 
        if (/^SOURCE\s+(.*)/) {
 
1087
        if (/^SOURCE\s+(.*)/o) {
826
1088
            # FIXME this is probably mostly wrong (e.g., it yields things like
827
1089
            # Homo sapiens adult placenta cDNA to mRNA
828
1090
            # which is certainly not what you want)
829
1091
            $common = $1;
830
1092
            $common =~ s/\.$//; # remove trailing dot
831
 
        } elsif (/^\s+ORGANISM/) {
 
1093
        } elsif (/^\s{2}ORGANISM/o) {
832
1094
            my @spflds = split(' ', $_);
 
1095
            ($ns_name) = $_ =~ /\w+\s+(.*)/o;
833
1096
            shift(@spflds); # ORGANISM
834
1097
            if(grep { $_ =~ /^$spflds[0]/i; } @organell_names) {
835
1098
                $organelle = shift(@spflds);
841
1104
                $species = "sp.";
842
1105
            }
843
1106
            $sub_species = shift(@spflds) if(@spflds);
844
 
        } elsif (/^\s+(.+)/) {
 
1107
        } elsif (/^\s+(.+)/o) {
845
1108
            # only split on ';' or '.' so that 
846
1109
            # classification that is 2 words will 
847
1110
            # still get matched
857
1120
    $$buffer = $_;
858
1121
    
859
1122
    # Don't make a species object if it's empty or "Unknown" or "None"
860
 
    return unless $genus and  $genus !~ /^(Unknown|None)$/i;
 
1123
    return unless $genus and  $genus !~ /^(Unknown|None)$/oi;
861
1124
    
862
1125
    # Bio::Species array needs array in Species -> Kingdom direction
863
 
    if ($class[$#class] eq $genus) {
 
1126
    if ($class[0] eq 'Viruses') {
 
1127
        push( @class, $ns_name );
 
1128
    }
 
1129
    elsif ($class[$#class] eq $genus) {
864
1130
        push( @class, $species );
865
1131
    } else {
866
1132
        push( @class, $genus, $species );
868
1134
    @class = reverse @class;
869
1135
    
870
1136
    my $make = Bio::Species->new();
871
 
    $make->classification( @class );
 
1137
    $make->classification( \@class, "FORCE" ); # no name validation please
872
1138
    $make->common_name( $common      ) if $common;
873
 
    $make->sub_species( $sub_species ) if $sub_species;
 
1139
    unless ($class[-1] eq 'Viruses') {
 
1140
        $make->sub_species( $sub_species ) if $sub_species;
 
1141
    }
874
1142
    $make->organelle($organelle) if $organelle;
875
1143
    return $make;
876
1144
}
895
1163
        );
896
1164
    my @qual = (); # An arrray of lines making up the qualifiers
897
1165
    
898
 
    if ($$buffer =~ /^     (\S+)\s+(\S+)/) {
 
1166
    if ($$buffer =~ /^\s{5}(\S+)\s+(.+?)\s*$/o) {
899
1167
        $key = $1;
900
1168
        $loc = $2;
901
1169
        # Read all the lines up to the next feature
902
1170
        while ( defined($_ = $self->_readline) ) {
903
 
            if (/^(\s+)(.+?)\s*$/) {
904
 
                # Lines inside features are preceeded by 21 spaces
905
 
                # A new feature is preceeded by 5 spaces
 
1171
            if (/^(\s+)(.+?)\s*$/o) {
 
1172
                # Lines inside features are preceded by 21 spaces
 
1173
                # A new feature is preceded by 5 spaces
906
1174
                if (length($1) > 6) {
907
1175
                    # Add to qualifiers if we're in the qualifiers, or if it's
908
1176
                    # the first qualifier
909
 
                    if (($#qual >= 0) || (substr($2, 0, 1) eq '/')) {
 
1177
                    if (@qual || (index($2,'/') == 0)) {
910
1178
                        push(@qual, $2);
911
1179
                    }
912
1180
                    # We're still in the location line, so append to location
925
1193
    } else {
926
1194
        # No feature key
927
1195
        $self->debug("no feature key!\n");
 
1196
        # change suggested by JDiggans to avoid infinite loop- 
 
1197
        # see bugreport 1062.
 
1198
        # reset buffer to prevent infinite loop
 
1199
        $$buffer = $self->_readline();
928
1200
        return;
929
1201
    } 
930
1202
    
954
1226
                # and the quotes are balanced
955
1227
                while ($value !~ /\"$/ or $value =~ tr/"/"/ % 2) {
956
1228
                    if($i >= $#qual) {
957
 
                        # We haven't found the closing douple quote ...
958
 
                        # Even though this should be considered as a malformed
959
 
                        # entry, we allow for a single exception, namely if the
960
 
                        # value ends exactly at the last char position of a
961
 
                        # GenBank line.
962
 
                        # At least 78 chars required, of which 21 are spaces
963
 
                        #Currently disabled, let's wait for an actual sequence
964
 
                        #entry that really requires this.
965
 
                        #if(length($qual[$i]) >= 57) {
966
 
                        #    $self->warn("unbalanced quotes in feature $key ".
967
 
                        #               "(location: $loc), ".
968
 
                        #               "qualifier $qualifier, ".
969
 
                        #               "accepting though");
970
 
                        #    last;
971
 
                        #} else {
972
 
                            $self->warn("Unbalanced quote in:\n" .
973
 
                                        join('', map("$_\n", @qual)) .
974
 
                                        "No further qualifiers will " .
975
 
                                        "be added for this feature");
976
 
                            last QUAL;
977
 
                        #}
 
1229
                       $self->warn("Unbalanced quote in:\n" .
 
1230
                                   join("\n", @qual) .
 
1231
                                   "No further qualifiers will " .
 
1232
                                   "be added for this feature");
 
1233
                       last QUAL;
978
1234
                    }
979
1235
                    $i++; # modifying a for-loop variable inside of the loop
980
1236
                          # is not the best programming style ...
981
1237
                    my $next = $qual[$i];
982
1238
 
983
 
                    # Join to value with space if value or next line contains a space
984
 
                    $value .= (grep /\s/, ($value, $next)) ? " $next" : $next;
 
1239
                    # add to value with a space unless the value appears
 
1240
                    # to be a sequence (translation for example)
 
1241
                    if(($value.$next) =~ /[^A-Za-z\"\-]/o) {
 
1242
                        $value .= " ";
 
1243
                    }
 
1244
                    $value .= $next;
985
1245
                }
986
1246
                # Trim leading and trailing quotes
987
1247
                $value =~ s/^"|"$//g;
988
1248
                # Undouble internal quotes
989
1249
                $value =~ s/""/\"/g;
990
 
            }
 
1250
            } elsif ( $value =~ /^\(/ ) { # values quoted by ()s
 
1251
                # Keep addingto value until we find the trailing bracket
 
1252
                # and the ()s are balanced
 
1253
                my $left = ($value =~ tr/\(/\(/); # count left parens
 
1254
                my $right = ($value =~ tr/\)/\)/); # count right parens
 
1255
                while( $value !~ /\)$/ or $left != $right ) {
 
1256
                    if( $i >= $#qual) {
 
1257
                        $self->warn("Unbalanced parens in:\n".
 
1258
                                    join("\n", @qual).
 
1259
                                    "No further qualifiers will ".
 
1260
                                    "be added for this feature");
 
1261
                        last QUAL;
 
1262
                    }
 
1263
                    $i++;
 
1264
                    my $next = $qual[$i];                   
 
1265
                    $value .= $next;
 
1266
                    $left += ($next =~ tr/\(/\(/);
 
1267
                    $right += ($next =~ tr/\)/\)/);
 
1268
                }
 
1269
            }
991
1270
        } else {
992
1271
            $value = '_no_value';
993
1272
        }
1052
1331
 
1053
1332
   $length || $self->throw( "Miscalled write_line_GenBank without length. Programming error!");
1054
1333
 
1055
 
   if( length $pre1 != length $pre2 ) {
1056
 
       $self->throw( "Programming error - cannot called write_line_GenBank_regex with different length pre1 and pre2 tags!");
1057
 
   }
 
1334
#   if( length $pre1 != length $pre2 ) {
 
1335
#       $self->throw( "Programming error - cannot called write_line_GenBank_regex with different length pre1 and pre2 tags!");
 
1336
#   }
1058
1337
 
1059
1338
   my $subl = $length - (length $pre1) - 2;
1060
 
   my @lines;
 
1339
   my @lines = ();
1061
1340
 
1062
 
   while($line =~ m/(.{1,$subl})($regex)/g) {
1063
 
       # be strict about not padding spaces according to 
1064
 
       # genbank format
1065
 
       my $l = $1.$2;
1066
 
       $l =~ s/\s+$//;
1067
 
       push(@lines, $l);
 
1341
   CHUNK: while($line) {
 
1342
       foreach my $pat ($regex, '[,;\.\/-]\s|'.$regex, '[,;\.\/-]|'.$regex) {
 
1343
           if($line =~ m/^(.{1,$subl})($pat)(.*)/) {
 
1344
               $line = $3;
 
1345
               # be strict about not padding spaces according to 
 
1346
               # genbank format
 
1347
               my $l = $1.$2;
 
1348
               $l =~ s/\s+$//;
 
1349
               push(@lines, $l);
 
1350
               next CHUNK;
 
1351
           }
 
1352
       }
 
1353
       # if we get here none of the patterns matched $subl or less chars
 
1354
       $self->warn("trouble dissecting \"$line\" into chunks ".
 
1355
                   "of $subl chars or less - this tag won't print right");
 
1356
       # insert a space char to prevent infinite loops
 
1357
       $line = substr($line,0,$subl) . " " . substr($line,$subl);
1068
1358
   }
1069
1359
   
1070
1360
   my $s = shift @lines;