← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/SeqIO/embl.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/SeqIO/embl.pm
1
 
# $Id: embl.pm,v 1.69 2003/12/22 18:33:15 heikki Exp $
 
1
# $Id: embl.pm,v 1.92.4.5 2006/10/31 15:48:58 cjfields Exp $
2
2
#
3
3
# BioPerl module for Bio::SeqIO::EMBL
4
4
#
17
17
=head1 SYNOPSIS
18
18
 
19
19
It is probably best not to use this object directly, but
20
 
rather go through the AnnSeqIO handler system. Go:
 
20
rather go through the SeqIO handler system. Go:
21
21
 
22
22
    $stream = Bio::SeqIO->new(-file => $filename, -format => 'EMBL');
23
23
 
24
24
    while ( (my $seq = $stream->next_seq()) ) {
25
 
        # do something with $seq
 
25
            # do something with $seq
26
26
    }
27
27
 
28
28
=head1 DESCRIPTION
30
30
This object can transform Bio::Seq objects to and from EMBL flat
31
31
file databases.
32
32
 
33
 
There is alot of flexibility here about how to dump things which I need
34
 
to document fully.
 
33
There is a lot of flexibility here about how to dump things which
 
34
should be documented more fully.
35
35
 
36
 
There should be a common object that this and genbank share (probably
37
 
with swissprot). Too much of the magic is identical. 
 
36
There should be a common object that this and Genbank share (probably
 
37
with Swissprot). Too much of the magic is identical.
38
38
 
39
39
=head2 Optional functions
40
40
 
51
51
 
52
52
=item _id_generation_func()
53
53
 
54
 
This is function which is called as 
 
54
This is function which is called as
55
55
 
56
56
   print "ID   ", $func($annseq), "\n";
57
57
 
58
58
To generate the ID line. If it is not there, it generates a sensible ID
59
59
line using a number of tools.
60
60
 
61
 
If you want to output annotations in embl format they need to be
 
61
If you want to output annotations in EMBL format they need to be
62
62
stored in a Bio::Annotation::Collection object which is accessible
63
 
through the Bio::SeqI interface method L<annotation()|annotation>.  
 
63
through the Bio::SeqI interface method L<annotation()|annotation>.
64
64
 
65
65
The following are the names of the keys which are polled from a
66
66
L<Bio::Annotation::Collection> object.
67
67
 
68
 
reference       - Should contain Bio::Annotation::Reference objects
69
 
comment         - Should contain Bio::Annotation::Comment objects
70
 
dblink          - Should contain Bio::Annotation::DBLink objects
 
68
 reference  - Should contain Bio::Annotation::Reference objects
 
69
 comment    - Should contain Bio::Annotation::Comment objects
 
70
 dblink     - Should contain Bio::Annotation::DBLink objects
71
71
 
72
72
=back
73
73
 
75
75
 
76
76
=head2 Mailing Lists
77
77
 
78
 
User feedback is an integral part of the evolution of this
79
 
and other Bioperl modules. Send your comments and suggestions preferably
80
 
 to one of the Bioperl mailing lists.
81
 
Your participation is much appreciated.
 
78
User feedback is an integral part of the evolution of this and other
 
79
Bioperl modules. Send your comments and suggestions preferably to one
 
80
of the Bioperl mailing lists.  Your participation is much appreciated.
82
81
 
83
 
  bioperl-l@bioperl.org                 - General discussion
84
 
  http://www.bioperl.org/MailList.shtml - About the mailing lists
 
82
  bioperl-l@bioperl.org                  - General discussion
 
83
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
85
84
 
86
85
=head2 Reporting Bugs
87
86
 
88
87
Report bugs to the Bioperl bug tracking system to help us keep track
89
 
 the bugs and their resolution.
90
 
 Bug reports can be submitted via email or the web:
 
88
the bugs and their resolution. Bug reports can be submitted via
 
89
the web:
91
90
 
92
 
  bioperl-bugs@bio.perl.org
93
 
  http://bugzilla.bioperl.org/
 
91
  http://bugzilla.open-bio.org/
94
92
 
95
93
=head1 AUTHOR - Ewan Birney
96
94
 
97
95
Email birney@ebi.ac.uk
98
96
 
99
 
Describe contact details here
100
 
 
101
97
=head1 APPENDIX
102
98
 
103
99
The rest of the documentation details each of the object
110
106
 
111
107
 
112
108
package Bio::SeqIO::embl;
113
 
use vars qw(@ISA %FTQUAL_NO_QUOTE);
 
109
use vars qw(%FTQUAL_NO_QUOTE);
114
110
use strict;
115
111
use Bio::SeqIO::FTHelper;
116
112
use Bio::SeqFeature::Generic;
121
117
use Bio::Annotation::Reference;
122
118
use Bio::Annotation::DBLink;
123
119
 
124
 
@ISA = qw(Bio::SeqIO);
 
120
use base qw(Bio::SeqIO);
125
121
 
