← Back to branch summary

~ubuntu-branches/ubuntu/lucid/bioperl/lucid

« back to all changes in this revision

Viewing changes to Bio/SeqIO/swiss.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/swiss.pm
1
 
# $Id: swiss.pm,v 1.48 2002/02/14 18:04:22 jason Exp $
 
1
# $Id: swiss.pm,v 1.75 2003/12/22 18:33:15 heikki Exp $
2
2
#
3
3
# BioPerl module for Bio::SeqIO::swiss
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
#
30
30
This object can transform Bio::Seq objects to and from swissprot flat
31
31
file databases.
32
32
 
33
 
There is alot of flexibility here about how to dump things which I need
 
33
There is a lot of flexibility here about how to dump things which I need
34
34
to document fully.
35
35
 
36
36
 
56
56
To generate the ID line. If it is not there, it generates a sensible ID
57
57
line using a number of tools.
58
58
 
 
59
If you want to output annotations in swissprot format they need to be
 
60
stored in a Bio::Annotation::Collection object which is accessible
 
61
through the Bio::SeqI interface method L<annotation()|annotation>.  
 
62
 
 
63
The following are the names of the keys which are polled from a
 
64
L<Bio::Annotation::Collection> object.
 
65
 
 
66
reference       - Should contain Bio::Annotation::Reference objects
 
67
comment         - Should contain Bio::Annotation::Comment objects
 
68
dblink          - Should contain Bio::Annotation::DBLink objects
 
69
gene_name       - Should contain Bio::Annotation::SimpleValue object
59
70
 
60
71
=back
61
72
 
78
89
 Bug reports can be submitted via email or the web:
79
90
 
80
91
  bioperl-bugs@bio.perl.org
81
 
  http://bio.perl.org/bioperl-bugs/
 
92
  http://bugzilla.bioperl.org/
82
93
 
83
94
=head1 AUTHOR - Elia Stupka
84
95
 
85
 
Email elia@ebi.ac.uk
 
96
Email elia@tll.org.sg
86
97
 
87
98
Describe contact details here
88
99
 
99
110
package Bio::SeqIO::swiss;
100
111
use vars qw(@ISA);
101
112
use strict;
102
 
use Bio::Seq::RichSeq;
103
113
use Bio::SeqIO;
104
114
use Bio::SeqIO::FTHelper;
105
115
use Bio::SeqFeature::Generic;
106
116
use Bio::Species;
107
117
use Bio::Tools::SeqStats;
 
118
use Bio::Seq::SeqFactory;
 
119
use Bio::Annotation::Collection;
 
120
use Bio::Annotation::Comment;
 
121
use Bio::Annotation::Reference;
 
122
use Bio::Annotation::DBLink;
 
123
use Bio::Annotation::SimpleValue;
 
124
use Bio::Annotation::StructuredValue;
108
125
 
109
126
@ISA = qw(Bio::SeqIO);
110
127
 
116
133
  # hash for functions for decoding keys.
117
134
  $self->{'_func_ftunit_hash'} = {};
118
135
  $self->_show_dna(1); # sets this to one by default. People can change it
 
136
  if( ! defined $self->sequence_factory ) {
 
137
      $self->sequence_factory(new Bio::Seq::SeqFactory
 
138
                              (-verbose => $self->verbose(), 
 
139
                               -type => 'Bio::Seq::RichSeq'));      
 
140
  }
119
141
}
120
142
 
121
143
=head2 next_seq
132
154
sub next_seq {
133
155
   my ($self,@args) = @_;
134
156
   my ($pseq,$c,$line,$name,$desc,$acc,$seqc,$mol,$div,
135
 
       $date,$comment,@date_arr, @sec);
136
 
   my ($keywords,$acc_string);
137
 
   my $seq = Bio::Seq::RichSeq->new(-verbose =>$self->verbose());
 
157
       $date,$comment,@date_arr);
 
158
   
 
159
   my $genename = "";
 
160
   my ($annotation, %params, @features) = ( new Bio::Annotation::Collection);
 
161
 
138
162
   $line = $self->_readline;
139
163
 
140
164
   if( !defined $line) {
141
165
       return undef; # no throws - end of file
142
166
   }
143
 
   
 
167
 
