File Coverage

Bio/SeqIO/swiss.pm
Criterion Covered Total %
statement 471 591 79.7
branch 222 330 67.2
condition 64 100 64.0
subroutine 28 32 87.5
pod 2 2 100.0
total 787 1055 74.6


line stmt bran cond sub pod time code
1             #
2             # BioPerl module for Bio::SeqIO::swiss
3             #
4             # Copyright Elia Stupka
5             #
6             # You may distribute this module under the same terms as perl itself
7              
8             # POD documentation - main docs before the code
9              
10             =head1 NAME
11              
12             Bio::SeqIO::swiss - Swissprot sequence input/output stream
13              
14             =head1 SYNOPSIS
15              
16             It is probably best not to use this object directly, but
17             rather go through the SeqIO handler system:
18              
19             use Bio::SeqIO;
20              
21             $stream = Bio::SeqIO->new(-file => $filename,
22             -format => 'swiss');
23              
24             while ( my $seq = $stream->next_seq() ) {
25             # do something with $seq
26             }
27              
28             =head1 DESCRIPTION
29              
30             This object can transform Bio::Seq objects to and from Swiss-Pprot flat
31             file databases.
32              
33             There is a lot of flexibility here about how to dump things which needs
34             to be documented.
35              
36             =head2 GN (Gene name) line management details
37              
38             A Uniprot/Swiss-Prot entry holds information on one protein
39             sequence. If that sequence is identical across genes and species, they
40             are all merged into one entry. This creates complex needs for several
41             annotation fields in swiss-prot format.
42              
43             The latest syntax for GN line is described in the user manual:
44              
45             http://www.expasy.ch/sprot/userman.html#GN_line
46              
47             Each of the possibly multiple genes in an entry can have Name,
48             Synonyms (only if there is a name), OrderedLocusNames (names from
49             genomic sequences) and ORFNames (temporary or cosmid names). "Name"
50             here really means "symbol". This complexity is now dealt with the
51             following way:
52              
53             A new Bio::AnnotationI class was created in order to store the
54             data in tag-value pairs. This class (Bio::Annotation::TagTree)
55             is stored in the Bio::Annotation::Collection object and is
56             accessed like all other annotations. The tag name is 'gene_name'.
57              
58             There is a single Bio::Annotation::TagTree per sequence record, which
59             corresponds to the original class that stored this data
60             (Bio::Annotation::StructuredValue). Depending on how we progress
61             this may change to represent each group of gene names.
62              
63             For now, to access the gene name tree annotation, one uses the below method:
64              
65             my ($gene) = $seq->annotation->get_Annotations('gene_name');
66              
67             If you are only interested in displaying the values, value() returns a
68             string with similar formatting.
69              
70             There are several ways to get directly at the information you want if you
71             know the element (tag) for the data. For gene names all data is stored with
72             the element-tag pairs:
73              
74             "element1=tag1, tag2, tag3; element2=tag4, tag5;"
75              
76             This normally means the element will be 'Name', 'Synonyms', etc. and the
77             gene names the values. Using findval(), you can do the following:
78              
79             # grab a flattened list of all gene names
80             my @names = $ann->findval('Name');
81              
82             # or iterated through the nodes and grab the name for each group
83             for my $node ($ann->findnode('gene_name')) {
84             my @names = $node->findval('Name');
85             }
86              
87             The current method for parsing gene name data (and reconstructing gene name
88             output) is very generic. This is somewhat preemptive if, for instance, UniProt
89             decides to update and add another element name to the current ones using the
90             same formatting layout. Under those circumstances, one can iterate through the
91             tag tree in a safe way and retrieve all node data like so.
92              
93             # retrieve the gene name nodes (groups like names, synonyms, etc).
94             for my $ann ($seq->annotation->get_Annotations('gene_name')) {
95              
96             # each gene name group
97             for my $node ($ann->findnode('gene_name')) {
98             print "Gene name:\n";
99              
100             # each gene name node (tag => value pair)
101             for my $n ($node->children) {
102             print "\t".$n->element.": ".$n->children."\n";
103             }
104             }
105             }
106              
107             For more uses see Bio::Annotation::TagTree.
108              
109             Since Uniprot/Swiss-Prot format have been around for quite some time, the
110             parser is also able to read in the older GN line syntax where genes
111             are separated by AND and various symbols by OR. The first symbol is
112             taken to be the 'Name' and the remaining ones are stored as 'Synonyms'.
113              
114             Also, for UniProt output we support using other Bio::AnnotationI, but in this
115             case we only use the stringified version of the annotation. This is to allow for
116             backwards compatibility with code that previously used
117             Bio::Annotation::SimpleValue or other Bio::AnnotationI classes.
118              
119             =head2 Optional functions
120              
121             =over 3
122              
123             =item _show_dna()
124              
125             (output only) shows the dna or not
126              
127             =item _post_sort()
128              
129             (output only) provides a sorting func which is applied to the FTHelpers
130             before printing
131              
132             =item _id_generation_func()
133              
134             This is function which is called as
135              
136             print "ID ", $func($seq), "\n";
137              
138             To generate the ID line. If it is not there, it generates a sensible ID
139             line using a number of tools.
140              
141             If you want to output annotations in Swissprot format they need to be
142             stored in a Bio::Annotation::Collection object which is accessible
143             through the Bio::SeqI interface method L.
144              
145             The following are the names of the keys which are polled from a
146             L object.
147              
148             reference - Should contain Bio::Annotation::Reference objects
149             comment - Should contain Bio::Annotation::Comment objects
150             dblink - Should contain Bio::Annotation::DBLink objects
151             gene_name - Should contain Bio::Annotation::SimpleValue object
152              
153             =back
154              
155             =head1 FEEDBACK
156              
157             =head2 Mailing Lists
158              
159             User feedback is an integral part of the evolution of this
160             and other Bioperl modules. Send your comments and suggestions,
161             preferably to one of the Bioperl mailing lists.
162             Your participation is much appreciated.
163              
164             bioperl-l@bioperl.org - General discussion
165             http://bioperl.org/wiki/Mailing_lists - About the mailing lists
166              
167             =head2 Support
168              
169             Please direct usage questions or support issues to the mailing list:
170              
171             I
172              
173             rather than to the module maintainer directly. Many experienced and
174             reponsive experts will be able look at the problem and quickly
175             address it. Please include a thorough description of the problem
176             with code and data examples if at all possible.
177              
178             =head2 Reporting Bugs
179              
180             Report bugs to the Bioperl bug tracking system to help us keep track
181             the bugs and their resolution.
182             Bug reports can be submitted via the web:
183              
184             https://github.com/bioperl/bioperl-live/issues
185              
186             =head1 AUTHOR - Elia Stupka
187              
188             Email elia@tll.org.sg
189              
190             =head1 APPENDIX
191              
192             The rest of the documentation details each of the object methods.
193             Internal methods are usually preceded with a _
194              
195             =cut
196              
197             # Let the code begin...
198              
199             package Bio::SeqIO::swiss;
200 5     5   621 use vars qw(@Unknown_names @Unknown_genus);
  5         6  
  5         314  
201 5     5   20 use strict;
  5         6  
  5         97  
202 5     5   798 use Bio::SeqIO::FTHelper;
  5         8  
  5         112  
203 5     5   22 use Bio::SeqFeature::Generic;
  5         5  
  5         83  
204 5     5   886 use Bio::Species;
  5         7  
  5         102  
205 5     5   2687 use Bio::Tools::SeqStats;
  5         10  
  5         123  
206 5     5   26 use Bio::Seq::SeqFactory;
  5         7  
  5         91  
207 5     5   16 use Bio::Annotation::Collection;
  5         6  
  5         77  
208 5     5   764 use Bio::Annotation::Comment;
  5         6  
  5         92  
209 5     5   840 use Bio::Annotation::Reference;
  5         5  
  5         96  
210 5     5   17 use Bio::Annotation::DBLink;
  5         6  
  5         70  
211 5     5   14 use Bio::Annotation::SimpleValue;
  5         7  
  5         80  