126
122
%FTQUAL_NO_QUOTE=(
127
123
  'anticodon'=>1,
144
140
sub _initialize {
145
141
  my($self,@args) = @_;
146
142
 
147
 
  $self->SUPER::_initialize(@args);  
 
143
  $self->SUPER::_initialize(@args);
148
144
  # hash for functions for decoding keys.
149
 
  $self->{'_func_ftunit_hash'} = {}; 
150
 
  $self->_show_dna(1); # sets this to one by default. People can change it 
 
145
  $self->{'_func_ftunit_hash'} = {};
 
146
  $self->_show_dna(1); # sets this to one by default. People can change it
151
147
  if( ! defined $self->sequence_factory ) {
152
148
      $self->sequence_factory(new Bio::Seq::SeqFactory
153
 
                              (-verbose => $self->verbose(), 
 
149
                              (-verbose => $self->verbose(),
154
150
                               -type => 'Bio::Seq::RichSeq'));
155
151
  }
156
152
}
166
162
=cut
167
163
 
168
164
sub next_seq {
169
 
   my ($self,@args) = @_;
170
 
   my ($pseq,$c,$line,$name,$desc,$acc,$seqc,$mol,$div, 
171
 
       $date, $comment, @date_arr);
172
 
 
173
 
   my ($annotation, %params, @features) = ( new Bio::Annotation::Collection);
174
 
 
175
 
   $line = $self->_readline;   # This needs to be before the first eof() test
176
 
 
177
 
   if( !defined $line ) {
178
 
       return undef; # no throws - end of file
179
 
   }
180
 
 
181
 
   if( $line =~ /^\s+$/ ) {
182
 
       while( defined ($line = $self->_readline) ) {
183
 
           $line =~/^\S/ && last;
184
 
       }
185
 
   }   
186
 
   if( !defined $line ) {
187
 
       return undef; # end of file
188
 
   }
189
 
   $line =~ /^ID\s+\S+/ || $self->throw("EMBL stream with no ID. Not embl in my book");
190
 
   $line =~ /^ID\s+(\S+)\s+\S+\;\s+([^;]+)\;\s+(\S+)\;/;
191
 
   $name = $1;
192
 
   $mol = $2;
193
 
   $div = $3;
194
 
   if(! $name) {
195
 
       $name = "unknown id";
196
 
   }
197
 
   my $alphabet;
198
 
 
199
 
    # this is important to have the id for display in e.g. FTHelper, otherwise
200
 
    # you won't know which entry caused an error
201
 
   if($mol) {
202
 
       if ( $mol =~ /circular/ ) {
203
 
           $params{'-is_circular'} = 1;
204
 
           $mol =~  s|circular ||;
205
 
       }
206
 
       if (defined $mol ) {        
207
 
           if ($mol =~ /DNA/) {
208
 
               $alphabet='dna';
209
 
           }
210
 
           elsif ($mol =~ /RNA/) {
211
 
               $alphabet='rna';
212
 
           }
213
 
           elsif ($mol =~ /AA/) {
214
 
               $alphabet='protein';
215
 
           }
216
 
       }
217
 
 
218
 
   }
219
 
   
220
 
#   $self->warn("not parsing upper annotation in EMBL file yet!");
 
165
    my ($self,@args) = @_;
 
166
    my ($pseq,$c,$line,$name,$desc,$acc,$seqc,$mol,$div,
 
167
        $date, $comment, @date_arr);
 
168
 
 
169
    my ($annotation, %params, @features) =
 
170
       new Bio::Annotation::Collection;
 
171
 
 
172
    $line = $self->_readline;
 
173
    # This needs to be before the first eof() test
 
174
 
 
175
    if( !defined $line ) {
 
176
        return; # no throws - end of file
 
177
    }
 
178
 
 
179
    if( $line =~ /^\s+$/ ) {
 
180
        while( defined ($line = $self->_readline) ) {
 
181
            $line =~/^\S/ && last;
 
182
        }
 
183
        # return without error if the whole next sequence was just a single
 
184
        # blank line and then eof
 
185
        return unless $line;
 
186
    }
 
187
 
 
188
    # no ID as 1st non-blank line, need short circuit and exit routine
 
189
    $self->throw("EMBL stream with no ID. Not embl in my book")
 
190
       unless $line =~ /^ID\s+\S+/;
 
191
 
 
192
        # At this point we are sure that $line contains an ID header line
 
193
        my $alphabet;
 
194
    if ( $line =~ tr/;/;/ == 6) {   # New style headers contain exactly six semicolons.
 
195
 
 
196
        # New style header (EMBL Release >= 87, after June 2006)
 
197
        my $topology;
 
198
        my $sv;
 
199
 
 
200
        # ID   DQ299383; SV 1; linear; mRNA; STD; MAM; 431 BP.
 
201
                # This regexp comes from the new2old.pl conversion script, from EBI
 
202
        $line =~ m/^ID   (\w+);\s+SV (\d+); (\w+); ([^;]+); (\w{3}); (\w{3}); (\d+) BP./;
 
203
        ($name, $sv, $topology, $mol, $div) = ($1, $2, $3, $4, $6);
 
204
 
 
205
        if (defined($sv)) {
 
206
            $params{'-seq_version'} = $sv;
 
207
            $params{'-version'} = $sv;
 
208
        }
 
209
 
 
210
        if ($topology eq "circular") {
 
211
            $params{'-is_circular'} = 1;
 
212
        }
 
213
        
 
214
        if (defined $mol ) {
 
215
            if ($mol =~ /DNA/) {
 
216
                $alphabet='dna';
 
217
            }
 
218
            elsif ($mol =~ /RNA/) {
 
219
                $alphabet='rna';
 
220
            }
 
221
            elsif ($mol =~ /AA/) {
 
222
                $alphabet='protein';
 
223
            }
 
224
        }
 
225
    }
 
226
    else {
 
227
        
 
228
        # Old style header (EMBL Release < 87, before June 2006)
 
229
        ($name, $mol, $div) = ($line =~ /^ID\s+(\S+)[^;]*;\s+(\S+)[^;]*;\s+(\S+)[^;]*;/);
 
230
        
 
231
        if($mol) {
 
232
            if ( $mol =~ /circular/ ) {
 
233
                $params{'-is_circular'} = 1;
 
234
                $mol =~  s|circular ||;
 
235
            }
 
236
            if (defined $mol ) {
 
237
                if ($mol =~ /DNA/) {
 
238
                    $alphabet='dna';
 
239
                }
 
240
                elsif ($mol =~ /RNA/) {
 
241
                    $alphabet='rna';
 
242
                }
 
243
                elsif ($mol =~ /AA/) {
 
244
                    $alphabet='protein';
 
245
                }
 
246
            }
 
247
        }
 
248
    }
 
249
 
 
250
    unless( defined $name && length($name) ) {
 
251
        $name = "unknown_id";
 
252
    }
 
253
 
 
254
        # $self->warn("not parsing upper annotation in EMBL file yet!");
221
255
   my $buffer = $line;
222
 
 
223
 
   local $_;
224
 
 
225
 
   BEFORE_FEATURE_TABLE :
226
 
   until( !defined $buffer ) {
227
 
       $_ = $buffer;
228
 
 
229
 
       # Exit at start of Feature table
230
 
       last if /^F[HT]/;
231
 
 
232
 
       # Description line(s)
233
 
       if (/^DE\s+(\S.*\S)/) {
234
 
           $desc .= $desc ? " $1" : $1;
235
 
       }
236
 
 
237
 
       #accession number
238
 
       if( /^AC\s+(.*)?/ ) {
239
 
           my @accs = split(/[; ]+/, $1); # allow space in addition
240
 
           $params{'-accession_number'} = shift @accs 
241
 
               unless defined $params{'-accession_number'};
242
 
           push @{$params{'-secondary_accessions'}}, @accs;
243
 
       }
244
 
       
245
 
       #version number
246
 
       if( /^SV\s+\S+\.(\d+);?/ ) {
247
 
           my $sv = $1;
248
 
           #$sv =~ s/\;//;
249
 
           $params{'-seq_version'} = $sv;
250
 
           $params{'-version'} = $sv;
251
 
       }
252
 
 
253
 
       #date (NOTE: takes last date line)
254
 
       if( /^DT\s+(.+)$/ ) {
255
 
           my $date = $1;
256
 
           push @{$params{'-dates'}},$date;
257
 
       }
258
 
       
259
 
       #keywords
260
 
       if( /^KW   (.*)\S*$/ ) {
261
 
           my @kw = split(/\s*\;\s*/,$1);
262
 
           
263
 
           push @{$params{'-keywords'}}, @kw;
264
 
       }
265
 
 
266
 
       # Organism name and phylogenetic information
267
 
       elsif (/^O[SC]/) {
268
 
           my $species = $self->_read_EMBL_Species(\$buffer);
269
 
           $params{'-species'}= $species;
270
 
       }
271
 
 
272
 
       # References
273
 
       elsif (/^R/) {
274
 
           my @refs = $self->_read_EMBL_References(\$buffer);
275
 
           foreach my $ref ( @refs ) {
276
 
               $annotation->add_Annotation('reference',$ref);
277
 
           }
278
 
       }
279
 
       
280
 
       # DB Xrefs
281
 
       elsif (/^DR/) {
282
 
           my @links = $self->_read_EMBL_DBLink(\$buffer);
283
 
           foreach my $dblink ( @links ) {
284
 
               $annotation->add_Annotation('dblink',$dblink);
285
 
           }
286
 
       }
287
 
       
288
 
       # Comments
289
 
       elsif (/^CC\s+(.*)/) {
290
 
           $comment .= $1;
291
 
           $comment .= " ";
292
 
           while (defined ($_ = $self->_readline) ) {
293
 
               if (/^CC\s+(.*)/) {
294
 
                   $comment .= $1;
295
 
                   $comment .= " ";
296
 
               }
297
 
               else { 
298
 
                   last;
299
 
               }
300
 
           }
301
 
           my $commobj = Bio::Annotation::Comment->new();
302
 
           $commobj->text($comment);
303
 
           $annotation->add_Annotation('comment',$commobj);
304
 
           $comment = "";
305
 
       }
306
 
 
307
 
       # Get next line.
308
 
       $buffer = $self->_readline;
309
 
   }
 
256
    local $_;
 
257
    BEFORE_FEATURE_TABLE :
 
258
        until( !defined $buffer ) {
 
259
            $_ = $buffer;
 
260
            # Exit at start of Feature table
 
261
            if( /^(F[HT]|SQ)/ ) {
 
262
                $self->_pushback($_) if( $1 eq 'SQ' || $1 eq 'FT');
 
263
                last;
 
264
            }
 
265
            # Description line(s)
 
266
            if (/^DE\s+(\S.*\S)/) {
 
267
                $desc .= $desc ? " $1" : $1;
 
268
            }
 
269
 
 
270
            #accession number
 
271
            if( /^AC\s+(.*)?/ ) {
 
272
                my @accs = split(/[; ]+/, $1); # allow space in addition
 
273
                $params{'-accession_number'} = shift @accs
 
274
                    unless defined $params{'-accession_number'};
 
275
                push @{$params{'-secondary_accessions'}}, @accs;
 
276
            }
 
277
 
 
278
            #version number
 
279
            if( /^SV\s+\S+\.(\d+);?/ ) {
 
280
                my $sv = $1;
 
281
                #$sv =~ s/\;//;
 
282
                $params{'-seq_version'} = $sv;
 
283
                $params{'-version'} = $sv;
 
284
            }
 
285
 
 
286
            #date (NOTE: takes last date line)
 
287
            if( /^DT\s+(.+)$/ ) {
 
288
                my $line = $1;
 
289
                my ($date, $version) = split(' ', $line, 2);
 
290
                $date =~ tr/,//d; # remove comma if new version
 
291
                if ($version =~ /\(Rel\. (\d+), Created\)/xms ) {
 
292
                    my $release = Bio::Annotation::SimpleValue->new(
 
293
                                                                    -tagname    => 'creation_release',
 
294
                                                                    -value      => $1
 
295
                                                                    );
 
296
                    $annotation->add_Annotation($release);
 
297
                } elsif ($version =~ /\(Rel\. (\d+), Last updated, Version (\d+)\)/xms ) {
 
298
                    my $release = Bio::Annotation::SimpleValue->new(
 
299
                                                                    -tagname    => 'update_release',
 
300
                                                                    -value      => $1
 
301
                                                                    );
 
302
                    $annotation->add_Annotation($release);
 
303
 
 
304
                    my $update = Bio::Annotation::SimpleValue->new(
 
305
                                                                   -tagname    => 'update_version',
 
306
                                                                   -value      => $2
 
307
                                                                   );
 
308
                    $annotation->add_Annotation($update);
 
309
                }
 
310
                push @{$params{'-dates'}}, $date;
 
311
            }
 
312
 
 
313
            #keywords
 
314
            if( /^KW   (.*)\S*$/ ) {
 
315
                my @kw = split(/\s*\;\s*/,$1);
 
316
                push @{$params{'-keywords'}}, @kw;
 
317
            }
 
318
 
 
319
            # Organism name and phylogenetic information
 
320
            elsif (/^O[SC]/) {
 
321
                # pass the accession number so we can give an informative throw message if necessary
 
322
                my $species = $self->_read_EMBL_Species(\$buffer, $params{'-accession_number'});
 
323
                $params{'-species'}= $species;
 
324
            }
 
325
 
 
326
            # References
 
327
            elsif (/^R/) {
 
328
                my @refs = $self->_read_EMBL_References(\$buffer);
 
329
                foreach my $ref ( @refs ) {
 
330
                    $annotation->add_Annotation('reference',$ref);
 
331
                }
 
332
            }
 
333
 
 
334
            # DB Xrefs
 
335
            elsif (/^DR/) {
 
336
                my @links = $self->_read_EMBL_DBLink(\$buffer);
 
337
                foreach my $dblink ( @links ) {
 
338
                    $annotation->add_Annotation('dblink',$dblink);
 
339
                }
 
340
            }
 
341
 
 
342
            # Comments
 
343
            elsif (/^CC\s+(.*)/) {
 
344
                $comment .= $1;
 
345
                $comment .= " ";
 
346
                while (defined ($_ = $self->_readline) ) {
 
347
                    if (/^CC\s+(.*)/) {
 
348
                        $comment .= $1;
 
349
                        $comment .= " ";
 
350
                    }
 
351
                    else {
 
352
                        last;
 
353
                    }
 
354
                }
 
355
                my $commobj = Bio::Annotation::Comment->new();
 
356
                $commobj->text($comment);
 
357
                $annotation->add_Annotation('comment',$commobj);
 
358
                $comment = "";
 
359
            }
 
360
 
 
361
            # Get next line.
 
362
            $buffer = $self->_readline;
 
363
        }
310
364
 
311
365
   while( defined ($_ = $self->_readline) ) {
312
 
       /^FT   \w/ && last;
313
 
       /^SQ / && last;
314
 
       /^CO / && last;
 
366
                /^FT\s{3}\w/ && last;
 
367
                /^SQ / && last;
 
368
                /^CO / && last;
315
369
   }
316
370
   $buffer = $_;
317
371
 
318
372
   if (defined($buffer) && $buffer =~ /^FT /) {
319
 
     until( !defined ($buffer) ) {
320
 
         my $ftunit = $self->_read_FTHelper_EMBL(\$buffer);
321
 
         # process ftunit
322
 
 
323
 
         push(@features,
324
 
              $ftunit->_generic_seqfeature($self->location_factory(), $name));
325
 
 
326
 
         if( $buffer !~ /^FT/ ) {
327
 
             last;
328
 
         }
329
 
     }
 
373
                until( !defined ($buffer) ) {
 
374
                        my $ftunit = $self->_read_FTHelper_EMBL(\$buffer);
 
375
 
 
376
                        # process ftunit
 
377
         my $feat =
 
378
                          $ftunit->_generic_seqfeature($self->location_factory(), $name);
 
379
 
 
380
         # add taxon_id from source if available
 
381
         if($params{'-species'} && ($feat->primary_tag eq 'source')
 
382
            && $feat->has_tag('db_xref')
 
383
            && (! $params{'-species'}->ncbi_taxid())) {
 
384
                                foreach my $tagval ($feat->get_tag_values('db_xref')) {
 
385
                                        if(index($tagval,"taxon:") == 0) {
 
386
                                                $params{'-species'}->ncbi_taxid(substr($tagval,6));
 
387
                                                last;
 
388
                                        }
 
389
                                }
 
390
         }
 
391
 
 
392
         # add feature to list of features
 
393
                        push(@features, $feat);
 
394
 
 
395
                        if( $buffer !~ /^FT/ ) {
 
396
                                last;
 
397
                        }
 
398
                }
330
399
   }
331
 
 
332
 
 
333
400
   # skip comments
334
 
   while( defined ($buffer) && $buffer =~ /^XX/ ) { 
 
401
   while( defined ($buffer) && $buffer =~ /^XX/ ) {
335
402
       $buffer = $self->_readline();
336
 
   }
337
 
   
338
 
   if( $buffer =~ /^CO/  ) {       
339
 
       until( !defined ($buffer) ) {
340
 
           my $ftunit = $self->_read_FTHelper_EMBL(\$buffer);
341
 
           # process ftunit
342
 
           push(@features,
343
 
                $ftunit->_generic_seqfeature($self->location_factory(),
344
 
                                             $name));
345
 
           
346
 
           if( $buffer !~ /^CO/ ) {
347
 
               last;
348
 
           }
349
 
       }       
 
403
         }
 
404
 
 
405
   if( $buffer =~ /^CO/  ) {
 
406
                until( !defined ($buffer) ) {
 
407
                        my $ftunit = $self->_read_FTHelper_EMBL(\$buffer);
 
408
                        # process ftunit
 
409
                        push(@features,
 
410
                             $ftunit->_generic_seqfeature($self->location_factory(),
 
411
                                                          $name));
 
412
 
 
413
                        if( $buffer !~ /^CO/ ) {
 
414
                                last;
 
415
                        }
 
416
                }
350
417
   }
351
418
   if( $buffer !~ /^SQ/  ) {
352
 
       while( defined ($_ = $self->_readline) ) {
353
 
           /^SQ/ && last;
354
 
       }
 
419
                while( defined ($_ = $self->_readline) ) {
 
420
                        /^SQ/ && last;
 
421
                }
355
422
   }
356
 
   
357
 
   $seqc = "";         
 
423
   $seqc = "";
358
424
   while( defined ($_ = $self->_readline) ) {
359
 
       /^\/\// && last;
360
 
       $_ = uc($_);
361
 
       s/[^A-Za-z]//g;
362
 
       $seqc .= $_;
 
425
                m{^//} && last;
 
426
                $_ = uc($_);
 
427
                s/[^A-Za-z]//g;
 
428
                $seqc .= $_;
363
429
   }
364
430
   my $seq = $self->sequence_factory->create
365
 
       (-verbose => $self->verbose(),
366
 
        -division => $div,
367
 
        -seq => $seqc,
368
 
        -desc => $desc,
369
 
        -display_id => $name,
370
 
        -annotation => $annotation,
371
 
        -molecule => $mol,
372
 
        -alphabet => $alphabet,
373
 
        -features => \@features,
374
 
        %params);
375
 
 
 
431
          (-verbose => $self->verbose(),
 
432
                -division => $div,
 
433
                -seq => $seqc,
 
434
                -desc => $desc,
 
435
                -display_id => $name,
 
436
                -annotation => $annotation,
 
437
                -molecule => $mol,
 
438
                -alphabet => $alphabet,
 
439
                -features => \@features,
 
440
                %params);
376
441
   return $seq;
377
442
}
378
443
 
 
444
 
 
445
 
 
446
=head2 _write_ID_line
 
447
 
 
448
 Title   : _write_ID_line
 
449
 Usage   : $self->_write_ID_line($seq);
 
450
 Function: Writes the EMBL Release 87 format ID line to the stream, unless
 
451
         : there is a user-supplied ID line generation function in which
 
452
         : case that is used instead.
 
453
         : ( See Bio::SeqIO::embl::_id_generation_function(). )
 
454
 Returns : nothing
 
455
 Args    : Bio::Seq object
 
456
 
 
457
=cut
 
458
 
 
459
sub _write_ID_line {
 
460
 
 
461
        my ($self, $seq) = @_;
 
462
 
 
463
        my $id_line;
 
464
        # If there is a user-supplied ID generation function, use it.
 
465
        if( $self->_id_generation_func ) {
 
466
                $id_line = "ID   " . &{$self->_id_generation_func}($seq) . "\nXX\n";
 
467
        }
 
468
        # Otherwise, generate a standard EMBL release 87 (June 2006) ID line.
 
469
        else {
 
470
 
 
471
                # The sequence name is supposed to be the primary accession number,
 
472
                my $name = $seq->accession_number();
 
473
                if (!$name) {
 
474
                        # but if it is not present, use the sequence ID.
 
475
                        $name = $seq->id();
 
476
                }
 
477
 
 
478
                $self->warn("No whitespace allowed in EMBL id [". $name. "]") if $name =~ /\s/;
 
479
 
 
480
                # Use the sequence version, or default to 1.
 
481
                my $version = $seq->version() || 1;
 
482
 
 
483
                my $len = $seq->length();
 
484
 
 
485
                # Taxonomic division.
 
486
                my $div;
 
487
                if ( $seq->can('division') && defined($seq->division) && $self->_is_valid_division($seq->division) ) {
 
488
                        $div = $seq->division();
 
489
                }
 
490
                else {
 
491
                        $div ||= 'UNC';                 # 'UNC' is the EMBL division code for 'unclassified'.
 
492
                }
 
493
 
 
494
                my $mol;
 
495
                # If the molecule type is a valid EMBL type, use it.
 
496
                if (  $seq->can('molecule')
 
497
                      && defined($seq->molecule)
 
498
                      && $self->_is_valid_molecule_type($seq->molecule)
 
499
                    )
 
500
                {
 
501
                        $mol = $seq->molecule();
 
502
                }
 
503
                # Otherwise, choose unassigned DNA or RNA based on the alphabet.
 
504
                elsif ($seq->can('primary_seq') && defined $seq->primary_seq->alphabet) {
 
505
                        my $alphabet =$seq->primary_seq->alphabet;
 
506
                        if ($alphabet eq 'dna') {
 
507
                                $mol ='unassigned DNA';
 
508
                        }
 
509
                        elsif ($alphabet eq 'rna') {
 
510
                                $mol='unassigned RNA';
 
511
                        }
 
512
                        elsif ($alphabet eq 'protein') {
 
513
                                $self->warn("Protein sequence found; EMBL is a nucleotide format.");
 
514
                                $mol='AA';      # AA is not a valid EMBL molecule type.
 
515
                        }
 
516
                }
 
517
 
 
518
                my $topology = 'linear';
 
519
                if ($seq->is_circular) {
 
520
                        $topology = 'circular';
 
521
                }
 
522
 
 
523
        $mol ||= '';# 'unassigned'; ?
 
524
                $id_line = "ID   $name; SV $version; $topology; $mol; STD; $div; $len BP.\nXX\n";
 
525
                $self->_print($id_line);
 
526
        }
 
527
}
 
528
 
 
529
=head2 _is_valid_division
 
530
 
 
531
 Title   : _is_valid_division
 
532
 Usage   : $self->_is_valid_division($div)
 
533
 Function: tests division code for validity
 
534
 Returns : true if $div is a valid EMBL release 87 taxonomic division.
 
535
 Args    : taxonomic division code string
 
536
 
 
537
=cut
 
538
 
 
539
sub _is_valid_division {
 
540
        my ($self, $division) = @_;
 
541
 
 
542
        my %EMBL_divisions = (
 
543
                "PHG"    => 1,                  # Bacteriophage
 
544
                "ENV"    => 1,                  # Environmental Sample
 
545
                "FUN"    => 1,                  # Fungal
 
546
                "HUM"    => 1,                  # Human
 
547
                "INV"    => 1,                  # Invertebrate
 
548
                "MAM"    => 1,                  # Other Mammal
 
549
                "VRT"    => 1,                  # Other Vertebrate
 
550
                "MUS"    => 1,                  # Mus musculus
 
551
                "PLN"    => 1,                  # Plant
 
552
                "PRO"    => 1,                  # Prokaryote
 
553
            "ROD"    => 1,                      # Other Rodent
 
554
            "SYN"    => 1,                      # Synthetic
 
555
            "UNC"    => 1,                      # Unclassified
 
556
            "VRL"    => 1                       # Viral
 
557
        );
 
558
 
 
559
        return exists($EMBL_divisions{$division});
 
560
}
 
561
 
 
562
=head2 _is_valid_molecule_type
 
563
 
 
564
 Title   : _is_valid_molecule_type
 
565
 Usage   : $self->_is_valid_molecule_type($mol)
 
566
 Function: tests molecule type for validity
 
567
 Returns : true if $mol is a valid EMBL release 87 molecule type.
 
568
 Args    : molecule type string
 
569
 
 
570
=cut
 
571
 
 
572
sub _is_valid_molecule_type {
 
573
        my ($self, $moltype) = @_;
 
574
 
 
575
        my %EMBL_molecule_types = (
 
576
                "genomic DNA"    => 1,
 
577
                "genomic RNA"    => 1,
 
578
                "mRNA"           => 1,
 
579
                "tRNA"           => 1,
 
580
                "rRNA"           => 1,
 
581
                "snoRNA"         => 1,
 
582
                "snRNA"          => 1,
 
583
                "scRNA"          => 1,
 
584
                "pre-RNA"        => 1,
 
585
                "other RNA"      => 1,
 
586
            "other DNA"      => 1,
 
587
            "unassigned DNA" => 1,
 
588
            "unassigned RNA" => 1
 
589
        );
 
590
 
 
591
        return exists($EMBL_molecule_types{$moltype});
 
592
}
 
593
 
379
594
=head2 write_seq
380
595
 
381
596
 Title   : write_seq
382
597
 Usage   : $stream->write_seq($seq)
383
598
 Function: writes the $seq object (must be seq) to the stream
384
 
 Returns : 1 for success and 0 for error
 
599
 Returns : 1 for success and undef for error
385
600
 Args    : array of 1 to n Bio::SeqI objects
386
601
 
387
602
 
388
603
=cut
389
604
 
390
605
sub write_seq {
391
 
    my ($self,@seqs) = @_;
392
 
 
393
 
    foreach my $seq ( @seqs ) { 
394
 
        $self->throw("Attempting to write with no seq!") unless defined $seq;
395
 
        unless ( ref $seq && $seq->isa('Bio::SeqI' ) ) {
396
 
            $self->warn("$seq is not a SeqI compliant sequence object!") 
397
 
                if $self->verbose >= 0;
398
 
            unless ( ref $seq && $seq->isa('Bio::PrimarySeqI' ) ) {
399
 
                $self->throw("$seq is not a PrimarySeqI compliant sequence object!");
400
 
            }
401
 
        }
402
 
        my $str = $seq->seq || '';
403
 
 
404
 
        my $mol;
405
 
        my $div;
406
 
        my $len = $seq->length();
407
 
 
408
 
        if ($seq->can('division') && defined $seq->division) {
409
 
            $div = $seq->division();
410
 
        }
411
 
        $div ||= 'UNK';
412
 
 
413
 
        if ($seq->can('molecule')) {
414
 
            $mol = $seq->molecule();
415
 
            $mol = 'RNA' if defined $mol && $mol =~ /RNA/; # no 'mRNA' 
416
 
        }
417
 
        elsif ($seq->can('primary_seq') && defined $seq->primary_seq->alphabet) {
418
 
            my $alphabet =$seq->primary_seq->alphabet;
419
 
            if ($alphabet eq 'dna') {
420
 
                $mol ='DNA';
421
 
            }
422
 
            elsif ($alphabet eq 'rna') {
423
 
                $mol='RNA';
424
 
            }
425
 
            elsif ($alphabet eq 'protein') {
426
 
                $mol='AA';
427
 
            }
428
 
        }
429
 
        $mol ||= 'XXX';
430
 
        $mol = "circular $mol" if $seq->is_circular;
431
 
 
432
 
        my $temp_line;
433
 
        if( $self->_id_generation_func ) {
434
 
            $temp_line = &{$self->_id_generation_func}($seq);
435
 
        } else {
436
 
 
437
 
            $self->warn("No whitespace allowed in EMBL display id [". $seq->display_id. "]")
438
 
                if $seq->display_id =~ /\s/;
439
 
 
440
 
            $temp_line = sprintf("%-11sstandard; $mol; $div; %d BP.", $seq->id(), $len);
441
 
        } 
442
 
 
443
 
        $self->_print( "ID   $temp_line\n","XX\n");
444
 
 
445
 
        # Write the accession line if present
446
 
        my( $acc );
447
 
        {
448
 
            if( my $func = $self->_ac_generation_func ) {
449
 
                $acc = &{$func}($seq);
450
 
            } elsif( $seq->isa('Bio::Seq::RichSeqI') && 
451
 
                     defined($seq->accession_number) ) {
452
 
                $acc = $seq->accession_number;
453
 
                $acc = join("; ", $acc, $seq->get_secondary_accessions);
454
 
            } elsif ( $seq->can('accession_number') ) {
455
 
                $acc = $seq->accession_number;
456
 
            }
457
 
 
458
 
            if (defined $acc) {
459
 
                $self->_print("AC   $acc;\n",
460
 
                              "XX\n");
461
 
            }
462
 
        }
463
 
 
464
 
        # Write the sv line if present
465
 
        {
466
 
            my( $sv );
467
 
            if (my $func = $self->_sv_generation_func) {
468
 
                $sv = &{$func}($seq);
469
 
            } elsif($seq->isa('Bio::Seq::RichSeqI') && 
470
 
                    defined($seq->seq_version)) {
471
 
                my ($prim_acc) = $acc =~ /([^;]+)/;
472
 
                $sv = "$prim_acc.". $seq->seq_version();
473
 
            }   
474
 
            if (defined $sv) {
475
 
                $self->_print( "SV   $sv\n",
476
 
                               "XX\n");
477
 
            }
478
 
        }
479
 
 
480
 
        # Date lines
481
 
        my $switch=0;
482
 
        if( $seq->can('get_dates') ) {
483
 
            foreach my $dt ( $seq->get_dates() ) {
484
 
                $self->_write_line_EMBL_regex("DT   ","DT   ",$dt,'\s+|$',80);#'
485
 
                    $switch=1;
486
 
            }
487
 
            if ($switch == 1) {
488
 
                $self->_print("XX\n");
489
 
            }
490
 
        }
491
 
 
492
 
        # Description lines
493
 
        $self->_write_line_EMBL_regex("DE   ","DE   ",$seq->desc(),'\s+|$',80); #'
494
 
        $self->_print( "XX\n");
495
 
 
496
 
        # if there, write the kw line
497
 
        {
498
 
            my( $kw );
499
 
            if( my $func = $self->_kw_generation_func ) {
500
 
                $kw = &{$func}($seq);
501
 
            } elsif( $seq->can('keywords') ) {
502
 
                $kw = $seq->keywords;
503
 
            }
504
 
            if (defined $kw) {
505
 
                $self->_write_line_EMBL_regex("KW   ", "KW   ", $kw, '\s+|$', 80); #'
506
 
                $self->_print( "XX\n");
507
 
            }
508
 
        }
509
 
 
510
 
        # Organism lines
511
 
 
512
 
        if ($seq->can('species') && (my $spec = $seq->species)) {
513
 
            my($species, @class) = $spec->classification();
514
 
            my $genus = $class[0];
515
 
            my $OS = "$genus $species";
516
 
            if (my $ssp = $spec->sub_species) {
517
 
                $OS .= " $ssp";
518
 
            }
519
 
            if (my $common = $spec->common_name) {
520
 
                $OS .= " ($common)";
521
 
            }
522
 
            $self->_print("OS   $OS\n");
523
 
            my $OC = join('; ', reverse(@class)) .'.';
524
 
            $self->_write_line_EMBL_regex("OC   ","OC   ",$OC,'; |$',80); #'
525
 
            if ($spec->organelle) {
526
 
                $self->_write_line_EMBL_regex("OG   ","OG   ",$spec->organelle,'; |$',80); #'
527
 
            }
528
 
            $self->_print("XX\n");
529
 
        }
530
 
 
531
 
        # Reference lines
532
 
        my $t = 1;
533
 
        if ( $seq->can('annotation') && defined $seq->annotation ) {
534
 
            foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
535
 
                $self->_print( "RN   [$t]\n");
536
 
 
537
 
                # Having no RP line is legal, but we need both
538
 
                # start and end for a valid location.
539
 
                my $start = $ref->start;
540
 
                my $end   = $ref->end;
541
 
                if ($start and $end) {
542
 
                    $self->_print( "RP   $start-$end\n");
543
 
                } elsif ($start or $end) {
544
 
                    $self->throw("Both start and end are needed for a valid RP line.  Got: start='$start' end='$end'");
545
 
                }
546
 
 
547
 
                if (my $med = $ref->medline) {
548
 
                    $self->_print( "RX   MEDLINE; $med.\n");
549
 
                }
550
 
                if (my $pm = $ref->pubmed) {
551
 
                    $self->_print( "RX   PUBMED; $pm.\n");
552
 
                }
553
 
                $self->_write_line_EMBL_regex("RA   ", "RA   ", 
554
 
                                              $ref->authors . ";",
555
 
                                              '\s+|$', 80); #'
556
 
 
557
 
                # If there is no title to the reference, it appears
558
 
                # as a single semi-colon.  All titles must end in
559
 
                # a semi-colon.
560
 
                my $ref_title = $ref->title || '';
561
 
                $ref_title =~ s/[\s;]*$/;/;
562
 
                $self->_write_line_EMBL_regex("RT   ", "RT   ", $ref_title,    '\s+|$', 80); #'
563
 
                $self->_write_line_EMBL_regex("RL   ", "RL   ", $ref->location, '\s+|$', 80); #'
564
 
                if ($ref->comment) {
565
 
                    $self->_write_line_EMBL_regex("RC   ", "RC   ", $ref->comment, '\s+|$', 80); #' 
566
 
                }
567
 
                $self->_print("XX\n");
568
 
                $t++;
569
 
            }
570
 
 
571
 
 
572
 
            # DB Xref lines
573
 
            if (my @db_xref = $seq->annotation->get_Annotations('dblink') ) {
574
 
                foreach my $dr (@db_xref) {
575
 
                    my $db_name = $dr->database;
576
 
                    my $prim    = $dr->primary_id;
577
 
                    my $opt     = $dr->optional_id || '';
578
 
 
579
 
                    my $line = "$db_name; $prim; $opt.";
580
 
                    $self->_write_line_EMBL_regex("DR   ", "DR   ", $line, '\s+|$', 80); #'
581
 
                }
582
 
                $self->_print("XX\n");
583
 
            }
584
 
 
585
 
            # Comment lines
586
 
            foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
587
 
                $self->_write_line_EMBL_regex("CC   ", "CC   ", $comment->text, '\s+|$', 80); #'
588
 
                $self->_print("XX\n");
589
 
            }
590
 
        }
591
 
        # "\\s\+\|\$"
592
 
 
593
 
        ## FEATURE TABLE
594
 
 
595
 
        $self->_print("FH   Key             Location/Qualifiers\n");
596
 
        $self->_print("FH\n");
597
 
 
598
 
        my @feats = $seq->can('top_SeqFeatures') ? $seq->top_SeqFeatures : ();
599
 
        if ($feats[0]) {
600
 
            if( defined $self->_post_sort ) {
601
 
                # we need to read things into an array. 
602
 
                # Process. Sort them. Print 'em
603
 
 
604
 
                my $post_sort_func = $self->_post_sort();
605
 
                my @fth;
606
 
 
607
 
                foreach my $sf ( @feats ) {
608
 
                    push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
609
 
                }
610
 
 
611
 
                @fth = sort { &$post_sort_func($a,$b) } @fth;
612
 
 
613
 
                foreach my $fth ( @fth ) {
614
 
                    $self->_print_EMBL_FTHelper($fth);
615
 
                }
616
 
            } else {
617
 
                # not post sorted. And so we can print as we get them.
618
 
                # lower memory load...
619
 
 
620
 
                foreach my $sf ( @feats ) {
621
 
                    my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
622
 
                    foreach my $fth ( @fth ) {
623
 
                        if( $fth->key eq 'CONTIG') {
624
 
                            $self->_show_dna(0);
625
 
                        }
626
 
                        $self->_print_EMBL_FTHelper($fth);
627
 
                    }
628
 
                }
629
 
            }
630
 
        }
631
 
 
632
 
        if( $self->_show_dna() == 0 ) {
633
 
            $self->_print( "//\n");
634
 
            return;
635
 
        }
636
 
        $self->_print( "XX\n");
637
 
 
638
 
        # finished printing features.
639
 
        
640
 
        $str =~ tr/A-Z/a-z/;
641
 
 
642
 
        # Count each nucleotide
643
 
        my $alen = $str =~ tr/a/a/;
644
 
        my $clen = $str =~ tr/c/c/;
645
 
        my $glen = $str =~ tr/g/g/;
646
 
        my $tlen = $str =~ tr/t/t/;
647
 
 
648
 
        my $olen = $len - ($alen + $tlen + $clen + $glen);
649
 
        if( $olen < 0 ) {
650
 
            $self->warn("Weird. More atgc than bases. Problem!");
651
 
        }
652
 
 
653
 
        $self->_print("SQ   Sequence $len BP; $alen A; $clen C; $glen G; $tlen T; $olen other;\n");
654
 
 
655
 
        my $nuc = 60;           # Number of nucleotides per line
656
 
        my $whole_pat = 'a10' x 6; # Pattern for unpacking a whole line
657
 
        my $out_pat   = 'A11' x 6; # Pattern for packing a line
658
 
        my $length = length($str);
659
 
 
660
 
        # Calculate the number of nucleotides which fit on whole lines
661
 
        my $whole = int($length / $nuc) * $nuc;
662
 
 
663
 
        # Print the whole lines
664
 
        my( $i );
665
 
        for ($i = 0; $i < $whole; $i += $nuc) {
666
 
            my $blocks = pack $out_pat,
667
 
            unpack $whole_pat,
668
 
            substr($str, $i, $nuc);
669
 
            $self->_print(sprintf("     $blocks%9d\n", $i + $nuc));
670
 
        }
671
 
 
672
 
        # Print the last line
673
 
        if (my $last = substr($str, $i)) {
674
 
            my $last_len = length($last);
675
 
            my $last_pat = 'a10' x int($last_len / 10) .'a'. $last_len % 10;
676
 
            my $blocks = pack $out_pat,
677
 
            unpack($last_pat, $last);
678
 
            $self->_print(sprintf("     $blocks%9d\n", $length)); # Add the length to the end
679
 
        }
680
 
 
681
 
        $self->_print( "//\n");
682
 
 
683
 
        $self->flush if $self->_flush_on_write && defined $self->_fh;
 
606
        my ($self,@seqs) = @_;
 
607
 
 
608
        foreach my $seq ( @seqs ) {
 
609
                $self->throw("Attempting to write with no seq!") unless defined $seq;
 
610
                unless ( ref $seq && $seq->isa('Bio::SeqI' ) ) {
 
611
                        $self->warn("$seq is not a SeqI compliant sequence object!")
 
612
                          if $self->verbose >= 0;
 
613
                        unless ( ref $seq && $seq->isa('Bio::PrimarySeqI' ) ) {
 
614
                                $self->throw("$seq is not a PrimarySeqI compliant sequence object!");
 
615
                        }
 
616
                }
 
617
                my $str = $seq->seq || '';
 
618
 
 
619
                # Write the ID line.
 
620
                $self->_write_ID_line($seq);
 
621
 
 
622
 
 
623
                # Write the accession line if present
 
624
                my( $acc );
 
625
                {
 
626
                        if( my $func = $self->_ac_generation_func ) {
 
627
                                $acc = &{$func}($seq);
 
628
                        } elsif( $seq->isa('Bio::Seq::RichSeqI') &&
 
629
                                                defined($seq->accession_number) ) {
 
630
                                $acc = $seq->accession_number;
 
631
                                $acc = join("; ", $acc, $seq->get_secondary_accessions);
 
632
                        } elsif ( $seq->can('accession_number') ) {
 
633
                                $acc = $seq->accession_number;
 
634
                        }
 
635
 
 
636
                        if (defined $acc) {
 
637
                                $self->_print("AC   $acc;\n",
 
638
                                                                  "XX\n") || return;
 
639
                        }
 
640
                }
 
641
 
 
642
                # Date lines
 
643
                my $switch=0;
 
644
                if( $seq->can('get_dates') ) {
 
645
            my @dates =  $seq->get_dates();
 
646
            my $ct = 1;
 
647
            my $date_flag = 0;
 
648
            my ($cr) = $seq->get_Annotations("creation_release");
 
649
            my ($ur) = $seq->get_Annotations("update_release");
 
650
            my ($uv) = $seq->get_Annotations("update_version");
 
651
 
 
652
            unless ($cr && $ur && $ur) {
 
653
                $date_flag = 1;
 
654
            }
 
655
 
 
656
            foreach my $dt (@dates){
 
657
                if (!$date_flag) {
 
658
                    $self->_write_line_EMBL_regex("DT   ","DT   ",
 
659
                            $dt." (Rel. $cr, Created)",
 
660
                            '\s+|$',80) if $ct == 1;
 
661
                    $self->_write_line_EMBL_regex("DT   ","DT   ",
 
662
                            $dt." (Rel. $ur, Last updated, Version $uv)",
 
663
                            '\s+|$',80) if $ct == 2;
 
664
                } else { # other formats?
 
665
                    $self->_write_line_EMBL_regex("DT   ","DT   ",
 
666
                            $dt,'\s+|$',80);
 
667
                }
 
668
                $switch =1;
 
669
                $ct++;
 
670
            }
 
671
                        if ($switch == 1) {
 
672
                                $self->_print("XX\n") || return;
 
673
                        }
 
674
                }
 
675
 
 
676
                # Description lines
 
677
                $self->_write_line_EMBL_regex("DE   ","DE   ",$seq->desc(),'\s+|$',80) || return; #'
 
678
                $self->_print( "XX\n") || return;
 
679
 
 
680
                # if there, write the kw line
 
681
                {
 
682
                        my( $kw );
 
683
                        if( my $func = $self->_kw_generation_func ) {
 
684
                                $kw = &{$func}($seq);
 
685
                        } elsif( $seq->can('keywords') ) {
 
686
                                $kw = $seq->keywords;
 
687
                        }
 
688
                        if (defined $kw) {
 
689
                                $self->_write_line_EMBL_regex("KW   ", "KW   ", $kw, '\s+|$', 80) || return; #'
 
690
                                $self->_print( "XX\n") || return;
 
691
                        }
 
692
                }
 
693
 
 
694
                # Organism lines
 
695
 
 
696
                if ($seq->can('species') && (my $spec = $seq->species)) {
 
697
                        my @class = $spec->classification();
 
698
            shift @class; # get rid of species name. Some embl files include
 
699
                          # the species name in the OC lines, but this seems
 
700
                          # more like an error than something we need to
 
701
                          # emulate
 
702
                        my $OS = $spec->scientific_name;
 
703
            if ($spec->common_name) {
 
704
                $OS .= ' ('.$spec->common_name.')';
 
705
            }
 
706
                        $self->_print("OS   $OS\n") || return;
 
707
                        my $OC = join('; ', reverse(@class)) .'.';
 
708
                        $self->_write_line_EMBL_regex("OC   ","OC   ",$OC,'; |$',80) || return;
 
709
                        if ($spec->organelle) {
 
710
                                $self->_write_line_EMBL_regex("OG   ","OG   ",$spec->organelle,'; |$',80) || return;
 
711
                        }
 
712
                        $self->_print("XX\n") || return;
 
713
                }
 
714
 
 
715
                # Reference lines
 
716
                my $t = 1;
 
717
                if ( $seq->can('annotation') && defined $seq->annotation ) {
 
718
                        foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
 
719
                                $self->_print( "RN   [$t]\n") || return;
 
720
 
 
721
                                # Having no RP line is legal, but we need both
 
722
                                # start and end for a valid location.
 
723
                                my $start = $ref->start;
 
724
                                my $end   = $ref->end;
 
725
                                if ($start and $end) {
 
726
                                        $self->_print( "RP   $start-$end\n") || return;
 
727
                                } elsif ($start or $end) {
 
728
                                        $self->throw("Both start and end are needed for a valid RP line.  Got: start='$start' end='$end'");
 
729
                                }
 
730
 
 
731
                                if (my $med = $ref->medline) {
 
732
                                        $self->_print( "RX   MEDLINE; $med.\n") || return;
 
733
                                }
 
734
                                if (my $pm = $ref->pubmed) {
 
735
                                        $self->_print( "RX   PUBMED; $pm.\n") || return;
 
736
                                }
 
737
                                $self->_write_line_EMBL_regex("RA   ", "RA   ",
 
738
                                                                                                                $ref->authors . ";",
 
739
                                                                                                                '\s+|$', 80) || return; #'
 
740
 
 
741
                                # If there is no title to the reference, it appears
 
742
                                # as a single semi-colon.  All titles must end in
 
743
                                # a semi-colon.
 
744
                                my $ref_title = $ref->title || '';
 
745
                                $ref_title =~ s/[\s;]*$/;/;
 
746
                                $self->_write_line_EMBL_regex("RT   ", "RT   ", $ref_title,    '\s+|$', 80) || return; #'
 
747
                                $self->_write_line_EMBL_regex("RL   ", "RL   ", $ref->location, '\s+|$', 80) || return; #'
 
748
                                if ($ref->comment) {
 
749
                                        $self->_write_line_EMBL_regex("RC   ", "RC   ", $ref->comment, '\s+|$', 80) || return; #'
 
750
                                }
 
751
                                $self->_print("XX\n") || return;
 
752
                                $t++;
 
753
                        }
 
754
 
 
755
                        # DB Xref lines
 
756
                        if (my @db_xref = $seq->annotation->get_Annotations('dblink') ) {
 
757
                            for my $dr (@db_xref) {
 
758
                                my $db_name = $dr->database;
 
759
                                my $prim    = $dr->primary_id;
 
760
 
 
761
                                my $opt     = $dr->optional_id || '';
 
762
                                my $line = $opt ? "$db_name; $prim; $opt." : "$db_name; $prim.";
 
763
                                $self->_write_line_EMBL_regex("DR   ", "DR   ", $line, '\s+|$', 80) || return; #'
 
764
                            }
 
765
                            $self->_print("XX\n") || return;
 
766
                        }
 
767
                        
 
768
                        # Comment lines
 
769
                        foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
 
770
                                $self->_write_line_EMBL_regex("CC   ", "CC   ", $comment->text, '\s+|$', 80) || return; #'
 
771
                                $self->_print("XX\n") || return;
 
772
                        }
 
773
                }
 
774
                # "\\s\+\|\$"
 
775
 
 
776
                ## FEATURE TABLE
 
777
 
 
778
                $self->_print("FH   Key             Location/Qualifiers\n") || return;
 
779
                $self->_print("FH\n") || return;
 
780
 
 
781
                my @feats = $seq->can('top_SeqFeatures') ? $seq->top_SeqFeatures : ();
 
782
                if ($feats[0]) {
 
783
                        if( defined $self->_post_sort ) {
 
784
                                # we need to read things into an array.
 
785
                                # Process. Sort them. Print 'em
 
786
 
 
787
                                my $post_sort_func = $self->_post_sort();
 
788
                                my @fth;
 
789
 
 
790
                                foreach my $sf ( @feats ) {
 
791
                                        push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
 
792
                                }
 
793
 
 
794
                                @fth = sort { &$post_sort_func($a,$b) } @fth;
 
795
 
 
796
                                foreach my $fth ( @fth ) {
 
797
                                        $self->_print_EMBL_FTHelper($fth) || return;
 
798
                                }
 
799
                        } else {
 
800
                                # not post sorted. And so we can print as we get them.
 
801
                                # lower memory load...
 
802
 
 
803
                                foreach my $sf ( @feats ) {
 
804
                                        my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
 
805
                                        foreach my $fth ( @fth ) {
 
806
                                                if( $fth->key eq 'CONTIG') {
 
807
                                                        $self->_show_dna(0);
 
808
                                                }
 
809
                                                $self->_print_EMBL_FTHelper($fth) || return;
 
810
                                        }
 
811
                                }
 
812
                        }
 
813
                }
 
814
 
 
815
                if( $self->_show_dna() == 0 ) {
 
816
                        $self->_print( "//\n") || return;
 
817
                        return;
 
818
                }
 
819
                $self->_print( "XX\n") || return;
 
820
 
 
821
                # finished printing features.
 
822
 
 
823
                $str =~ tr/A-Z/a-z/;
 
824
 
 
825
                # Count each nucleotide
 
826
                my $alen = $str =~ tr/a/a/;
 
827
                my $clen = $str =~ tr/c/c/;
 
828
                my $glen = $str =~ tr/g/g/;
 
829
                my $tlen = $str =~ tr/t/t/;
 
830
 
 
831
                my $len = $seq->length();
 
832
                my $olen = $seq->length() - ($alen + $tlen + $clen + $glen);
 
833
                if( $olen < 0 ) {
 
834
                        $self->warn("Weird. More atgc than bases. Problem!");
 
835
                }
 
836
 
 
837
                $self->_print("SQ   Sequence $len BP; $alen A; $clen C; $glen G; $tlen T; $olen other;\n") || return;
 
838
 
 
839
                my $nuc = 60;           # Number of nucleotides per line
 
840
                my $whole_pat = 'a10' x 6; # Pattern for unpacking a whole line
 
841
                my $out_pat   = 'A11' x 6; # Pattern for packing a line
 
842
                my $length = length($str);
 
843
 
 
844
                # Calculate the number of nucleotides which fit on whole lines
 
845
                my $whole = int($length / $nuc) * $nuc;
 
846
 
 
847
                # Print the whole lines
 
848
                my( $i );
 
849
                for ($i = 0; $i < $whole; $i += $nuc) {
 
850
                        my $blocks = pack $out_pat,
 
851
                          unpack $whole_pat,
 
852
                                 substr($str, $i, $nuc);
 
853
                        $self->_print(sprintf("     $blocks%9d\n", $i + $nuc)) || return;
 
854
                }
 
855
 
 
856
                # Print the last line
 
857
                if (my $last = substr($str, $i)) {
 
858
                        my $last_len = length($last);
 
859
                        my $last_pat = 'a10' x int($last_len / 10) .'a'. $last_len % 10;
 
860
                        my $blocks = pack $out_pat,
 
861
                          unpack($last_pat, $last);
 
862
                        $self->_print(sprintf("     $blocks%9d\n", $length)) ||
 
863
                          return; # Add the length to the end
 
864
                }
 
865
 
 
866
                $self->_print( "//\n") || return;
 
867
 
 
868
                $self->flush if $self->_flush_on_write && defined $self->_fh;
 
869
    }
684
870
        return 1;
685
 
    }
686
871
}
687
872
 