144
168
   if( $line =~ /^\s+$/ ) {
145
169
       while( defined ($line = $self->_readline) ) {
146
170
           $line =~ /\S/ && last;
152
176
   
153
177
   # fixed to allow _DIVISION to be optional for bug #946
154
178
   # see bug report for more information
155
 
   $line =~ /^ID\s+([^\s_]+)(_([^\s_]+))?\s+([^\s;]+);\s+([^\s;]+);/ 
156
 
     || $self->throw("swissprot stream with no ID. Not swissprot in my book");
 
179
   if( !($line =~ /^ID\s+([^\s_]+)(_([^\s_]+))?\s+([^\s;]+);\s+([^\s;]+);/)
 
180
       && !($line =~ /^ID\s+([^\s_]+)(_([^\s_]+))?/ ) )  {
 
181
       $self->throw("swissprot stream with no ID. Not swissprot in my book");
 
182
   }
 
183
 
157
184
   if( $3 ) {
158
185
       $name = "$1$2";
159
 
       $seq->division($3);
 
186
       $params{'-division'} = $3;
160
187
   } else {
161
188
       $name = $1;
162
 
       $seq->division('UNK');       
 
189
       $params{'-division'} = 'UNK';
 
190
       $params{'-primary_id'} = $1;
163
191
   }
164
 
   
165
 
   $seq->primary_id($1);
166
 
   $seq->alphabet('protein');
 
192
   $params{'-alphabet'} = 'protein';
167
193
    # this is important to have the id for display in e.g. FTHelper, otherwise
168
194
    # you won't know which entry caused an error
169
 
   $seq->display_id($name);
 
195
   $params{'-display_id'} = $name;
170
196
   
171
197
   my $buffer = $line;
172
198
 
184
210
           $desc .= $desc ? " $1" : $1;
185
211
       }
186
212
       #Gene name
187
 
       elsif(/^GN\s+([^\.]+)/) {
188
 
           # Drop trailing spaces and dots
189
 
           s/[\. ]*$//;  # imported from swiss knife by cjm
190
 
           if (/^GN\s+(.*)/) {
191
 
               foreach my $gn (split(/ OR /,$1)) {
192
 
                   my $sim = Bio::Annotation::SimpleValue->new();
193
 
                   $sim->value($gn);
194
 
                   $seq->annotation->add_Annotation('gene_name',$sim);
 
213
       elsif(/^GN\s+(.*)/) {
 
214
           $genename .= " " if $genename;
 
215
           $genename .= $1;
 
216
           # has GN terminated yet?
 
217
           if($genename =~ s/[\. ]+$//) {
 
218
               my $gn = Bio::Annotation::StructuredValue->new();
 
219
               foreach my $gene (split(/ AND /, $genename)) {
 
220
                   $gene =~ s/^\(//;
 
221
                   $gene =~ s/\)$//;
 
222
                   $gn->add_value([-1,-1], split(/ OR /, $gene));
195
223
               }
 
224
               $annotation->add_Annotation('gene_name',$gn,
 
225
                                           "Bio::Annotation::SimpleValue");
196
226
           }
197
227
       }
198
228
       #accession number(s)
199
229
       elsif( /^AC\s+(.+)/) {
200
 
           $acc_string .= $acc_string ? " $1" : $1;
 
230
           my @accs = split(/[; ]+/, $1); # allow space in addition
 
231
           $params{'-accession_number'} = shift @accs 
 
232
               unless defined $params{'-accession_number'};
 
233
           push @{$params{'-secondary_accessions'}}, @accs;
201
234
       }
202
235
       #version number
203
236
       elsif( /^SV\s+(\S+);?/ ) {
204
237
           my $sv = $1;
205
238
           $sv =~ s/\;//;
206
 
           $seq->seq_version($sv);
 
239
           $params{'-seq_version'} = $sv;
207
240
       }
208
241
       #date
209
242
       elsif( /^DT\s+(.*)/ ) {
210
243
           my $date = $1;
211
244
           $date =~ s/\;//;
212
245
           $date =~ s/\s+$//;
213
 
           $seq->add_date($date);
 
246
           push @{$params{'-dates'}}, $date;
214
247
       }
215
248
       # Organism name and phylogenetic information
216
249
       elsif (/^O[SCG]/) {
217
250
           my $species = $self->_read_swissprot_Species(\$buffer);
218
 
           $seq->species( $species );
 
251
           $params{'-species'}= $species;
219
252
           # now we are one line ahead -- so continue without reading the next
220
253
           # line   HL 05/11/2000
221
254
           next;
225
258
           my $refs = $self->_read_swissprot_References(\$buffer);
226
259
 
227
260
           foreach my $r (@$refs) {
228
 
               $seq->annotation->add_Annotation('reference',$r);
 
261
               $annotation->add_Annotation('reference',$r);
229
262
           }
230
263
           # now we are one line ahead -- so continue without reading the next
231
264
           # line   HL 05/11/2000
248
281
           # note: don't try to process comments here -- they may contain
249
282
           # structure. LP 07/30/2000
250
283
           $commobj->text($comment);
251
 
           $seq->annotation->add_Annotation('comment',$commobj);
 
284
           $annotation->add_Annotation('comment',$commobj);
252
285
           $comment = "";
253
286
           # now we are one line ahead -- so continue without reading the next
254
287
           # line   HL 05/11/2000
269
302
                   $dblinkobj->comment($comment);
270
303
               }
271
304
           }
272
 
           $seq->annotation->add_Annotation('dblink',$dblinkobj);
 
305
           $annotation->add_Annotation('dblink',$dblinkobj);
273
306
       }
274
307
       #keywords
275
308
       elsif( /^KW\s+(.*)$/ ) {
276
 
           $keywords .= $keywords ? " $1" : $1;
 
309
           my @kw = split(/\s*\;\s*/,$1);
 
310
           defined $kw[-1] && $kw[-1] =~ s/\.$//;
 
311
           push @{$params{'-keywords'}}, @kw;   
277
312
       }
278
313
 
279
314
       # Get next line. Getting here assumes that we indeed need to read the
292
327
       # process ftunit
293
328
       # when parsing of the line fails we get undef returned
294
329
       if($ftunit) {
295
 
           $ftunit->_generic_seqfeature($seq, "SwissProt");
 
330
           push(@features,
 
331
                $ftunit->_generic_seqfeature($self->location_factory(),
 
332
                                             $params{'-seqid'}, "SwissProt"));
296
333
       } else {
297
334
           $self->warn("failed to parse feature table line for seq " .
298
 
                       $seq->display_id());
 
335
                       $params{'-display_id'});
299
336
       }
300
337
   }
301
338
   if( $buffer !~ /^SQ/  ) {
310
347
       s/[^A-Za-z]//g;
311
348
       $seqc .= $_;
312
349
   }
313
 
   $seq->seq($seqc);