212 5     5   1331 use Bio::Annotation::TagTree;
  5         9  
  5         123  
213              
214 5     5   19 use base qw(Bio::SeqIO);
  5         106  
  5         24467  
215              
216             our $LINE_LENGTH = 76;
217              
218             # this is for doing species name parsing
219             @Unknown_names=('other', 'unidentified',
220             'unknown organism', 'not specified',
221             'not shown', 'Unspecified', 'Unknown',
222             'None', 'unclassified', 'unidentified organism',
223             'not supplied'
224             );
225             # dictionary of synonyms for taxid 32644
226             # all above can be part of valid species name
227             @Unknown_genus = qw(unknown unclassified uncultured unidentified);
228              
229             # if there are any other gene name tags, they are added to the end
230             our @GENE_NAME_ORDER = qw(Name Synonyms OrderedLocusNames ORFNames);
231              
232             sub _initialize {
233 23     23   50 my($self,@args) = @_;
234 23         89 $self->SUPER::_initialize(@args);
235             # hash for functions for decoding keys.
236 23         69 $self->{'_func_ftunit_hash'} = {};
237             # sets this to one by default. People can change it
238 23         83 $self->_show_dna(1);
239 23 50       71 if ( ! defined $self->sequence_factory ) {
240 23         85 $self->sequence_factory(Bio::Seq::SeqFactory->new
241             (-verbose => $self->verbose(),
242             -type => 'Bio::Seq::RichSeq'));
243             }
244             }
245              
246             =head2 next_seq
247              
248             Title : next_seq
249             Usage : $seq = $stream->next_seq()
250             Function: returns the next sequence in the stream
251             Returns : Bio::Seq object
252             Args :
253              
254             =cut
255              
256             sub next_seq {
257 36     36 1 2024 my ($self,@args) = @_;
258 36         39 my ($pseq,$c,$line,$name,$desc,$acc,$seqc,$mol,$div, $sptr,$seq_div,
259             $date,$comment,@date_arr);
260 36         52 my $genename = "";
261 36         186 my ($annotation, %params, @features) = ( Bio::Annotation::Collection->new());
262              
263 36         54 local $_;
264              
265 36   100     124 1 while defined($_ = $self->_readline) && /^\s+$/;
266 36 100 66     222 return unless defined $_ && /^ID\s/;
267              
268             # fixed to allow _DIVISION to be optional for bug #946
269             # see bug report for more information
270             #
271             # 9/6/06 Note: Swiss/TrEMBL sequences have no division acc. to UniProt
272             # release notes; this is fixed to simplify the regex parsing
273             # STANDARD (SwissProt) and PRELIMINARY (TrEMBL) added to namespace()
274 31 50       208 unless( m{^
275             ID \s+ #
276             (\S+) \s+ # $1 entryname
277             ([^\s;]+); \s+ # $2 DataClass
278             (?:PRT;)? \s+ # Molecule Type (optional)
279             [0-9]+[ ]AA \. # Sequencelength (capture?)
280             $
281             }ox ) {
282             # I couldn't find any new current UniProt sequences
283             # that matched this format:
284             # || m/^ID\s+(\S+)\s+(_([^\s_]+))? /ox ) {
285 0         0 $self->throw("swissprot stream with no ID. Not swissprot in my book");
286             }
287 31         98 ($name, $seq_div) = ($1, $2);
288 31 100 100     202 $params{'-namespace'} =
    100 100        
289             ($seq_div eq 'Reviewed' || $seq_div eq 'STANDARD') ? 'Swiss-Prot' :
290             ($seq_div eq 'Unreviewed' || $seq_div eq 'PRELIMINARY') ? 'TrEMBL' :
291             $seq_div;
292             # we shouldn't be setting the division, but for now...
293 31         100 my ($junk, $division) = split q(_), $name;
294 31         54 $params{'-division'} = $division;
295 31         66 $params{'-alphabet'} = 'protein';
296             # this is important to have the id for display in e.g. FTHelper, otherwise
297             # you won't know which entry caused an error
298 31         54 $params{'-display_id'} = $name;
299              
300             BEFORE_FEATURE_TABLE :
301 31         86 while ( defined($_ = $self->_readline) ) {
302             # Exit at start of Feature table and at the sequence at the
303             # latest HL 05/11/2000
304 981 100       2279 last if( /^(FT|SQ)/ );
305              
306             # Description line(s)
307 950 100       6317 if (/^DE\s+(\S.*\S)/) {
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    50          
308 54 100       229 $desc .= $desc ? " $1" : $1;
309             }
310             #Gene name
311             elsif (/^GN\s+(.*)/) {
312 36 100       67 $genename .= " " if $genename;
313 36         89 $genename .= $1;
314             }
315             #accession number(s)
316             elsif ( /^AC\s+(.+)/) {
317 31         130 my @accs = split(/[; ]+/, $1); # allow space in addition
318             $params{'-accession_number'} = shift @accs
319 31 50       106 unless defined $params{'-accession_number'};
320 31         39 push @{$params{'-secondary_accessions'}}, @accs;
  31         110  
321             }
322             #date and sequence version
323             elsif ( /^DT\s+(.*)/ ) {
324 90         135 my $line = $1;
325 90         186 my ($date, $version) = split(' ', $line, 2);
326 90         106 $date =~ tr/,//d; # remove comma if new version
327 90 100 100     616 if ($version =~ /\(Rel\. (\d+), Last sequence update\)/ || # old
    100 100        
328             /sequence version (\d+)/) { #new
329 30         176 my $update = Bio::Annotation::SimpleValue->new
330             (-tagname => 'seq_update',
331             -value => $1
332             );
333 30         112 $annotation->add_Annotation($update);
334             } elsif ($version =~ /\(Rel\. (\d+), Last annotation update\)/ || #old
335             /entry version (\d+)/) { #new
336 30         72 $params{'-version'} = $1;
337             }
338 90         85 push @{$params{'-dates'}}, $date;
  90         263  
339             }
340             # Evidence level
341             elsif ( /^PE\s+(.*)/ ) {
342 1         3 my $line = $1;
343 1         4 $line =~ s/;\s*//; # trim trailing semicolon and any spaces at the end of the line
344 1         8 my $evidence = Bio::Annotation::SimpleValue->new
345             (-tagname => 'evidence',
346             -value => $line
347             );
348 1         4 $annotation->add_Annotation($evidence);
349             }
350             # Organism name and phylogenetic information
351             elsif (/^O[SCG]/) {
352 30         96 my $species = $self->_read_swissprot_Species($_);
353 30         146 $params{'-species'}= $species;
354             # now we are one line ahead -- so continue without reading the next
355             # line HL 05/11/2000
356             }
357             # References
358             elsif (/^R/) {
359 30         85 my $refs = $self->_read_swissprot_References($_);
360 30         54 foreach my $r (@$refs) {
361 227         368 $annotation->add_Annotation('reference',$r);
362             }
363             }
364             # Comments
365             elsif (/^CC\s{3}(.*)/) {
366 30         88 $comment .= $1;
367 30         46 $comment .= "\n";
368 30   66     77 while (defined ($_ = $self->_readline) && /^CC\s{3}(.*)/ ) {
369 474         1080 $comment .= $1 . "\n";
370             }
371              
372             # note: don't try to process comments here -- they may contain
373             # structure. LP 07/30/2000
374              
375 30         241 my $commobj = Bio::Annotation::Comment->new(-tagname => 'comment',
376             -text => $comment);
377 30         88 $annotation->add_Annotation('comment',$commobj);
378 30         37 $comment = "";
379 30         117 $self->_pushback($_);
380             }
381             #DBLinks
382             # old regexp
383             # /^DR\s+(\S+)\;\s+(\S+)\;\s+(\S+)[\;\.](.*)$/) {
384             # new regexp from Andreas Kahari bug #1584
385             elsif (/^DR\s+(\S+)\;\s+(\S+)\;\s+([^;]+)[\;\.](.*)$/) {
386 599         1486 my ($database,$primaryid,$optional,$comment) = ($1,$2,$3,$4);
387              
388             # drop leading and training spaces and trailing .
389 599         1042 $comment =~ s/\.\s*$//;
390 599         690 $comment =~ s/^\s+//;
391              
392 599         2047 my $dblinkobj = Bio::Annotation::DBLink->new
393             (-database => $database,
394             -primary_id => $primaryid,
395             -optional_id => $optional,
396             -comment => $comment,
397             -tagname => 'dblink',
398             );
399              
400 599         1357 $annotation->add_Annotation('dblink',$dblinkobj);
401             }
402             #keywords
403             elsif ( /^KW\s+(.*)$/ ) {
404 49         332 my @kw = split(/\s*\;\s*/,$1);
405 49 50       175 defined $kw[-1] && $kw[-1] =~ s/\.$//;
406 49         57 push @{$params{'-keywords'}}, @kw;
  49         205  
407             }
408             }
409             # process and parse the gene name line if there was one (note: we
410             # can't do this above b/c GN may be multi-line and we can't
411             # unequivocally determine whether we've seen the last GN line in
412             # the new format)
413 31 100       82 if ($genename) {
414 30         31 my @stags;
415 30 100       84 if ($genename =~ /\w=\w/) {
416             # new format (e.g., Name=RCHY1; Synonyms=ZNF363, CHIMP)
417 12         47 for my $n (split(m{\s+and\s+},$genename)) {
418 14         16 my @genenames;
419 14         61 for my $section (split(m{\s*;\s*},$n)) {
420 19         48 my ($tag, $rest) = split("=",$section);
421 19   50     41 $rest ||= '';
422 19         33 for my $val (split(m{\s*,\s*},$rest)) {
423 23         54 push @genenames, [$tag => $val];
424             }
425             }
426 14         45 push @stags, ['gene_name' => \@genenames];
427             }
428             } else {
429             # old format
430 18         66 for my $section (split(/ AND /, $genename)) {
431 22         19 my @genenames;
432 22         74 $section =~ s/[\(\)\.]//g;
433 22         73 my @names = split(m{\s+OR\s+}, $section);
434 22         57 push @genenames, ['Name' => shift @names];
435 22         36 push @genenames, map {['Synonyms' => $_]} @names;
  25         43  
436 22         75 push @stags, ['gene_name' => \@genenames]
437             }
438             } #use Data::Dumper; print Dumper $gn, $genename;# exit;
439 30         258 my $gn = Bio::Annotation::TagTree->new(-tagname => 'gene_name',
440             -value => ['gene_names' => \@stags]);
441 30         104 $annotation->add_Annotation('gene_name', $gn);
442             }
443              
444             FEATURE_TABLE :
445             # if there is no feature table, or if we've got beyond, exit loop or don't
446             # even enter HL 05/11/2000
447 31   66     209 while (defined $_ && /^FT/ ) {
448 443         835 my $ftunit = $self->_read_FTHelper_swissprot($_);
449              
450             # process ftunit
451             # when parsing of the line fails we get undef returned
452 443 50       567 if ($ftunit) {
453             push(@features,
454             $ftunit->_generic_seqfeature($self->location_factory(),
455 443         939 $params{'-seqid'}, "SwissProt"));
456             } else {
457             $self->warn("failed to parse feature table line for seq " .
458 0         0 $params{'-display_id'}. "\n$_");
459             }
460 443         1198 $_ = $self->_readline;
461             }
462 31   33     213 while ( defined($_) && ! /^SQ/ ) {
463 0         0 $_ = $self->_readline;
464             }
465 31         64 $seqc = "";
466 31         69 while ( defined ($_ = $self->_readline) ) {
467 242 100       439 last if m{^//};
468 211         1105 s/[^A-Za-z]//g;
469 211         428 $seqc .= uc($_);
470             }
471              
472 31         112 my $seq= $self->sequence_factory->create
473             (-verbose => $self->verbose,
474             %params,
475             -seq => $seqc,
476             -desc => $desc,
477             -features => \@features,
478             -annotation => $annotation,
479             );
480              
481             # The annotation doesn't get added by the contructor
482 31         103 $seq->annotation($annotation);
483              
484 31         218 return $seq;
485             }
486              
487             =head2 write_seq
488              
489             Title : write_seq
490             Usage : $stream->write_seq($seq)
491             Function: writes the $seq object (must be seq) to the stream
492             Returns : 1 for success and 0 for error
493             Args : array of 1 to n Bio::SeqI objects
494              
495              
496             =cut
497              
498             sub write_seq {
499 4     4 1 49 my ($self,@seqs) = @_;
500 4         10 foreach my $seq ( @seqs ) {
501 4 50       12 $self->throw("Attempting to write with no seq!") unless defined $seq;
502              
503 4 50 33     30 if ( ! ref $seq || ! $seq->isa('Bio::SeqI') ) {
504 0         0 $self->warn(" $seq is not a SeqI compliant module. Attempting to dump, but may fail!");
505             }
506              
507 4         4 my $i;
508 4         20 my $str = $seq->seq;
509              
510 4         6 my $div;
511 4   66     51 my $ns = ($seq->can('namespace')) && $seq->namespace();
512 4         17 my $len = $seq->length();
513              
514 4 100 66     29 if ( !$seq->can('division') || ! defined ($div = $seq->division()) ) {
515 1         3 $div = 'UNK';
516             }
517              
518             # namespace dictates database, takes precedent over division. Sorry!
519 4 100 66     26 if (defined($ns) && $ns ne '') {
520 3 0       11 $div = ($ns eq 'Swiss-Prot') ? 'Reviewed' :
    50          
521             ($ns eq 'TrEMBL') ? 'Unreviewed' :
522             $ns;
523             } else {
524 1         3 $ns = 'Swiss-Prot';
525             # division not reset; acts as fallback
526             }
527              
528 4 50       16 $self->warn("No whitespace allowed in SWISS-PROT display id [". $seq->display_id. "]")
529             if $seq->display_id =~ /\s/;
530              
531 4         6 my $temp_line;
532 4 50       13 if ( $self->_id_generation_func ) {
533 0         0 $temp_line = &{$self->_id_generation_func}($seq);
  0         0  
534             } else {
535             #$temp_line = sprintf ("%10s STANDARD; %3s; %d AA.",
536             # $seq->primary_id()."_".$div,$mol,$len);
537             # Reconstructing the ID relies heavily upon the input source having
538             # been in a format that is parsed as this routine expects it -- that is,
539             # by this module itself. This is bad, I think, and immediately breaks
540             # if e.g. the Bio::DB::GenPept module is used as input.
541             # Hence, switch to display_id(); _every_ sequence is supposed to have
542             # this. HL 2000/09/03
543             # Changed to reflect ID line changes in UniProt
544             # Oct 2006 - removal of molecule type - see bug 2134
545 4         12 $temp_line = sprintf ("%-24s%-12s%9d AA.",
546             $seq->display_id(), $div.';', $len);
547             }
548              
549 4         27 $self->_print( "ID $temp_line\n");
550              
551             # if there, write the accession line
552 4         22 local($^W) = 0; # supressing warnings about uninitialized fields
553              
554 4 50       14 if ( $self->_ac_generation_func ) {
    50          
555 0         0 $temp_line = &{$self->_ac_generation_func}($seq);
  0         0  
556 0         0 $self->_print( "AC $temp_line\n");
557             }
558             elsif ($seq->can('accession_number') ) {
559 4         14 my $ac_line = $seq->accession_number;
560 4 100       24 if ($seq->can('get_secondary_accessions') ) {
561 3         10 foreach my $sacc ($seq->get_secondary_accessions) {
562 0         0 $ac_line .= "; ". $sacc;;
563             }
564 3         9 $ac_line .= ";";
565             }
566              
567 4         16 $self->_write_line_swissprot_regex("AC ","AC ",$ac_line,
568             "\\s\+\|\$",$LINE_LENGTH);
569             }
570             # otherwise - cannot print
571              
572              
573             # Date lines and sequence versions (changed 6/15/2006)
574             # This is rebuilt from scratch using the current SwissProt/UniProt format
575 4 100       25 if ( $seq->can('get_dates') ) {
576 3         12 my @dates = $seq->get_dates();
577 3         7 my $ct = 1;
578 3         13 my $seq_version = $seq->version;
579 3         8 my ($update_version) = $seq->annotation->get_Annotations("seq_update");
580 3         6 foreach my $dt (@dates) {
581 9 100       27 $self->_write_line_swissprot_regex("DT ","DT ",
582             $dt.', integrated into UniProtKB/'.$ns.'.',
583             "\\s\+\|\$",$LINE_LENGTH) if $ct == 1;
584 9 100       30 $self->_write_line_swissprot_regex("DT ","DT ",
585             $dt.", sequence version ".$update_version->display_text.'.',
586             "\\s\+\|\$",$LINE_LENGTH) if $ct == 2;
587 9 100       25 $self->_write_line_swissprot_regex("DT ","DT ",
588             $dt.", entry version $seq_version.",
589             "\\s\+\|\$",$LINE_LENGTH) if $ct == 3;
590 9         9 $ct++;
591             }
592             }
593              
594             #Definition lines
595 4         22 $self->_write_line_swissprot_regex("DE ","DE ",$seq->desc(),"\\s\+\|\$",$LINE_LENGTH);
596              
597             #Gene name; print out new format
598 4         13 foreach my $gene ( my @genes = $seq->annotation->get_Annotations('gene_name') ) {
599             # gene is a Bio::Annotation::TagTree;
600 3 100       19 if ($gene->isa('Bio::Annotation::TagTree')) {
601 2         4 my @genelines;
602 2         8 for my $node ($gene->findnode('gene_name')) {
603             # check for Name and Synonym first, then the rest get tagged on
604 2         205 my $geneline = "GN ";
605 2         10 my %genedata = $node->hash;
606 2         53 for my $tag (@GENE_NAME_ORDER) {
607 8 100       16 if (exists $genedata{$tag}) {
608             $geneline .= (ref $genedata{$tag} eq 'ARRAY') ?
609 2 50       9 "$tag=".join(', ',@{$genedata{$tag}})."; " :
  0         0  
610             "$tag=$genedata{$tag}; ";
611 2         5 delete $genedata{$tag};
612             }
613             }
614             # add rest
615 2         8 for my $tag (sort keys %genedata) {
616             $geneline .= (ref $genedata{$tag} eq 'ARRAY') ?
617 0 0       0 "$tag=".join(', ',@{$genedata{$tag}})."; " :
  0         0  
618             "$tag=$genedata{$tag}; ";
619 0         0 delete $genedata{$tag};
620             }
621 2         7 push @genelines, "$geneline\n";
622             }
623 2         10 $self->_print(join("GN and\n",@genelines));
624             } else { # fall back to getting stringified output
625 1         5 $self->_write_line_swissprot_regex("GN ","GN ",
626             $gene->display_text,
627             "\\s\+\|\$",
628             $LINE_LENGTH);
629             }
630             }
631              
632             # Organism lines
633 4 100 66     31 if ($seq->can('species') && (my $spec = $seq->species)) {
634 3         14 my @class = $spec->classification();
635 3         7 shift(@class);
636 3         12 my $species = $spec->species;
637 3         7 my $genus = $spec->genus;
638 3         10 my $OS = $spec->scientific_name;
639 3 50       13 if ($class[-1] =~ /viruses/i) {
640 0         0 $OS = $species;
641 0 0       0 $OS .= " ". $spec->sub_species if $spec->sub_species;
642             }
643 3         11 foreach (($spec->variant, $spec->common_name)) {
644 3 50       9 $OS .= " ($_)" if $_;
645             }
646 3         15 $self->_print( "OS $OS.\n");
647 3         22 my $OC = join('; ', reverse(@class)) .'.';
648 3         9 $self->_write_line_swissprot_regex("OC ","OC ",$OC,"\; \|\$",$LINE_LENGTH);
649 3 50       11 if ($spec->organelle) {
650 0         0 $self->_write_line_swissprot_regex("OG ","OG ",$spec->organelle,"\; \|\$",$LINE_LENGTH);
651             }
652 3 50       10 if ($spec->ncbi_taxid) {
653 3         10 $self->_print("OX NCBI_TaxID=".$spec->ncbi_taxid.";\n");
654             }
655             }
656              
657             # Reference lines
658 4         8 my $t = 1;
659 4         12 foreach my $ref ( $seq->annotation->get_Annotations('reference') ) {
660 3         21 $self->_print( "RN [$t]\n");
661             # changed by lorenz 08/03/00
662             # j.gilbert and h.lapp agreed that the rp line in swissprot seems
663             # more like a comment than a parseable value, so print it as is
664 3 50       10 if ($ref->rp) {
665 3         11 $self->_write_line_swissprot_regex("RP ","RP ",$ref->rp,
666             "\\s\+\|\$",$LINE_LENGTH);
667             }
668 3 100       11 if ($ref->comment) {
669 2         6 $self->_write_line_swissprot_regex("RC ","RC ",$ref->comment,
670             "\\s\+\|\$",$LINE_LENGTH);
671             }
672 3 50 33     12 if ($ref->medline or $ref->pubmed or $ref->doi) {
      33        
673             # new RX format in swissprot LP 09/17/00
674             # RX line can now have a DOI, Heikki 13 Feb 2008
675              
676 0         0 my $line;
677 0 0       0 $line .= "MEDLINE=". $ref->medline. '; ' if $ref->medline;
678 0 0       0 $line .= "PubMed=". $ref->pubmed. '; ' if $ref->pubmed;
679 0 0       0 $line .= "DOI=". $ref->doi. '; ' if $ref->doi;
680 0         0 chop $line;
681              
682 0         0 $self->_write_line_swissprot_regex("RX ","RX ",
683             $line,
684             "\\s\+\|\$",$LINE_LENGTH);
685              
686             }
687 3 50       9 my $author = $ref->authors .';' if($ref->authors);
688 3 50       11 my $title = $ref->title .';' if( $ref->title);
689 3 50       13 my $rg = $ref->rg . ';' if $ref->rg;
690 3         28 $author =~ s/([\w\.]) (\w)/$1#$2/g; # add word wrap protection char '#'
691              
692 3 50       8 $self->_write_line_swissprot_regex("RG ","RG ",$rg,"\\s\+\|\$",$LINE_LENGTH) if $rg;
693 3 50       14 $self->_write_line_swissprot_regex("RA ","RA ",$author,"\\s\+\|\$",$LINE_LENGTH) if $author;
694 3 50       9 $self->_write_line_swissprot_regex("RT ","RT ",$title,'[\s\-]+|$',$LINE_LENGTH) if $title;
695 3         10 $self->_write_line_swissprot_regex("RL ","RL ",$ref->location,"\\s\+\|\$",$LINE_LENGTH);
696 3         7 $t++;
697             }
698              
699             # Comment lines
700              
701 4         13 foreach my $comment ( $seq->annotation->get_Annotations('comment') ) {
702 3         14 foreach my $cline (split ("\n", $comment->text)) {
703 48         62 while (length $cline > 74) {
704 0         0 $self->_print("CC ",(substr $cline,0,74),"\n");
705 0         0 $cline = substr $cline,74;
706             }
707 48         58 $self->_print("CC ",$cline,"\n");
708             }
709             }
710              
711             # Database xref lines
712              
713 4         12 foreach my $dblink ( $seq->annotation->get_Annotations('dblink') ) {
714 27         37 my ($primary_id) = $dblink->primary_id;
715            
716 27 100 66     30 if (defined($dblink->comment) && ($dblink->comment) ) {
    50          
717 12         15 $self->_print("DR ",$dblink->database,"; ",$primary_id,"; ",
718             $dblink->optional_id,"; ",$dblink->comment,".\n");
719             } elsif ($dblink->optional_id) {
720 15         21 $self->_print("DR ",$dblink->database,"; ",
721             $primary_id,"; ",
722             $dblink->optional_id,".\n");
723             } else {
724 0         0 $self->_print("DR ",$dblink->database,
725             "; ",$primary_id,"; ","-.\n");
726             }
727             }
728              
729             # Evidence lines
730              
731 4         13 foreach my $evidence ( $seq->annotation->get_Annotations('evidence') ) {
732 0         0 $self->_print("PE ",$evidence->value,";\n");
733             }
734              
735             # if there, write the kw line
736             {
737 4         4 my $kw;
  4         6  
738 4 50       14 if ( my $func = $self->_kw_generation_func ) {
    100          
739 0         0 $kw = &{$func}($seq);
  0         0  
740             } elsif ( $seq->can('keywords') ) {
741 3         10 $kw = $seq->keywords;
742 3 50       13 if ( ref($kw) =~ /ARRAY/i ) {
743 0         0 $kw = join("; ", @$kw);
744             }
745 3 50 33     30 $kw .= '.' if $kw and $kw !~ /\.$/ ;
746             }
747 4         34 $kw =~ s/([\w\.]) (\w)/$1#$2/g; # add word wrap protection char '#'
748 4 100       19 $self->_write_line_swissprot_regex("KW ","KW ",
749             $kw, "\\s\+\|\$",$LINE_LENGTH)
750             if $kw;
751             }
752              
753             #Check if there is seqfeatures before printing the FT line
754 4 50       29 my @feats = $seq->can('top_SeqFeatures') ? $seq->top_SeqFeatures : ();
755 4 100       12 if ($feats[0]) {
756 3 50       9 if ( defined $self->_post_sort ) {
757              
758             # we need to read things into an array. Process. Sort them. Print 'em
759              
760 0         0 my $post_sort_func = $self->_post_sort();
761 0         0 my @fth;
762              
763 0         0 foreach my $sf ( @feats ) {
764 0         0 push(@fth,Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq));
765             }
766 0         0 @fth = sort { &$post_sort_func($a,$b) } @fth;
  0         0  
767              
768 0         0 foreach my $fth ( @fth ) {
769 0         0 $self->_print_swissprot_FTHelper($fth);
770             }
771             } else {
772             # not post sorted. And so we can print as we get them.
773             # lower memory load...
774              
775 3         9 foreach my $sf ( @feats ) {
776 9         23 my @fth = Bio::SeqIO::FTHelper::from_SeqFeature($sf,$seq);
777 9         10 foreach my $fth ( @fth ) {
778 9 50       22 if ( ! $fth->isa('Bio::SeqIO::FTHelper') ) {
779 0         0 $sf->throw("Cannot process FTHelper... $fth");
780             }
781 9         23 $self->_print_swissprot_FTHelper($fth);
782             }
783             }
784             }
785              
786 3 50       10 if ( $self->_show_dna() == 0 ) {
787 0         0 return;
788             }
789             }
790             # finished printing features.
791              
792             # molecular weight
793 4         7 my $mw = ${Bio::Tools::SeqStats->get_mol_wt($seq->primary_seq)}[0];
  4         17  
794             # checksum
795             # was crc32 checksum, changed it to crc64
796 4         18 my $crc64 = $self->_crc64(\$str);
797 4         69 $self->_print( sprintf("SQ SEQUENCE %4d AA; %d MW; %16s CRC64;\n",
798             $len,$mw,$crc64));
799 4         12 $self->_print( " ");
800 4         5 my $linepos;
801 4         16 for ($i = 0; $i < length($str); $i += 10) {
802 164         259 $self->_print( " ", substr($str,$i,10));
803 164         152 $linepos += 11;
804 164 100 66     424 if ( ($i+10)%60 == 0 && (($i+10) < length($str))) {
805 24         33 $self->_print( "\n ");
806             }
807             }
808 4         11 $self->_print( "\n//\n");
809              
810 4 50 33     14 $self->flush if $self->_flush_on_write && defined $self->_fh;
811 4         47 return 1;
812             }
813             }
814              
815             # Thanks to James Gilbert for the following two. LP 08/01/2000
816              
817             =head2 _generateCRCTable
818              
819             Title : _generateCRCTable
820             Usage :
821             Function:
822             Example :
823             Returns :
824             Args :
825              
826              
827             =cut
828              
829             sub _generateCRCTable {
830             # 10001000001010010010001110000100
831             # 32
832 0     0   0 my $poly = 0xEDB88320;
833 0         0 my ($self) = shift;
834              
835 0         0 $self->{'_crcTable'} = [];
836 0         0 foreach my $i (0..255) {
837 0         0 my $crc = $i;
838 0         0 for (my $j=8; $j > 0; $j--) {
839 0 0       0 if ($crc & 1) {
840 0         0 $crc = ($crc >> 1) ^ $poly;
841             } else {
842 0         0 $crc >>= 1;
843             }
844             }
845 0         0 ${$self->{'_crcTable'}}[$i] = $crc;
  0         0  
846             }
847             }
848              
849              
850             =head2 _crc32
851              
852             Title : _crc32
853             Usage :
854             Function:
855             Example :
856             Returns :
857             Args :
858              
859              
860             =cut
861              
862             sub _crc32 {
863 0     0   0 my( $self, $str ) = @_;
864              
865 0 0       0 $self->throw("Argument to crc32() must be ref to scalar")
866             unless ref($str) eq 'SCALAR';
867              
868 0 0       0 $self->_generateCRCTable() unless exists $self->{'_crcTable'};
869              
870 0         0 my $len = length($$str);
871              
872 0         0 my $crc = 0xFFFFFFFF;
873 0         0 for (my $i = 0; $i < $len; $i++) {
874             # Get upper case value of each letter
875 0         0 my $int = ord uc substr $$str, $i, 1;
876             $crc = (($crc >> 8) & 0x00FFFFFF) ^
877 0         0 ${$self->{'_crcTable'}}[ ($crc ^ $int) & 0xFF ];
  0         0  
878             }
879 0         0 return $crc;
880             }
881              
882             =head2 _crc64
883              
884             Title : _crc64
885             Usage :
886             Function:
887             Example :
888             Returns :
889             Args :
890              
891              
892             =cut
893              
894             sub _crc64{
895 4     4   8 my ($self, $sequence) = @_;
896 4         6 my $POLY64REVh = 0xd8000000;
897 4         9 my @CRCTableh = 256;
898 4         7 my @CRCTablel = 256;
899 4         3 my $initialized;
900              
901 4         5 my $seq = $$sequence;
902              
903 4         6 my $crcl = 0;
904 4         4 my $crch = 0;
905 4 50       12 if (!$initialized) {
906 4         7 $initialized = 1;
907 4         14 for (my $i=0; $i<256; $i++) {
908 1024         593 my $partl = $i;
909 1024         607 my $parth = 0;
910 1024         1145 for (my $j=0; $j<8; $j++) {
911 8192         5115 my $rflag = $partl & 1;
912 8192         4521 $partl >>= 1;
913 8192 50       8360 $partl |= (1 << 31) if $parth & 1;
914 8192         4543 $parth >>= 1;
915 8192 100       13120 $parth ^= $POLY64REVh if $rflag;
916             }
917 1024         778 $CRCTableh[$i] = $parth;
918 1024         1278 $CRCTablel[$i] = $partl;
919             }
920             }
921              
922 4         307 foreach (split '', $seq) {
923 1636         1017 my $shr = ($crch & 0xFF) << 24;
924 1636         950 my $temp1h = $crch >> 8;
925 1636         962 my $temp1l = ($crcl >> 8) | $shr;
926 1636         1241 my $tableindex = ($crcl ^ (unpack "C", $_)) & 0xFF;
927 1636         1069 $crch = $temp1h ^ $CRCTableh[$tableindex];
928 1636         1241 $crcl = $temp1l ^ $CRCTablel[$tableindex];
929             }
930 4         94 my $crc64 = sprintf("%08X%08X", $crch, $crcl);
931 4         51 return $crc64;
932             }
933              
934             =head2 _print_swissprot_FTHelper
935              
936             Title : _print_swissprot_FTHelper
937             Usage :
938             Function:
939             Example :
940             Returns :
941             Args :
942              
943              
944             =cut
945              
946             sub _print_swissprot_FTHelper {
947 9     9   9 my ($self,$fth,$always_quote) = @_;
948 9   50     27 $always_quote ||= 0;
949 9         11 my ($start,$end) = ('?', '?');
950              
951 9 50 33     39 if ( ! ref $fth || ! $fth->isa('Bio::SeqIO::FTHelper') ) {
952 0         0 $fth->warn("$fth is not a FTHelper class. ".
953             "Attempting to print, but there could be tears!");
954             }
955 9         9 my $desc = "";
956              
957 9         14 for my $tag ( qw(description gene note product) ) {
958 9 50       13 if ( exists $fth->field->{$tag} ) {
959 9         7 $desc = @{$fth->field->{$tag}}[0].".";
  9         14  
960 9         12 last;
961             }
962             }
963 9         30 $desc =~ s/\.$//;
964              
965 9         9 my $ftid = "";
966 9 50       13 if ( exists $fth->field->{'FTId'} ) {
967 0         0 $ftid = @{$fth->field->{'FTId'}}[0]. '.';
  0         0  
968             }
969              
970 9         17 my $key =substr($fth->key,0,8);
971 9         17 my $loc = $fth->loc;
972 9 100       47 if ( $loc =~ /(\?|\d+|\>\d+|<\d+)?\.\.(\?|\d+|<\d+|>\d+)?/ ) {
    50          
973 6 50       20 $start = $1 if defined $1;
974 6 50       17 $end = $2 if defined $2;
975              
976             # to_FTString only returns one value when start == end, #JB955
977             # so if no match is found, assume it is both start and end #JB955
978             } elsif ( $loc =~ /join\((\d+)((?:,\d+)+)?\)/) {
979 0         0 my @y = ($1);
980 0 0       0 if ( defined( my $m = $2) ) {
981 0         0 $m =~ s/^\,//;
982 0         0 push @y, split(/,/,$m);
983             }
984 0         0 for my $x ( @y ) {
985 0         0 $self->_write_line_swissprot_regex(
986             sprintf("FT %-8s %6s %6s ",
987             $key,
988             $x ,$x),
989             "FT ",
990             $desc.'.','\s+|$',$LINE_LENGTH);
991             }
992 0         0 return;
993             } else {
994 3         7 $start = $end = $fth->loc;
995             }
996 9 50       17 if ($desc) {
997 9         56 $self->_write_line_swissprot_regex(sprintf("FT %-8s %6s %6s ",
998             $key,
999             $start ,$end),
1000             "FT ",
1001             $desc. '.', '\s+|$', $LINE_LENGTH);
1002             } else { #HELIX and STRAND do not have descriptions
1003 0         0 $self->_write_line_swissprot_regex(sprintf("FT %-8s %6s %6s",
1004             $key,
1005             $start ,$end),
1006             "FT ",
1007             ' ', '\s+|$', $LINE_LENGTH);
1008             }
1009              
1010              
1011 9 50       43 if ($ftid) {
1012 0         0 $self->_write_line_swissprot_regex("FT ",
1013             "FT ",
1014             "/FTId=$ftid",'.|$',$LINE_LENGTH);
1015              
1016             }
1017              
1018             }
1019             #'
1020              
1021             =head2 _read_swissprot_References
1022              
1023             Title : _read_swissprot_References
1024             Usage :
1025             Function: Reads references from swissprot format. Internal function really
1026             Example :
1027             Returns :
1028             Args :
1029              
1030              
1031             =cut
1032              
1033             sub _read_swissprot_References{
1034 30     30   41 my ($self,$line) = @_;
1035 30         38 my ($b1, $b2, $rp, $rg, $title, $loc, $au, $med, $com, $pubmed, $doi);
1036 0         0 my @refs;
1037 30         46 local $_ = $line;
1038 30         82 while ( defined $_ ) {
1039 1813 100 100     10514 if ( /^[^R]/ || /^RN/ ) {
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    50          
1040 257 100       360 if ( $rp ) {
1041 227 100       354 $rg =~ s/;\s*$//g if defined($rg);
1042 227 100       246 if (defined($au)) {
1043 223         554 $au =~ s/;\s*$//;
1044             } else {
1045 4         7 $au = $rg;
1046             }
1047 227 100       610 $title =~ s/;\s*$//g if defined($title);
1048 227         1239 push @refs, Bio::Annotation::Reference->new
1049             (-title => $title,
1050             -start => $b1,
1051             -end => $b2,
1052             -authors => $au,
1053             -location=> $loc,
1054             -medline => $med,
1055             -pubmed => $pubmed,
1056             -doi => $doi,
1057             -comment => $com,
1058             -rp => $rp,
1059             -rg => $rg,
1060             -tagname => 'reference',
1061             );
1062             # reset state for the next reference
1063 227         372 $rp = '';
1064             }
1065 257 100       559 if (index($_,'R') != 0) {
1066 30         83 $self->_pushback($_); # want this line to go back on the list
1067 30         66 last; # may be the safest exit point HL 05/11/2000
1068             }
1069             # don't forget to reset the state for the next reference
1070 227         263 $b1 = $b2 = $rg = $med = $com = $pubmed = $doi = undef;
1071 227         213 $title = $loc = $au = undef;
1072             } elsif ( /^RP\s{3}(.+? OF (\d+)-(\d+).*)/) {
1073 35         86 $rp .= $1;
1074 35         46 $b1 = $2;
1075 35         37 $b2 = $3;
1076             } elsif ( /^RP\s{3}(.*)/) {
1077 198 100       251 if ($rp) {
1078 6         18 $rp .= " ".$1;
1079             } else {
1080 192         359 $rp = $1;
1081             }
1082             } elsif (/^RX\s{3}(.*)/) { # each reference can have only one RX line
1083 177         252 my $line = $1;
1084 177 100       550 $med = $1 if $line =~ /MEDLINE=(\d+);/;
1085 177 100       506 $pubmed = $1 if $line =~ /PubMed=(\d+);/;
1086 177 100       316 $doi = $1 if $line =~ /DOI=(.+);/;
1087             } elsif ( /^RA\s{3}(.*)/ ) {
1088 382 100       887 $au .= $au ? " $1" : $1;
1089             } elsif ( /^RG\s{3}(.*)/ ) {
1090 16 50       51 $rg .= $rg ? " $1" : $1;
1091             } elsif ( /^RT\s{3}(.*)/ ) {
1092 381 100       433 if ($title) {
1093 177         245 my $tline = $1;
1094 177 100       580 $title .= ($title =~ /[\w;,:\?!]$/) ? " $tline" : $tline;
1095             } else {
1096 204         289 $title = $1;
1097             }
1098             } elsif (/^RL\s{3}(.*)/ ) {
1099 227 50       500 $loc .= $loc ? " $1" : $1;
1100             } elsif ( /^RC\s{3}(.*)/ ) {
1101 140 50       346 $com .= $com ? " $1" : $1;
1102             }
1103 1783         2513 $_ = $self->_readline;
1104             }
1105 30         82 return \@refs;
1106             }
1107              
1108              
1109             =head2 _read_swissprot_Species
1110              
1111             Title : _read_swissprot_Species
1112             Usage :
1113             Function: Reads the swissprot Organism species and classification
1114             lines.
1115             Able to deal with unconventional species names.
1116             Example : OS Unknown prokaryotic organism
1117             $genus = undef ; $species = Unknown prokaryotic organism
1118             Returns : A Bio::Species object
1119             Args :
1120              
1121             =cut
1122              
1123             sub _read_swissprot_Species {
1124 30     30   42 my( $self,$line ) = @_;
1125 30         32 my $org;
1126 30         46 local $_ = $line;
1127              
1128 30         32 my( $sub_species, $species, $genus, $common, $variant, $ncbi_taxid, $sci_name, $class_lines, $descr );
1129 30         47 my $osline = "";
1130 30         38 my $do_genus_check = 1;
1131 30         83 while ( defined $_ ) {
1132 183 100       403 last unless /^O[SCGX]/;
1133             # believe it or not, but OS may come multiple times -- at this time
1134             # we can't capture multiple species
1135 153 100 100     855 if (/^OS\s+(\S.+)/ && (! defined($sci_name))) {
    100 100        
    50          
    100          
1136 30 50       63 $osline .= " " if $osline;
1137 30         67 $osline .= $1;
1138 30 50       235 if ($osline =~ s/(,|, and|\.)$//) {
1139             # OS lines are usually like:
1140             # Homo sapiens (human)
1141             # where we have $sci_name followed by $descr (common name) in
1142             # brackets, but we can also have:
1143             # Venerupis (Ruditapes) philippinarum
1144             # where we have brackets but they don't indicate a $descr
1145 30 50       100 if ($osline =~ /[^\(\)]+\(.+\)[^\(\)]+$/) {
1146             #*** Danger! no idea if this will pick up some syntaxes for
1147             # common names as well)
1148 0         0 $sci_name = $osline;
1149 0         0 $sci_name =~ s/\.$//;
1150 0         0 $descr = '';
1151 0         0 $do_genus_check = 0;
1152             } else {
1153 30         125 ($sci_name, $descr) = $osline =~ /(\S[^\(]+)(.*)/;
1154             }
1155 30         83 $sci_name =~ s/\s+$//;
1156              
1157 30         122 while ($descr =~ /\(([^\)]+)\)/g) {
1158 13         25 my $item = $1;
1159             # strain etc may not necessarily come first (yes, swissprot
1160             # is messy)
1161 13 50 33     150 if ((! defined($variant)) &&
    50 33        
    50          
1162             (($item =~ /(^|[^\(\w])([Ss]train|isolate|serogroup|serotype|subtype|clone)\b/) ||
1163             ($item =~ /^(biovar|pv\.|type\s+)/))) {
1164 0         0 $variant = $item;
1165             } elsif ($item =~ s/^subsp\.\s+//) {
1166 0 0       0 if (! $sub_species) {
    0          
1167 0         0 $sub_species = $item;
1168             } elsif (! $variant) {
1169 0         0 $variant = $item;
1170             }
1171             } elsif (! defined($common)) {
1172             # we're only interested in the first common name
1173 13         17 $common = $item;
1174 13 50 33     68 if ((index($common, '(') >= 0) &&
1175             (index($common, ')') < 0)) {
1176 0         0 $common .= ')';
1177             }
1178             }
1179             }
1180             }
1181             } elsif (s/^OC\s+(\S.+)$//) {
1182 62         116 $class_lines .= $1;
1183             } elsif (/^OG\s+(.*)/) {
1184 0         0 $org = $1;
1185             } elsif (/^OX\s+(.*)/ && (! defined($ncbi_taxid))) {
1186 29         49 my $taxstring = $1;
1187             # we only keep the first one and ignore all others
1188 29 50       102 if ($taxstring =~ /NCBI_TaxID=([\w\d]+)/) {
1189 29         50 $ncbi_taxid = $1;
1190             } else {
1191 0         0 $self->throw("$taxstring doesn't look like NCBI_TaxID");
1192             }
1193             }
1194 153         254 $_ = $self->_readline;
1195             }
1196 30         93 $self->_pushback($_); # pushback the last line because we need it
1197              
1198 30 50       64 $sci_name || return;
1199              
1200             # if the organism belongs to taxid 32644 then no Bio::Species object.
1201 30 50       62 return if grep { $_ eq $sci_name } @Unknown_names;
  330         346  
1202              
1203             # Convert data in classification lines into classification array.
1204             # Remove trailing . then split on ';' or '.;' so that classification that is 2
1205             # or more words will still get matched, use map() to remove trailing/leading/intervening
1206             # spaces
1207 30         111 $class_lines=~s/\.\s*$//;
1208 30         257 my @class = map { s/^\s+//; s/\s+$//; s/\s{2,}/ /g; $_; } split /[;\.]*;/, $class_lines;
  298         377  
  298         249  
  298         226  
  298         348  
1209              
1210 30 50       130 if ($class[0] =~ /viruses/i) {
    50          
1211             # viruses have different OS/OC syntax
1212 0         0 my @virusnames = split(/\s+/, $sci_name);
1213 0 0       0 $species = (@virusnames > 1) ? pop(@virusnames) : '';
1214 0         0 $genus = join(" ", @virusnames);
1215 0         0 $sub_species = $descr;
1216             } elsif ($do_genus_check) {
1217             # do we have a genus?
1218 30         43 my $possible_genus = $class[-1];
1219 30 50       82 $possible_genus .= "|$class[-2]" if $class[-2];
1220 30 50       1054 if ($sci_name =~ /^($possible_genus)/) {
1221 30         80 $genus = $1;
1222 30         308 ($species) = $sci_name =~ /^$genus\s+(.+)/;
1223             } else {
1224 0         0 $species = $sci_name;
1225             }
1226             # is this organism of rank species or is it lower?
1227             # (doesn't catch everything, but at least the guess isn't dangerous)
1228 30 50 33     200 if ($species && $species =~ /subsp\.|var\./) {
1229 0         0 ($species, $sub_species) = $species =~ /(.+)\s+((?:subsp\.|var\.).+)/;
1230             }
1231             }
1232              
1233             # Bio::Species array needs array in Species -> Kingdom direction
1234 30 50       75 unless ($class[-1] eq $sci_name) {
1235 30         47 push(@class, $sci_name);
1236             }
1237 30         46 @class = reverse @class;
1238              
1239 30         174 my $taxon = Bio::Species->new();
1240 30         140 $taxon->scientific_name($sci_name);
1241 30         100 $taxon->classification(@class);
1242 30 100       94 $taxon->common_name($common) if $common;
1243 30 50       67 $taxon->sub_species($sub_species) if $sub_species;
1244 30 50       70 $taxon->organelle($org) if $org;
1245 30 100       124 $taxon->ncbi_taxid($ncbi_taxid) if $ncbi_taxid;
1246 30 50       67 $taxon->variant($variant) if $variant;
1247              
1248             # done
1249 30         118 return $taxon;
1250             }
1251              
1252             =head2 _filehandle
1253              
1254             Title : _filehandle
1255             Usage : $obj->_filehandle($newval)
1256             Function:
1257             Example :
1258             Returns : value of _filehandle
1259             Args : newvalue (optional)
1260              
1261              
1262             =cut
1263              
1264             # inherited from SeqIO.pm ! HL 05/11/2000
1265              
1266             =head2 _read_FTHelper_swissprot
1267              
1268             Title : _read_FTHelper_swissprot
1269             Usage : _read_FTHelper_swissprot(\$buffer)
1270             Function: reads the next FT key line
1271             Example :
1272             Returns : Bio::SeqIO::FTHelper object
1273             Args :
1274              
1275              
1276             =cut
1277              
1278             sub _read_FTHelper_swissprot {
1279 443     443   430 my ($self,$line ) = @_;
1280             # initial version implemented by HL 05/10/2000
1281             # FIXME this may not be perfect, so please review
1282             # lots of cleaning up by JES 2004/07/01, still may not be perfect =)
1283             # FTId now sepated from description as a qualifier
1284              
1285 443         471 local $_ = $line;
1286 443         315 my ($key, # The key of the feature
1287             $loc, # The location line from the feature
1288             $desc, # The descriptive text
1289             $ftid, # feature Id is like a qualifier but there can be only one of them
1290             );
1291 443 50       1592 if ( m/^FT\s{3}(\w+)\s+([\d\?\<]+)\s+([\d\?\>]+)\s*(.*)$/ox) {
1292 443         726 $key = $1;
1293 443         454 my $loc1 = $2;
1294 443         382 my $loc2 = $3;
1295 443         557 $loc = "$loc1..$loc2";
1296 443 100 66     1584 if ($4 && (length($4) > 0)) {
1297 384         379 $desc = $4;
1298 384         502 chomp($desc);
1299             } else {
1300 59         71 $desc = "";
1301             }
1302             }
1303              
1304 443   66     707 while ( defined($_ = $self->_readline) && /^FT\s{20,}(\S.*)$/ ) {
1305 312         447 my $continuation_line = $1;
1306 312 100       797 if ( $continuation_line =~ /.FTId=(.*)\./ ) {
    50          
1307 272         321 $ftid=$1;
1308             }
1309             elsif ( $desc) {
1310 40         85 $desc .= " $continuation_line";
1311             } else {
1312 0         0 $desc = $continuation_line;
1313             }
1314 312         531 chomp $desc;
1315             }
1316 443         812 $self->_pushback($_);
1317 443 50       697 unless( $key ) {
1318             # No feature key. What's this?
1319 0         0 $self->warn("No feature key in putative feature table line: $line");
1320 0         0 return;
1321             }
1322              
1323             # Make the new FTHelper object
1324 443         926 my $out = Bio::SeqIO::FTHelper->new(-verbose => $self->verbose());
1325 443         736 $out->key($key);
1326 443         590 $out->loc($loc);
1327              
1328             # store the description if there is one
1329 443 100 66     1263 if ( $desc && length($desc) ) {
1330 384         1013 $desc =~ s/\.$//;
1331 384         338 push(@{$out->field->{"description"}}, $desc);
  384         634  
1332             }
1333             # Store the qualifier i.e. FTId
1334 443 100       720 if ( $ftid ) {
1335 272         200 push(@{$out->field->{"FTId"}}, $ftid);
  272         371  
1336             }
1337 443         658 return $out;
1338             }
1339              
1340              
1341             =head2 _write_line_swissprot
1342              
1343             Title : _write_line_swissprot
1344             Usage :
1345             Function: internal function
1346             Example :
1347             Returns :
1348             Args :
1349              
1350              
1351             =cut
1352              
1353             sub _write_line_swissprot{
1354 0     0   0 my ($self,$pre1,$pre2,$line,$length) = @_;
1355              
1356 0 0       0 $length || $self->throw( "Miscalled write_line_swissprot without length. Programming error!");
1357 0         0 my $subl = $length - length $pre2;
1358 0         0 my $linel = length $line;
1359 0         0 my $i;
1360              
1361 0         0 my $sub = substr($line,0,$length - length $pre1);
1362              
1363 0         0 $self->_print( "$pre1$sub\n");
1364              
1365 0         0 for ($i= ($length - length $pre1);$i < $linel;) {
1366 0         0 $sub = substr($line,$i,($subl));
1367 0         0 $self->_print( "$pre2$sub\n");
1368 0         0 $i += $subl;
1369             }
1370              
1371             }
1372              
1373             =head2 _write_line_swissprot_regex
1374              
1375             Title : _write_line_swissprot_regex
1376             Usage :
1377             Function: internal function for writing lines of specified
1378             length, with different first and the next line
1379             left hand headers and split at specific points in the
1380             text
1381             Example :
1382             Returns : nothing
1383             Args : file handle, first header, second header, text-line, regex for line breaks, total line length
1384              
1385              
1386             =cut
1387              
1388             sub _write_line_swissprot_regex {
1389 44     44   60 my ($self,$pre1,$pre2,$line,$regex,$length) = @_;
1390              
1391             #print STDOUT "Going to print with $line!\n";
1392              
1393 44 50       69 $length || $self->throw( "Miscalled write_line_swissprot without length. Programming error!");
1394              
1395 44 50       65 if ( length $pre1 != length $pre2 ) {
1396 0         0 $self->warn( "len 1 is ". length ($pre1) . " len 2 is ". length ($pre2) . "\n");
1397 0         0 $self->throw( "Programming error - cannot called write_line_swissprot_regex with different length \npre1 ($pre1) and \npre2 ($pre2) tags!");
1398             }
1399              
1400 44         48 my $subl = $length - (length $pre1) -1 ;
1401              
1402 44         32 my $first_line = 1;
1403 44         366 while ($line =~ m/(.{1,$subl})($regex)/g) {
1404 54         105 my $s = $1.$2;
1405 54 100 100     240 $s =~ s/([\w\.])#(\w)/$1 $2/g # remove word wrap protection char '#'
1406             if $pre1 eq "RA " or $pre1 eq "KW ";
1407             # remove annoying extra spaces at the end of the wrapped lines
1408 54 100       99 substr($s, -1, 1, '') if substr($s, -1, 1) eq ' ';
1409 54 100       60 if ($first_line) {
1410 44         108 $self->_print( "$pre1$s\n");
1411 44         176 $first_line = 0;
1412             } else {
1413 10         29 $self->_print( "$pre2$s\n");
1414             }
1415              
1416             }
1417             }
1418              
1419             =head2 _post_sort
1420              
1421             Title : _post_sort
1422             Usage : $obj->_post_sort($newval)
1423             Function:
1424             Returns : value of _post_sort
1425             Args : newvalue (optional)
1426              
1427              
1428             =cut
1429              
1430             sub _post_sort{
1431 3     3   7 my $obj = shift;
1432 3 50       10 if ( @_ ) {
1433 0         0 my $value = shift;
1434 0         0 $obj->{'_post_sort'} = $value;
1435             }
1436 3         8 return $obj->{'_post_sort'};
1437              
1438             }
1439              
1440             =head2 _show_dna
1441              
1442             Title : _show_dna
1443             Usage : $obj->_show_dna($newval)
1444             Function:
1445             Returns : value of _show_dna
1446             Args : newvalue (optional)
1447              
1448              
1449             =cut
1450              
1451             sub _show_dna{
1452 26     26   38 my $obj = shift;
1453 26 100       72 if ( @_ ) {
1454 23         23 my $value = shift;
1455 23         39 $obj->{'_show_dna'} = $value;
1456             }
1457 26         42 return $obj->{'_show_dna'};
1458              
1459             }
1460              
1461             =head2 _id_generation_func
1462              
1463             Title : _id_generation_func
1464             Usage : $obj->_id_generation_func($newval)
1465             Function:
1466             Returns : value of _id_generation_func
1467             Args : newvalue (optional)
1468              
1469              
1470             =cut
1471              
1472             sub _id_generation_func{
1473 4     4   7 my $obj = shift;
1474 4 50       11 if ( @_ ) {
1475 0         0 my $value = shift;
1476 0         0 $obj->{'_id_generation_func'} = $value;
1477             }
1478 4         12 return $obj->{'_id_generation_func'};
1479              
1480             }
1481              
1482             =head2 _ac_generation_func
1483              
1484             Title : _ac_generation_func
1485             Usage : $obj->_ac_generation_func($newval)
1486             Function:
1487             Returns : value of _ac_generation_func
1488             Args : newvalue (optional)
1489              
1490              
1491             =cut
1492              
1493             sub _ac_generation_func{
1494 4     4   6 my $obj = shift;
1495 4 50       15 if ( @_ ) {
1496 0         0 my $value = shift;
1497 0         0 $obj->{'_ac_generation_func'} = $value;
1498             }
1499 4         31 return $obj->{'_ac_generation_func'};
1500              
1501             }
1502              
1503             =head2 _sv_generation_func
1504              
1505             Title : _sv_generation_func
1506             Usage : $obj->_sv_generation_func($newval)
1507             Function:
1508             Returns : value of _sv_generation_func
1509             Args : newvalue (optional)
1510              
1511              
1512             =cut
1513              
1514             sub _sv_generation_func{
1515 0     0   0 my $obj = shift;
1516 0 0       0 if ( @_ ) {
1517 0         0 my $value = shift;
1518 0         0 $obj->{'_sv_generation_func'} = $value;
1519             }
1520 0         0 return $obj->{'_sv_generation_func'};
1521              
1522             }
1523              
1524             =head2 _kw_generation_func
1525              
1526             Title : _kw_generation_func
1527             Usage : $obj->_kw_generation_func($newval)
1528             Function:
1529             Returns : value of _kw_generation_func
1530             Args : newvalue (optional)
1531              
1532              
1533             =cut
1534              
1535             sub _kw_generation_func{
1536 4     4   8 my $obj = shift;
1537 4 50       11 if ( @_ ) {
1538 0         0 my $value = shift;
1539 0         0 $obj->{'_kw_generation_func'} = $value;
1540             }
1541 4         29 return $obj->{'_kw_generation_func'};
1542              
1543             }
1544              
1545             1;