688
873
=head2 _print_EMBL_FTHelper
690
875
 Title   : _print_EMBL_FTHelper
691
876
 Usage   :
692
877
 Function: Internal function
693
 
 Returns : 
 
878
 Returns : 1 if writing suceeded, otherwise undef
694
879
 Args    :
695
880
 
696
881
 
698
883
 
699
884
sub _print_EMBL_FTHelper {
700
885
   my ($self,$fth) = @_;
701
 
   
 
886
 
702
887
   if( ! ref $fth || ! $fth->isa('Bio::SeqIO::FTHelper') ) {
703
888
       $fth->warn("$fth is not a FTHelper class. Attempting to print, but there could be tears!");
704
889
   }
705
 
   
 
890
 
706
891
 
707
892
   #$self->_print( "FH   Key             Location/Qualifiers\n");
708
893
   #$self->_print( sprintf("FT   %-15s  %s\n",$fth->key,$fth->loc));
709
894
   # let
710
895
   if( $fth->key eq 'CONTIG' ) {
711
 
       $self->_print("XX\n");
 
896
       $self->_print("XX\n") || return;
712
897
       $self->_write_line_EMBL_regex("CO   ",
713
898
                                     "CO   ",$fth->loc,
714
 
                                     '\,|$',80); #'
715
 
       return;
716
 
   } 
 
899
                                     '\,|$',80) || return; #'
 