314
 
   $seq->desc($desc);
315
 
   $seq->keywords($keywords);
316
 
   $acc_string =~ s/\;\s*/ /g;
317
 
   ( $acc, @sec ) = split " ",$acc_string;
318
 
   $seq->accession_number($acc);
319
 
   foreach my $s (@sec) {
320
 
     $seq->add_secondary_accession($s);
321
 
   }
 
350
 
 
351
   my $seq=  $self->sequence_factory->create
 
352
       (-verbose  => $self->verbose,
 
353
        %params,
 
354
        -seq      => $seqc,
 
355
        -desc     => $desc,
 
356
        -features => \@features,
 
357
        -annotation => $annotation,
 
358
        );
 
359
 
 
360
   # The annotation doesn't get added by the contructor
 
361
   $seq->annotation($annotation);
322
362
 
323
363
   return $seq;
324
364
}
329
369
 Usage   : $stream->write_seq($seq)
330
370
 Function: writes the $seq object (must be seq) to the stream
331
371
 Returns : 1 for success and 0 for error
332
 
 Args    : Bio::Seq
 
372
 Args    : array of 1 to n Bio::SeqI objects
333
373
 
334
374
 
335
375
=cut
336
376
 
337
377
sub write_seq {
338
 
   my ($self,$seq) = @_;
339
 
 
340
 
   if( !defined $seq ) {
341
 
       $self->throw("Attempting to write with no seq!");
342
 
   }
343
 
 
344
 
   if( ! ref $seq || ! $seq->isa('Bio::SeqI') ) {
345
 
       $self->warn(" $seq is not a SeqI compliant module. Attempting to dump, but may fail!");
346
 
   }
347
 
 
348
 
   my $i;
349
 
   my $str = $seq->seq;
350
 
   
351
 
   my $mol;
352
 
   my $div;
353
 
   my $len = $seq->length();
354
 
 
355
 
   if ( !$seq->can('division') || ! defined ($div = $seq->division()) ) {
356
 
       $div = 'UNK';
357
 
   }
358
 
   
359
 
   if( ! $seq->can('alphabet') || ! defined ($mol = $seq->alphabet) ) {
360
 
       $mol = 'XXX';
361
 
   }
362
 
   
363
 
   my $temp_line;
364
 
   if( $self->_id_generation_func ) {
365
 
       $temp_line = &{$self->_id_generation_func}($seq);
366
 
   } else {
367
 
       #$temp_line = sprintf ("%10s     STANDARD;      %3s;   %d AA.",
368
 
       #                     $seq->primary_id()."_".$div,$mol,$len);
369
 
       # Reconstructing the ID relies heavily upon the input source having
370
 
       # been in a format that is parsed as this routine expects it -- that is,
371
 
       # by this module itself. This is bad, I think, and immediately breaks
372
 
       # if e.g. the Bio::DB::GenPept module is used as input.
373
 
       # Hence, switch to display_id(); _every_ sequence is supposed to have
374
 
       # this. HL 2000/09/03
375
 
       $mol =~ s/protein/PRT/;
376
 
       $temp_line = sprintf ("%10s     STANDARD;      %3s;   %d AA.",
377
 
                             $seq->display_id(), $mol, $len);
378
 
   }
379
 
 
380
 
   $self->_print( "ID   $temp_line\n");
381
 
 
382
 
   # if there, write the accession line
383
 
   local($^W) = 0;   # supressing warnings about uninitialized fields
384
 
 
385
 
   if( $self->_ac_generation_func ) {
386
 
       $temp_line = &{$self->_ac_generation_func}($seq);
387
 
       $self->_print( "AC   $temp_line\n");
388
 
   } else {
389
 
       if ($seq->can('accession_number') ) {
390
 
           $self->_print("AC   ",$seq->accession_number,";");
391
 
           if ($seq->can('get_secondary_accessions') ) {
392
 
             foreach my $sacc ($seq->get_secondary_accessions) {
393
 
               $self->_print(" ",$sacc,";");
394
 
             }
395
 
             $self->_print("\n");
396
 
           }
397
 
           else {
398
 
               $self->_print("\n");
399
 
           }
400
 
       }
401
 
       # otherwise - cannot print <sigh>
402
 
   }
403
 
 
404
 
   # Date lines
405
 
 
406
 
   if( $seq->can('get_dates') ) {
407
 
       foreach my $dt ( $seq->get_dates() ) {
408
 
           $self->_write_line_swissprot_regex("DT   ","DT   ",
409
 
                                              $dt,"\\s\+\|\$",80);
410
 
       }
411
 
   }
412
 
 
413
 
   #Definition lines
414
 
   $self->_write_line_swissprot_regex("DE   ","DE   ",$seq->desc(),"\\s\+\|\$",80);
415
 
 
416
 
   #Gene name
417
 
   if ($seq->annotation->can('get_Annotations') &&
418
 
       (my @genes = $seq->annotation->get_Annotations('gene_name') ) ) {
419
 
       $self->_print("GN   ",join(' OR ', map { $_->value } @genes),".\n");
420
 
   }
421
 
   
422
 
   # Organism lines
423
 
   if ($seq->can('species') && (my $spec = $seq->species)) {
424
 
        my($species, @class) = $spec->classification();
425
 
        my $genus = $class[0];
426
 
        my $OS = "$genus $species";
427
 
        if (my $ssp = $spec->sub_species) {
428
 
            $OS .= " $ssp";
429
 
        }
430
 
        if (my $common = $spec->common_name) {
431
 
            $OS .= " ($common).";
432
 
        }
433
 
        if ($class[$#class] =~ /viruses/i) { # different OS / OC syntax
434
 
          $OS = $spec->common_name;          # for viruses LP 09/16/2000
435
 
          unshift @class, $species;
436
 
        }
437
 
        $self->_print( "OS   $OS\n");
438
 
        my $OC = join('; ', reverse(@class)) .'.';
439
 
        $self->_write_line_swissprot_regex("OC   ","OC   ",$OC,"\; \|\$",80);
440
 
        if ($spec->organelle) {
441
 
            $self->_write_line_swissprot_regex("OG   ","OG   ",$spec->organelle,"\; \|\$",80);
442
 
        }
443
 
        if ($spec->ncbi_taxid) {
444
 
            $self->_print("OX   NCBI_TaxID=".$spec->ncbi_taxid.";\n");
445
 
        }
446
 
   }
447
 
   
448
 
   # Reference lines
449
 
   my $t = 1;
450
 
   foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
451
 
       $self->_print( "RN   [$t]\n");
452
 
#       if ($ref->start && $ref->end) {
453
 
#          # changed by jason <jason@chg.mc.duke.edu> on 3/16/00 
454
 
#          print "RP   SEQUENCE OF ",$ref->start,"-",$ref->end," FROM N.A.\n";
455
 
#          # to
456
 
#          $self->_print("RP   SEQUENCE OF ",$ref->start,"-",$ref->end," FROM N.A.\n");
457
 
#       }
458
 
#       elsif ($ref->rp) {
459
 
#          # changed by jason <jason@chg.mc.duke.edu> on 3/16/00 
460
 
#          print "RP   ",$ref->rp,"\n";
461
 
#          # to
462
 
 
463
 
       # changed by lorenz 08/03/00
464
 
       # j.gilbert and h.lapp agreed that the rp line in swissprot seems 
465
 
       # more like a comment than a parseable value, so print it as is
466
 
       if ($ref->rp) {
467
 
         $self->_write_line_swissprot_regex("RP   ","RP   ",$ref->rp,
468
 
                                            "\\s\+\|\$",80);
469
 
       }
470
 
       if ($ref->comment) {
471
 
         $self->_write_line_swissprot_regex("RC   ","RC   ",$ref->comment,
472
 
                                            "\\s\+\|\$",80);
473
 
       }
474
 
       if ($ref->medline) {
475
 
         # new RX format in swissprot LP 09/17/00
476
 
         if ($ref->pubmed) {
477
 
         $self->_write_line_swissprot_regex("RX   ","RX   ",
478
 
                                            "MEDLINE=".$ref->medline.
479
 
                                            "; PubMed=".$ref->pubmed.";",
480
 
                                            "\\s\+\|\$",80);
481
 
         } else {
482
 
         $self->_write_line_swissprot_regex("RX   MEDLINE; ","RX   MEDLINE; ",
483
 
                                            $ref->medline.".","\\s\+\|\$",80);
484
 
         }
485
 
       }
486
 
       my $author = $ref->authors .';' if($ref->authors);
487
 
       my $title = $ref->title .';' if( $ref->title);
488
 
       
489
 
       $self->_write_line_swissprot_regex("RA   ","RA   ",$author,"\\s\+\|\$",80);
490
 
       $self->_write_line_swissprot_regex("RT   ","RT   ",$title,"\\s\+\|\$",80);
491
 
       $self->_write_line_swissprot_regex("RL   ","RL   ",$ref->location,"\\s\+\|\$",80);
492
 
       $t++;
493
 
   }
494
 
   
495
 
   # Comment lines
496
 
 
497
 
   foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
498
 
       foreach my $cline (split ("\n", $comment->text)) {
499
 
         while (length $cline > 74) {
500
 
           $self->_print("CC   ",(substr $cline,0,74),"\n");
501
 
           $cline = substr $cline,74;
502
 
         }
503
 
         $self->_print("CC   ",$cline,"\n");
504
 
       }
505
 
   }
506
 
 
507
 
   foreach my $dblink ( $seq->annotation->get_Annotations('dblink') ) 
508
 
   {
509
 
     if (defined($dblink->comment)&&($dblink->comment)) {
510
 
         $self->_print("DR   ",$dblink->database,"; ",$dblink->primary_id,"; ",
511
 
                       $dblink->optional_id,"; ",$dblink->comment,".\n");
512
 
     } elsif($dblink->optional_id) {
513
 
         $self->_print("DR   ",$dblink->database,"; ",$dblink->primary_id,"; ",
514
 
                       $dblink->optional_id,".\n");
515
 
     }
516
 
     else {
517
 
         $self->_print("DR   ",$dblink->database,"; ",$dblink->primary_id,"; ",
518
 
                       "-.\n");
519
 
     }
520
 
   }   
521
 
 
522
 
   # if there, write the kw line
523
 
   
524
 
   if( $self->_kw_generation_func ) {
525
 
       $temp_line = &{$self->_kw_generation_func}($seq);
526
 
       $self->_print( "KW   $temp_line\n");
527
 
   } else {
528
 
       if( $seq->can('keywords') ) {
529
 
         $self->_write_line_swissprot_regex("KW   ","KW   ",
530
 
                                            $seq->keywords,"\\s\+\|\$",80);       
531
 
       }
532
 
   }
533
 
 
534
 
#Check if there is seqfeatures before printing the FT line
535
 
   my @feats = $seq->top_SeqFeatures;