900
       return 1;
 
901
   }
717
902
   $self->_write_line_EMBL_regex(sprintf("FT   %-15s ",$fth->key),
718
903
                                 "FT                   ",$fth->loc,
719
 
                                 '\,|$',80); #'
 
904
                                 '\,|$',80) || return; #'
720
905
   foreach my $tag ( keys %{$fth->field} ) {
721
 
       if( ! defined $fth->field->{$tag} ) { next; } 
 
906
       if( ! defined $fth->field->{$tag} ) { next; }
722
907
       foreach my $value ( @{$fth->field->{$tag}} ) {
723
908
           $value =~ s/\"/\"\"/g;
724
909
           if ($value eq "_no_value") {
725
910
               $self->_write_line_EMBL_regex("FT                   ",
726
911
                                             "FT                   ",
727
 
                                             "/$tag",'.|$',80); #'
 
912
                                             "/$tag",'.|$',80) || return; #'
728
913
           }
729
914
           # there are almost 3x more quoted qualifier values and they
730
915
           # are more common too so we take quoted ones first
732
917
              my $pat = $value =~ /\s/ ? '\s|\-|$' : '.|\-|$';
733
918
              $self->_write_line_EMBL_regex("FT                   ",
734
919
                                            "FT                   ",
735
 
                                            "/$tag=\"$value\"",$pat,80);
 
920
                                            "/$tag=\"$value\"",$pat,80) || return;