536
 
       if ($feats[0]) {
537
 
         if( defined $self->_post_sort ) {
538
 
             
539
 
             # we need to read things into an array. Process. Sort them. Print 'em
540
 
             
541
 
             my $post_sort_func = $self->_post_sort();
542
 
             my @fth;
543
 
             
544
 
           my @feats = $seq->top_SeqFeatures;
545
 
           
546
 
           foreach my $sf ( $seq->top_SeqFeatures ) {
547
 
               push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
548
 
           }
549
 
           @fth = sort { &$post_sort_func($a,$b) } @fth;
550
 
           
551
 
           foreach my $fth ( @fth ) {
552
 
               $self->_print_swissprot_FTHelper($fth);
553
 
           }
554
 
       } else {
555
 
           # not post sorted. And so we can print as we get them.
556
 
           # lower memory load...
557
 
           
558
 
           foreach my $sf ( $seq->top_SeqFeatures ) {
559
 
               my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
560
 
               foreach my $fth ( @fth ) {
561
 
                   if( ! $fth->isa('Bio::SeqIO::FTHelper') ) {
562
 
                       $sf->throw("Cannot process FTHelper... $fth");
563
 
                   }
564
 
                   
565
 
                   $self->_print_swissprot_FTHelper($fth);
566
 
               }
567
 
           }
568
 
       }
569
 
       
570
 
       if( $self->_show_dna() == 0 ) {
571
 
           return;
572
 
       }
573
 
   }
574
 
   # finished printing features.
575
 
 
576
 
   # molecular weight
577
 
   my $mw = ${Bio::Tools::SeqStats->get_mol_wt($seq->primary_seq)}[0];
578
 
   # checksum
579
 
   # was crc32 checksum, changed it to crc64 
580
 
   my $crc64 = $self->_crc64(\$str); 
581
 
   $self->_print( sprintf("SQ   SEQUENCE  %4d AA;  %d MW;  %16s CRC64;\n",
582
 
                          $len,$mw,$crc64));
583
 
   $self->_print( "     ");
584
 
   my $linepos;
585
 
   for ($i = 0; $i < length($str); $i += 10) {
586
 
       $self->_print( substr($str,$i,10), " ");
587
 
       $linepos += 11;
588
 
       if( ($i+10)%60 == 0 && (($i+10) < length($str))) {
589
 
           $self->_print( "\n     ");
590
 
      }
591
 
   }
592
 
   $self->_print( "\n//\n");
593
 
   return 1;
 
378
    my ($self,@seqs) = @_;
 
379
    foreach my $seq ( @seqs ) {
 
380
        $self->throw("Attempting to write with no seq!") unless defined $seq;
 
381
 
 
382
        if( ! ref $seq || ! $seq->isa('Bio::SeqI') ) {
 
383
            $self->warn(" $seq is not a SeqI compliant module. Attempting to dump, but may fail!");
 
384
        }
 
385
 
 
386
        my $i;
 
387
        my $str = $seq->seq;
 
388
 
 
389
        my $mol;
 
390
        my $div;
 
391
        my $len = $seq->length();
 
392
 
 
393
        if ( !$seq->can('division') || ! defined ($div = $seq->division()) ) {
 
394
            $div = 'UNK';
 
395
        }
 
396
 
 
397
        if( ! $seq->can('alphabet') || ! defined ($mol = $seq->alphabet) ) {
 
398
            $mol = 'XXX';
 
399
        }
 
400
 
 
401
        $self->warn("No whitespace allowed in SWISS-PROT display id [". $seq->display_id. "]")
 
402
            if $seq->display_id =~ /\s/;
 
403
 
 
404
        my $temp_line;
 
405
        if( $self->_id_generation_func ) {
 
406
            $temp_line = &{$self->_id_generation_func}($seq);
 
407
        } else {
 
408
            #$temp_line = sprintf ("%10s     STANDARD;      %3s;   %d AA.",
 
409
            #                $seq->primary_id()."_".$div,$mol,$len);
 
410
            # Reconstructing the ID relies heavily upon the input source having
 
411
            # been in a format that is parsed as this routine expects it -- that is,
 
412
            # by this module itself. This is bad, I think, and immediately breaks
 
413
            # if e.g. the Bio::DB::GenPept module is used as input.
 
414
            # Hence, switch to display_id(); _every_ sequence is supposed to have
 
415
            # this. HL 2000/09/03
 
416
            $mol =~ s/protein/PRT/;
 
417
            $temp_line = sprintf ("%-10s     STANDARD;      %3s;   %d AA.",
 
418
                                  $seq->display_id(), $mol, $len);
 
419
        }
 
420
 
 
421
        $self->_print( "ID   $temp_line\n");
 
422
 
 
423
        # if there, write the accession line
 
424
        local($^W) = 0; # supressing warnings about uninitialized fields
 
425
 
 
426
        if( $self->_ac_generation_func ) {
 
427
            $temp_line = &{$self->_ac_generation_func}($seq);
 
428
            $self->_print( "AC   $temp_line\n");
 
429
        } else {
 
430
            if ($seq->can('accession_number') ) {
 
431
                $self->_print("AC   ",$seq->accession_number,";");
 
432
                if ($seq->can('get_secondary_accessions') ) {
 
433
                    foreach my $sacc ($seq->get_secondary_accessions) {
 
434
                        $self->_print(" ",$sacc,";");
 
435
                    }
 
436
                    $self->_print("\n");
 
437
                }
 
438
                else {
 
439
                    $self->_print("\n");
 
440
                }
 
441
            }
 
442
            # otherwise - cannot print <sigh>
 
443
        }
 
444
 
 
445
        # Date lines
 
446
 
 
447
        if( $seq->can('get_dates') ) {
 
448
            foreach my $dt ( $seq->get_dates() ) {
 
449
                $self->_write_line_swissprot_regex("DT   ","DT   ",
 
450
                                                   $dt,"\\s\+\|\$",80);
 
451
            }
 
452
        }
 
453
 
 
454
        #Definition lines
 
455
        $self->_write_line_swissprot_regex("DE   ","DE   ",$seq->desc(),"\\s\+\|\$",80);
 
456
 
 
457
        #Gene name
 
458
        if ((my @genes = $seq->annotation->get_Annotations('gene_name') ) ) {
 
459
            $self->_print("GN   ",
 
460
                          join(' OR ',
 
461
                               map {
 
462
                                   $_->isa("Bio::Annotation::StructuredValue") ?
 
463
                                       $_->value(-joins => [" AND ", " OR "]) :
 
464
                                       $_->value();
 
465
                               } @genes),
 
466
                          ".\n");
 
467
        }
 
468
 
 
469
        # Organism lines
 
470
        if ($seq->can('species') && (my $spec = $seq->species)) {
 
471
            my($species, @class) = $spec->classification();
 
472
            my $genus = $class[0];
 
473
            my $OS = "$genus $species";
 
474
            if ($class[-1] eq 'Viruses') {
 
475
                $OS = $species;
 
476
                $OS .=  " ". $spec->sub_species if $spec->sub_species;
 
477
            } else {
 
478
                if ($class[$#class] =~ /viruses/i) {
 
479
                    # different OS / OC syntax for viruses LP 09/16/2000
 
480
                    shift @class;
 
481
                }
 
482
                if (my $ssp = $spec->sub_species) {
 
483
                    $OS .= " $ssp";
 
484
                }
 
485
                foreach (($spec->variant, $spec->common_name)) {
 
486
                    $OS .= " ($_)" if $_;
 
487
                }
 
488
            }
 
489
            $self->_print( "OS   $OS.\n");
 
490
            my $OC = join('; ', reverse(@class)) .'.';
 
491
            $self->_write_line_swissprot_regex("OC   ","OC   ",$OC,"\; \|\$",80);
 
492
            if ($spec->organelle) {
 
493
                $self->_write_line_swissprot_regex("OG   ","OG   ",$spec->organelle,"\; \|\$",80);
 
494
            }
 
495
            if ($spec->ncbi_taxid) {
 
496
                $self->_print("OX   NCBI_TaxID=".$spec->ncbi_taxid.";\n");
 
497
            }
 
498
        }
 
499
 
 
500
        # Reference lines
 
501
        my $t = 1;
 
502
        foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
 
503
            $self->_print( "RN   [$t]\n");
 
504
            # changed by lorenz 08/03/00
 
505
            # j.gilbert and h.lapp agreed that the rp line in swissprot seems 
 
506
            # more like a comment than a parseable value, so print it as is
 
507
            if ($ref->rp) {
 
508
                $self->_write_line_swissprot_regex("RP   ","RP   ",$ref->rp,
 
509
                                                   "\\s\+\|\$",80);
 
510
            }
 
511
            if ($ref->comment) {
 
512
                $self->_write_line_swissprot_regex("RC   ","RC   ",$ref->comment,
 
513
                                                   "\\s\+\|\$",80);
 
514
            }
 
515
            if ($ref->medline) {
 
516
                # new RX format in swissprot LP 09/17/00
 
517
                if ($ref->pubmed) {
 
518
                    $self->_write_line_swissprot_regex("RX   ","RX   ",
 
519
                                                       "MEDLINE=".$ref->medline.
 
520
                                                       "; PubMed=".$ref->pubmed.";",
 
521
                                                       "\\s\+\|\$",80);
 
522
                } else {
 
523
                    $self->_write_line_swissprot_regex("RX   MEDLINE; ","RX   MEDLINE; ",
 
524
                                                       $ref->medline.".","\\s\+\|\$",80);
 
525
                }
 
526
            }
 
527
            my $author = $ref->authors .';' if($ref->authors);
 
528
            my $title = $ref->title .';' if( $ref->title);
 
529
 
 
530
            $self->_write_line_swissprot_regex("RA   ","RA   ",$author,"\\s\+\|\$",80);
 
531
            $self->_write_line_swissprot_regex("RT   ","RT   ",$title,"\\s\+\|\$",80);
 
532
            $self->_write_line_swissprot_regex("RL   ","RL   ",$ref->location,"\\s\+\|\$",80);
 
533
            $t++;
 
534
        }
 
535
 
 
536
        # Comment lines
 
537
 
 
538
        foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
 
539
            foreach my $cline (split ("\n", $comment->text)) {
 
540
                while (length $cline > 74) {
 
541
                    $self->_print("CC   ",(substr $cline,0,74),"\n");
 
542
                    $cline = substr $cline,74;
 
543
                }
 
544
                $self->_print("CC   ",$cline,"\n");
 
545
            }
 
546
        }
 
547
 
 
548
        foreach my $dblink ( $seq->annotation->get_Annotations('dblink') ) 
 
549
        {
 
550
            if (defined($dblink->comment)&&($dblink->comment)) {
 
551
                $self->_print("DR   ",$dblink->database,"; ",$dblink->primary_id,"; ",
 
552
                              $dblink->optional_id,"; ",$dblink->comment,".\n");
 
553
            } elsif($dblink->optional_id) {
 
554
                $self->_print("DR   ",$dblink->database,"; ",
 
555
                              $dblink->primary_id,"; ",
 
556
                              $dblink->optional_id,".\n");
 
557
            }
 
558
            else {
 
559
                $self->_print("DR   ",$dblink->database,
 
560
                              "; ",$dblink->primary_id,"; ","-.\n");
 
561
            }
 
562
        }   
 