736
921
           } else {
737
922
              $self->_write_line_EMBL_regex("FT                   ",
738
923
                                            "FT                   ",
739
 
                                            "/$tag=$value",'.|$',80); #'
 
924
                                            "/$tag=$value",'.|$',80) || return; #'
740
925
           }
741
926
       }
742
927
   }
 
928
 
 
929
   return 1;
743
930
}
744
931
 
745
932
#'
749
936
 Usage   :
750
937
 Function: Reads references from EMBL format. Internal function really
751
938
 Example :
752
 
 Returns : 
 
939
 Returns :
753
940
 Args    :
754
941
 
755
942
 
758
945
sub _read_EMBL_References {
759
946
   my ($self,$buffer) = @_;
760
947
   my (@refs);
761
 
   
 
948
 
762
949
   # assumme things are starting with RN
763
950
 
764
951
   if( $$buffer !~ /^RN/ ) {
777
964
       /^R/ || last;
778
965
       /^RP   (\d+)-(\d+)/ && do {$b1=$1;$b2=$2;};
779
966
       /^RX   MEDLINE;\s+(\d+)/ && do {$med=$1};
780
 
       /^RX   PUBMED;\s+(\d+)/ && do {$pm=$1};       
 
967
       /^RX   PUBMED;\s+(\d+)/ && do {$pm=$1};
781
968
       /^RA   (.*)/ && do {
782
969
           $au = $self->_concatenate_lines($au,$1); next;
783
970
       };
791
978
           $com = $self->_concatenate_lines($com,$1); next;
792
979
       };
793
980
   }