563
 
 
564
        # if there, write the kw line
 
565
        {
 
566
            my( $kw );
 
567
            if( my $func = $self->_kw_generation_func ) {
 
568
                $kw = &{$func}($seq);
 
569
            } elsif( $seq->can('keywords') ) {          
 
570
                $kw = $seq->keywords;
 
571
                if( ref($kw) =~ /ARRAY/i ) {
 
572
                    $kw = join("; ", @$kw);
 
573
                }
 
574
                $kw .= '.' if( $kw !~ /\.$/ );
 
575
            }       
 
576
            $self->_write_line_swissprot_regex("KW   ","KW   ",
 
577
                                               $kw, "\\s\+\|\$",80);           
 
578
        }
 
579
 
 
580
        #Check if there is seqfeatures before printing the FT line
 
581
        my @feats = $seq->can('top_SeqFeatures') ? $seq->top_SeqFeatures : ();
 
582
        if ($feats[0]) {
 
583
            if( defined $self->_post_sort ) {
 
584
 
 
585
                # we need to read things into an array. Process. Sort them. Print 'em
 
586
                
 
587
                my $post_sort_func = $self->_post_sort();
 
588
                my @fth;
 
589
 
 
590
                foreach my $sf ( @feats ) {
 
591
                    push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
 
592
                }
 
593
                @fth = sort { &$post_sort_func($a,$b) } @fth;
 
594
 
 
595
                foreach my $fth ( @fth ) {
 
596
                    $self->_print_swissprot_FTHelper($fth);
 
597
                }
 
598
            } else {
 
599
                # not post sorted. And so we can print as we get them.
 
600
                # lower memory load...
 
601
 
 
602
                foreach my $sf ( @feats ) {
 
603
                    my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
 
604
                    foreach my $fth ( @fth ) {
 
605
                        if( ! $fth->isa('Bio::SeqIO::FTHelper') ) {
 
606
                            $sf->throw("Cannot process FTHelper... $fth");
 
607
                        }
 
608
 
 
609
                        $self->_print_swissprot_FTHelper($fth);
 
610
                    }
 
611
                }
 
612
            }
 
613
 
 
614
            if( $self->_show_dna() == 0 ) {
 
615
                return;
 
616
            }
 
617
        }
 
618
        # finished printing features.
 
619
 
 
620
        # molecular weight
 
621
        my $mw = ${Bio::Tools::SeqStats->get_mol_wt($seq->primary_seq)}[0];
 
622
        # checksum
 
623
        # was crc32 checksum, changed it to crc64 
 
624
        my $crc64 = $self->_crc64(\$str); 
 
625
        $self->_print( sprintf("SQ   SEQUENCE  %4d AA;  %d MW;  %16s CRC64;\n",
 
626
                               $len,$mw,$crc64));
 
627
        $self->_print( "     ");
 
628
        my $linepos;
 
629
        for ($i = 0; $i < length($str); $i += 10) {
 
630
            $self->_print( substr($str,$i,10), " ");
 
631
            $linepos += 11;
 
632
            if( ($i+10)%60 == 0 && (($i+10) < length($str))) {
 
633
                $self->_print( "\n     ");
 
634
            }
 
635
        }
 
636
        $self->_print( "\n//\n");
 
637
 
 
638
        $self->flush if $self->_flush_on_write && defined $self->_fh;
 
639
        return 1;
 
640
    }
594
641
}
595
642
 
596
643
# Thanks to James Gilbert for the following two. LP 08/01/2000
738
785
                  "Attempting to print, but there could be tears!");
739
786
   }
740
787
 
741
 
   if( $fth->loc =~ /(\?|\d+)?\.\.(\?|\d+)?/ ) {
 
788
   if( $fth->loc =~ /(\?|\d+|\>\d+|<\d+)?\.\.(\?|\d+|<\d+|>\d+)?/ ) {
742
789
       $start = $1 if defined $1;
743
790
       $end = $2 if defined $2;
744
791
 
755
802
                                              substr($fth->key,0,8),
756
803
                                              $start,$end),
757
804
                                      "FT                                ",
758
 
                                      $desc,'\s+|$',80);
 
805
                                      $desc.'.','\s+|$',80);
759
806
}
760
807
#'
761
808
 
843
890
    my $org;
844
891
 
845
892
    $_ = $$buffer;
846
 
    my( $sub_species, $species, $genus, $common, @class, $osline, $ncbi_taxid );
 
893
    my( $subspecies, $species, $genus, $common, $variant, $ncbi_taxid, $ns_name );
 
894
    my @class;
 
895
    my ($binomial, $descr);
 
896
    my $osline = "";
 
897
 
847
898
    while (defined( $_ ||= $self->_readline )) {
848
 
        if (/^OS\s+((\S+)(?:\s+([^\(]\S*))?(?:\s+([^\(]\S*))?(?:\s+\((.*)\))?.*)/) {
849
 
            $osline = $1;
850
 
            $genus   = $2;
851
 
            if ($3) {
852
 
                $species = $3;
853
 
                # remove trailing dot -- TrEMBL has that. HL 05/11/2000
854
 
                # forgot to escape the dot. LP 07/30/2000
855
 
                $species =~ s/\.$//;
856
 
            } else {
857
 
                $species = "sp.";
 
899
        last unless /^O[SCGX]/;
 
900
        # believe it or not, but OS may come multiple times -- at this time
 
901
        # we can't capture multiple species
 
902
        if(/^OS\s+(\S.+)/ && (! defined($binomial))) {
 
903
            $osline .= " " if $osline;
 
904
            $osline .= $1;
 
905
            if($osline =~ s/(,|, and|\.)$//) {
 
906
                ($binomial, $descr) = $osline =~ /(\S[^\(]+)(.*)/;
 
907
                ($ns_name) = $binomial;
 
908
                $ns_name =~ s/\s+$//; #####
 
909
                ($genus, $species, $subspecies) = split(/\s+/, $binomial);
 
910
                $species = "sp." unless $species;
 
911
                while($descr =~ /\(([^\)]+)\)/g) {
 
912
                    my $item = $1;
 
913
                    # strain etc may not necessarily come first (yes, swissprot
 
914
                    # is messy)
 
915
                    if((! defined($variant)) &&
 
916
                       (($item =~ /(^|[^\(\w])([Ss]train|isolate|serogroup|serotype|subtype|clone)\b/) ||
 
917
                        ($item =~ /^(biovar|pv\.|type\s+)/))) {
 
918
                        $variant = $item;
 
919
                    } elsif($item =~ s/^subsp\.\s+//) {
 
920
                        if(! $subspecies) {
 
921
                            $subspecies = $item;
 
922
                        } elsif(! $variant) {
 
923
                            $variant = $item;
 
924
                        }
 
925
                    } elsif(! defined($common)) {
 
926
                        # we're only interested in the first common name
 
927
                        $common = $item;
 
928
                        if((index($common, '(') >= 0) &&
 
929
                           (index($common, ')') < 0)) {
 
930
                            $common .= ')';
 
931
                        }
 
932
                    }
 
933
                }
858
934
            }
859
 
            $sub_species = $4 if $4;
860
 
            $common      = $5 if $5;
861
935
        }
862
936
        elsif (s/^OC\s+//) {
863
937
            push(@class, split /[\;\.]\s*/);
864
 
            if ($class[0] =~ /viruses/i) { # viruses have different OS / OC syntax
865
 
              $common = $osline;           # LP 09/16/2000
 
938
            if($class[0] =~ /viruses/i) { 
 
939
                # viruses have different OS/OC syntax
 
940
                my @virusnames = split(/\s+/, $binomial);
 
941
                $species = (@virusnames > 1) ? pop(@virusnames) : '';
 
942
                $genus = join(" ", @virusnames);
 
943
                $subspecies = $descr;
866
944
            }
867
945
        }
868
946
        elsif (/^OG\s+(.*)/) {
869
947
            $org = $1;
870
948
        }
871
 
        elsif (/^OX\s+(.*)\;/) {
 
949
        elsif (/^OX\s+(.*)/ && (! defined($ncbi_taxid))) {
872
950
            my $taxstring = $1;
873
 
            if ($taxstring =~ /NCBI_TaxID=(.*)/) {
 
951
            # we only keep the first one and ignore all others
 
952
            if ($taxstring =~ /NCBI_TaxID=([\w\d]+)/) {
874
953
                $ncbi_taxid = $1;
875
 
            }
876
 
            else {
 
954
            } else {
877
955
                $self->throw("$taxstring doesn't look like NCBI_TaxID");
878
956
            }
879
957
        }
880
 
        else {
881
 
            last;
882
 
        }
883
958
        
884
959
        $_ = undef; # Empty $_ to trigger read of next line
885
960
    }
889
964
    # Don't make a species object if it is "Unknown" or "None"
890
965
    return if $genus =~ /^(Unknown|None)$/i;
891
966
 
892
 
    if ($class[0] !~ /viruses/i) { # different OS / OC syntax for viruses
893
 
      # Bio::Species array needs array in Species -> Kingdom direction
894
 
      if ($class[$#class] eq $genus) {
 
967
    if ($class[0] eq 'Viruses') {
 
968
        push( @class, $ns_name );
 
969
    }
 
970
    elsif ($class[0] eq 'Viruses') {
 
971
        push( @class, $ns_name );
 
972
    }
 
973
    elsif ($class[$#class] eq $genus) {
895
974
        push( @class, $species );
896
 
      } else {
 
975
    } else {
897
976
        push( @class, $genus, $species );
898
 
      }
899
977
    }
 
978
 
900
979
    @class = reverse @class;
901
980
    
902
 
    my $make = Bio::Species->new();
903
 
    $make->classification( @class );
904
 
    $make->common_name( $common      ) if $common;
905
 
    $make->sub_species( $sub_species ) if $sub_species;
906
 
    $make->organelle  ( $org         ) if $org;
907
 
    $make->ncbi_taxid ( $ncbi_taxid  ) if $ncbi_taxid;
908
 
    return $make;
 
981
    my $taxon = Bio::Species->new();
 
982
    $taxon->classification( \@class, "FORCE" ); # no name validation please
 
983
    $taxon->common_name( $common      ) if $common;
 
984
    $taxon->sub_species( $subspecies ) if $subspecies;
 
985
    $taxon->organelle  ( $org         ) if $org;
 
986
    $taxon->ncbi_taxid ( $ncbi_taxid  ) if $ncbi_taxid;
 
987
    $taxon->variant($variant)           if $variant;
 
988
 
 
989
    # done
 
990
    return $taxon;
909
991
}
910
992
 
911
993
=head2 _filehandle