794
 
   
 
981
 
795
982
   my $ref = new Bio::Annotation::Reference;
796
983
   $au =~ s/;\s*$//g;
797
984
   $title =~ s/;\s*$//g;
807
994
 
808
995
   push(@refs,$ref);
809
996
   $$buffer = $_;
810
 
   
 
997
 
811
998
   return @refs;
812
999
}
813
1000
 
819
1006
           lines.
820
1007
 Example :
821
1008
 Returns : A Bio::Species object
822
 
 Args    : a reference to the current line buffer
 
1009
 Args    : a reference to the current line buffer, accession number
823
1010
 
824
1011
=cut
825
1012
 
826
1013
sub _read_EMBL_Species {
827
 
    my( $self, $buffer ) = @_;
 
1014
    my( $self, $buffer, $acc ) = @_;
828
1015
    my $org;
829
1016
 
830
1017
    $_ = $$buffer;
831
 
    my( $sub_species, $species, $genus, $common, @class, $ns_name );
 
1018
    my( $sub_species, $species, $genus, $common, $sci_name, $class_lines );
832
1019
    while (defined( $_ ||= $self->_readline )) {
833
 
 
834
 
        if (/^OS\s+((\S+)(?:\s+([^\(]\S*))?(?:\s+([^\(]\S*))?(?:\s+\((.*)\))?)/) {
835
 
            $ns_name = $1;
836
 
            $genus   = $2;
837
 
            $species = $3 || 'sp.';
838
 
            $sub_species = $4 if $4;
839
 
            $common      = $5 if $5;
840
 
        }
841
 
        elsif (s/^OC\s+//) {
842
 
            # only split on ';' or '.' so that 
843
 
            # classification that is 2 words will 
844
 
            # still get matched
845
 
            # use map to remove trailing/leading spaces
846
 
            chomp;
847
 
            push(@class,  map { s/^\s+//; s/\s+$//; $_; } split /[;\.]+/);
848
 
        }
849
 
        elsif (/^OG\s+(.*)/) {
850
 
            $org = $1;
851
 
        }
 
1020
        if (/^OS\s+(.+)/) {
 
1021
            $sci_name .= ($sci_name) ? ' '.$1 : $1;
 
1022
        }
 
1023
        elsif (s/^OC\s+(.+)$//) {
 
1024
            $class_lines .= $1;
 
1025
        }
 
1026
        elsif (/^OG\s+(.*)/) {
 
1027
            $org = $1;
 
1028
        }
852
1029
        else {
853
1030
            last;
854
1031
        }
858
1035
 
859
1036
    $$buffer = $_;
860
1037
 
 
1038
    $sci_name || return;
 
1039
 
 
1040
    # Convert data in classification lines into classification array.
 
1041
    # only split on ';' or '.' so that classification that is 2 or more words
 
1042
    # will still get matched, use map() to remove trailing/leading/intervening
 
1043
    # spaces
 
1044
    my @class = map { s/^\s+//; s/\s+$//; s/\s{2,}/ /g; $_; } split /[;\.]+/, $class_lines;
 
1045
 
 
1046
    # do we have a genus?
 
1047
    my $possible_genus = $class[-1];
 
1048
    $possible_genus .= "|$class[-2]" if $class[-2];
 
1049
    if ($sci_name =~ /^($possible_genus)/) {
 
1050
        $genus = $1;
 
1051
        ($species) = $sci_name =~ /^$genus\s+(.+)/;
 
1052
    }
 
1053
    else {
 
1054
        $species = $sci_name;
 
1055
    }
 
1056
 
861
1057
    # Don't make a species object if it is "Unknown" or "None"
862
 
    return if $genus =~ /^(Unknown|None)$/i;
 
1058
    if ($genus) {
 
1059
        return if $genus =~ /^(Unknown|None)$/i;
 
1060
    }
 
1061
 
 
1062
    # is this organism of rank species or is it lower?
 
1063
    # (doesn't catch everything, but at least the guess isn't dangerous)
 
1064
    if ($species =~ /subsp\.|var\./) {
 
1065
        ($species, $sub_species) = $species =~ /(.+)\s+((?:subsp\.|var\.).+)/;
 
1066
    }
 
1067
 
 
1068
    # sometimes things have common name in brackets, like
 
1069
    # Schizosaccharomyces pombe (fission yeast), so get rid of the common
 
1070
    # name bit. Probably dangerous if real scientific species name ends in
 
1071
    # bracketed bit.
 
1072
    unless ($class[-1] eq 'Viruses') {
 
1073
        ($species, $common) = $species =~ /^(.+)\s+\((.+)\)$/;
 
1074
        $sci_name =~ s/\s+\(.+\)$// if $common;
 
1075
    }
863
1076
 
864
1077
    # Bio::Species array needs array in Species -> Kingdom direction
865
 
    if ($class[0] eq 'Viruses') {
866
 
        push( @class, $ns_name );
867
 
    }
868
 
    elsif ($class[$#class] eq $genus) {
869
 
        push( @class, $species );
870
 
    } else {
871
 
        push( @class, $genus, $species );
 
1078
    unless ($class[-1] eq $sci_name) {
 
1079
        push(@class, $sci_name);
872
1080
    }
873
1081
    @class = reverse @class;
874
1082
 
 
1083
    # do minimal sanity checks before we hand off to Bio::Species which won't
 
1084
    # be able to give informative throw messages if it has to throw because
 
1085
    # of problems here
 
1086
    $self->throw("$acc seems to be missing its OS line: invalid.") unless $sci_name;
 
1087
    my %names;
 
1088
    foreach my $i (0..$#class) {
 
1089
        my $name = $class[$i];
 
1090
        $names{$name}++;
 
1091
        if ($names{$name} > 1 && $name ne $class[$i - 1]) {
 
1092
            $self->throw("$acc seems to have an invalid species classification.");
 
1093
        }
 
1094
    }
 
1095
 
875
1096
    my $make = Bio::Species->new();
876
 
    $make->classification( \@class, "FORCE" );  # no name validation please
877
 
    $make->common_name( $common      ) if $common;
 
1097
    $make->scientific_name($sci_name);
 
1098
    $make->classification(@class);
878
1099
    unless ($class[-1] eq 'Viruses') {
879
 
        $make->sub_species( $sub_species ) if $sub_species;
 
1100
        $make->genus($genus) if $genus;
 
1101
        $make->species($species) if $species;
 
1102
        $make->sub_species($sub_species) if $sub_species;
 
1103
        $make->common_name($common) if $common;
880
1104
    }
881
 
    $make->organelle  ( $org         ) if $org;
 
1105
    $make->organelle($org) if $org;
882
1106
    return $make;
883
1107
}
884
1108
 
899
1123
 
900
1124
    $_ = $$buffer;
901
1125
    while (defined( $_ ||= $self->_readline )) {
902
 
        
903
 
        if (my($databse, $prim_id, $sec_id)
904
 
                = /^DR   ([^\s;]+);\s*([^\s;]+);\s*([^\s;]+)?\.$/) {
905
 
            my $link = Bio::Annotation::DBLink->new();
906
 
            $link->database   ( $databse );
907
 
            $link->primary_id ( $prim_id );
908
 
            $link->optional_id( $sec_id  ) if $sec_id;
 
1126
        if( /^DR   ([^\s;]+);\s*([^\s;]+);?\s*([^\s;]+)?\.$/) {
 
1127
            my ($databse, $prim_id, $sec_id) = ($1,$2,$3);
 
1128
            my $link = Bio::Annotation::DBLink->new(-database    => $databse,
 
1129
                                                    -primary_id  => $prim_id,
 
1130
                                                    -optional_id => $sec_id);
 
1131
 
909
1132
            push(@db_link, $link);
910
 
        }
911
 
        else {
 
1133
        } else {
912
1134
            last;
913
1135
        }
914
 
        
915
 
        $_ = undef; # Empty $_ to trigger read of next line
 
1136
        $_ = undef;            # Empty $_ to trigger read of next line
916
1137
    }
917
1138
    
918
1139
    $$buffer = $_;
919
 
    
 
1140
 
920
1141
    return @db_link;
921
1142
}
922
1143
 
924
1145
 
925
1146
 Title   : _filehandle
926
1147
 Usage   : $obj->_filehandle($newval)
927
 
 Function: 
928
 
 Example : 
 
1148
 Function:
 
1149
 Example :
929
1150
 Returns : value of _filehandle
930
1151
 Args    : newvalue (optional)
931
1152
 
947
1168
 Usage   : _read_FTHelper_EMBL($buffer)
948
1169
 Function: reads the next FT key line
949
1170
 Example :
950
 
 Returns : Bio::SeqIO::FTHelper object 
 
1171
 Returns : Bio::SeqIO::FTHelper object
951
1172
 Args    : filehandle and reference to a scalar
952
1173
 
953
1174
 
955
1176
 
956
1177
sub _read_FTHelper_EMBL {
957
1178
    my ($self,$buffer) = @_;
958
 
    
 
1179
 
959
1180
    my ($key,   # The key of the feature
960
1181
        $loc,   # The location line from the feature
961
1182
        @qual,  # An arrray of lines making up the qualifiers
962
1183
        );
963
 
    
964
 
    if ($$buffer =~ /^FT\s{3}(\S+)\s+(\S+)/) {
 
1184
 
 
1185
    if ($$buffer =~ /^FT\s{3}(\S+)\s+(\S+)/ ) {
965
1186
        $key = $1;
966
1187
        $loc = $2;
967
1188
        # Read all the lines up to the next feature
1007
1228
        # No feature key
1008
1229
        return;
1009
1230
    }
1010
 
    
 
1231
 
1011
1232
    # Put the first line of the next feature into the buffer
1012
1233
    $$buffer = $_;
1013
1234
 
1029
1250
            if (substr($value, 0, 1) eq '"') {
1030
1251
                # Keep adding to value until we find the trailing quote
1031
1252
                # and the quotes are balanced
 
1253
                QUOTES:
1032
1254
                while ($value !~ /"$/ or $value =~ tr/"/"/ % 2) { #"
1033
1255
                    $i++;
1034
1256
                    my $next = $qual[$i];
1035
 
                    unless (defined($next)) {
1036
 
                        warn("Unbalanced quote in:\n", map("$_\n", @qual),
1037
 
                            "No further qualifiers will be added for this feature");
1038
 
                        last QUAL;
 
1257
                    if (!defined($next)) {
 
1258
                        $self->warn("Unbalanced quote in:\n".join("\n", @qual).
 
1259
                            "\nAdding quote to close...".
 
1260
                            "Check sequence quality!");
 
1261
                        $value .= '"';
 
1262
                        last QUOTES;
1039
1263
                    }
1040
1264
 
1041
1265
                    # Join to value with space if value or next line contains a space
1053
1277
        # Store the qualifier
1054
1278
        $out->field->{$qualifier} ||= [];
1055
1279
        push(@{$out->field->{$qualifier}},$value);
1056
 
    }   
 
1280
    }
1057
1281
 
1058
1282
    return $out;
1059
1283
}
1064
1288
 Usage   :
1065
1289
 Function: internal function
1066
1290
 Example :
1067
 
 Returns : 
 
1291
 Returns : 1 if writing suceeded, else undef
1068
1292
 Args    :
1069
1293
 
1070
1294
 
1073
1297
sub _write_line_EMBL {
1074
1298
   my ($self,$pre1,$pre2,$line,$length) = @_;
1075
1299
 
1076
 
   $length || die "Miscalled write_line_EMBL without length. Programming error!";
 
1300
   $length || $self->throw("Miscalled write_line_EMBL without length. Programming error!");
1077
1301
   my $subl = $length - length $pre2;
1078
1302
   my $linel = length $line;
1079
1303
   my $i;
1080
1304
 
1081
1305
   my $sub = substr($line,0,$length - length $pre1);
1082
1306
 
1083
 
   $self->_print( "$pre1$sub\n");
1084
 
   
 
1307
   $self->_print( "$pre1$sub\n") || return;
 
1308
 
1085
1309
   for($i= ($length - length $pre1);$i < $linel;) {
1086
1310
       $sub = substr($line,$i,($subl));
1087
 
       $self->_print( "$pre2$sub\n");
 
1311
       $self->_print( "$pre2$sub\n") || return;
1088
1312
       $i += $subl;
1089
1313
   }
1090
1314
 
 
1315
   return 1;
1091
1316
}
1092
1317
 
1093
1318
=head2 _write_line_EMBL_regex
1095
1320
 Title   : _write_line_EMBL_regex
1096
1321
 Usage   :
1097
1322
 Function: internal function for writing lines of specified
1098
 
           length, with different first and the next line 
 
1323
           length, with different first and the next line
1099
1324
           left hand headers and split at specific points in the
1100
1325
           text
1101
1326
 Example :
1110
1335
 
1111
1336
    #print STDOUT "Going to print with $line!\n";
1112
1337
 
1113
 
    $length || die "Programming error - called write_line_EMBL_regex without length.";
1114
 
 
1115
 
#    if( length $pre1 != length $pre2 ) {
1116
 
#      $self->throw(<<END);
1117
 
#Programming error - called write_line_EMBL_regex with different length pre1 and pre2 tags!
1118
 
#pre1 = $pre1
1119
 
#pre2 = $pre2
1120
 
#END
1121
 
#    }
 
1338
    $length || $self->throw("Programming error - called write_line_EMBL_regex without length.");
1122
1339
 
1123
1340
    my $subl = $length - (length $pre1) -1 ;
1124
 
 
1125
 
 
1126
 
 
1127
1341
    my( @lines );
1128
 
    while(defined $line && 
1129
 
          $line =~ m/(.{1,$subl})($regex)/g) {
1130
 
        push(@lines, $1.$2);
 
1342
 
 
1343
  CHUNK: while($line) {
 
1344
        foreach my $pat ($regex, '[,;\.\/-]\s|'.$regex, '[,;\.\/-]|'.$regex) {
 
1345
            if($line =~ m/^(.{1,$subl})($pat)(.*)/ ) {
 
1346
                my $l = $1.$2;
 
1347
                my $newl = $3;
 
1348
                $line = substr($line,length($l));
 
1349
                # be strict about not padding spaces according to
 
1350
                # genbank format
 
1351
                $l =~ s/\s+$//;
 
1352
                push(@lines, $l);
 
1353
                next CHUNK;
 
1354
            }
 
1355
        }
 
1356
        # if we get here none of the patterns matched $subl or less chars
 
1357
        $self->warn("trouble dissecting \"$line\"\n     into chunks ".
 
1358
                    "of $subl chars or less - this tag won't print right");
 
1359
        # insert a space char to prevent infinite loops
 
1360
        $line = substr($line,0,$subl) . " " . substr($line,$subl);
1131
1361
    }
1132
 
    foreach (@lines) { s/\s+$//; }
1133
 
    
1134
 
    # Print first line
1135
 
    my $s = shift(@lines) || '';    
1136
 
    $self->_print( "$pre1$s\n");
1137
 
    
1138
 
    # Print the rest
 
1362
    my $s = shift @lines;
 
1363
    ($self->_print("$pre1$s\n") || return) if $s;
1139
1364
    foreach my $s ( @lines ) {
1140
 
        $s = '' if( !defined $s );
1141
 
        $self->_print( "$pre2$s\n");
 
1365
        $self->_print("$pre2$s\n") || return;
1142
1366
    }
 
1367
 
 
1368
    return 1;
1143
1369
}
1144
1370
 
1145
1371
=head2 _post_sort
1146
1372
 
1147
1373
 Title   : _post_sort
1148
1374
 Usage   : $obj->_post_sort($newval)
1149
 
 Function: 
 
1375
 Function:
1150
1376
 Returns : value of _post_sort
1151
1377
 Args    : newvalue (optional)
1152
1378
 
1167
1393
 
1168
1394
 Title   : _show_dna
1169
1395
 Usage   : $obj->_show_dna($newval)
1170
 
 Function: 
 
1396
 Function:
1171
1397
 Returns : value of _show_dna
1172
1398
 Args    : newvalue (optional)
1173
1399
 
1188
1414
 
1189
1415
 Title   : _id_generation_func
1190
1416
 Usage   : $obj->_id_generation_func($newval)
1191
 
 Function: 
 
1417
 Function:
1192
1418
 Returns : value of _id_generation_func
1193
1419
 Args    : newvalue (optional)
1194
1420
 
1209
1435
 
1210
1436
 Title   : _ac_generation_func
1211
1437
 Usage   : $obj->_ac_generation_func($newval)
1212
 
 Function: 
 
1438
 Function:
1213
1439
 Returns : value of _ac_generation_func
1214
1440
 Args    : newvalue (optional)
1215
1441
 
1230
1456
 
1231
1457
 Title   : _sv_generation_func
1232
1458
 Usage   : $obj->_sv_generation_func($newval)
1233
 
 Function: 
 
1459
 Function:
1234
1460
 Returns : value of _sv_generation_func
1235
1461
 Args    : newvalue (optional)
1236
1462
 
1251
1477
 
1252
1478
 Title   : _kw_generation_func
1253
1479
 Usage   : $obj->_kw_generation_func($newval)
1254
 
 Function: 
 
1480
 Function:
1255
1481
 Returns : value of _kw_generation_func
1256
1482
 Args    : newvalue (optional)
1257
1483