Bio/SeqIO/interpro.pm

   1 #
   2 # BioPerl module for interpro
   3 # You may distribute this module under the same terms as perl itself
   4
   5 # POD documentation - main docs before the code
   6
   7 =head1 NAME
   8
   9 Bio::SeqIO::interpro - InterProScan XML input/output stream
  10
  11 =head1 SYNOPSIS
  12
  13   # do not call this module directly, use Bio::SeqIO
  14
  15   use strict;
  16   use Bio::SeqIO;
  17
  18   my $io = Bio::SeqIO->new(-format => "interpro",
  19                            -file   => $interpro_file);
  20
  21   while (my $seq = $io->next_seq) {
  22     # use the Sequence object
  23   }
  24
  25 =head1 DESCRIPTION
  26
  27 L<Bio::SeqIO::interpro> will parse Interpro scan XML (version 1.2) and
  28 create L<Bio::SeqFeature::Generic> objects based on the contents of the
  29 XML document.
  30
  31 L<Bio::SeqIO::interpro> will also attach the annotation given in the XML
  32 file to the L<Bio::SeqFeature::Generic> objects that it creates.
  33
  34 =head1 FEEDBACK
  35
  36 =head2 Mailing Lists
  37
  38 User feedback is an integral part of the evolution of this and other
  39 Bioperl modules. Send your comments and suggestions preferably to
  40 the Bioperl mailing list.  Your participation is much appreciated.
  41
  42   bioperl-l@bioperl.org                  - General discussion
  43   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  44
  45 =head2 Support
  46
  47 Please direct usage questions or support issues to the mailing list:
  48
  49 I<bioperl-l@bioperl.org>
  50
  51 rather than to the module maintainer directly. Many experienced and
  52 reponsive experts will be able look at the problem and quickly
  53 address it. Please include a thorough description of the problem
  54 with code and data examples if at all possible.
  55
  56 =head2 Reporting Bugs
  57
  58 Report bugs to the Bioperl bug tracking system to help us keep track
  59 of the bugs and their resolution. Bug reports can be submitted via
  60 the web:
  61
  62   https://github.com/bioperl/bioperl-live/issues
  63
  64 =head1 AUTHOR - Jared Fox
  65
  66 Email jaredfox@ucla.edu
  67
  68 =head1 CONTRIBUTORS
  69
  70 Allen Day allenday@ucla.edu
  71
  72 =head1 APPENDIX
  73
  74 The rest of the documentation details each of the object methods.
  75 Internal methods are usually preceded with a _
  76
  77 =cut
  78
  79 # Let the code begin...
  80
  81 package Bio::SeqIO::interpro;
  82 use strict;
  83 use Bio::SeqFeature::Generic;
  84 use XML::DOM;
  85 use XML::DOM::XPath;
  86 use Bio::Seq::SeqFactory;
  87 use Bio::Annotation::Collection;
  88 use Bio::Annotation::DBLink;
  89 use base qw(Bio::SeqIO);
  90
  91 my $idcounter = {};  # Used to generate unique id values
  92 my $nvtoken = ": ";  # The token used if a name/value pair has to be stuffed
  93                      # into a single line
  94
  95 =head2 next_seq
  96
  97  Title   : next_seq
  98  Usage   : my $seqobj = $stream->next_seq
  99  Function: Retrieves the next sequence from a SeqIO::interpro stream.
 100  Returns : A Bio::Seq::RichSeq object
 101  Args    :
 102
 103 =cut
 104
 105 sub next_seq {
 106         my $self = shift;
 107         my ($desc);
 108         my $bioSeq = $self->_sequence_factory->create(-verbose =>$self->verbose());
 109
 110         my $zinc = "(\"zincins\")";
 111         my $wing = "\"Winged helix\"";
 112         my $finger = "\"zinc finger\"";
 113
 114         my $xml_fragment = undef;
 115         while(my $line = $self->_readline()){
 116
 117                 my $where = index($line, $zinc);
 118                 my $wherefinger = index($line, $finger);
 119                 my $finishedline = $line;
 120                 my $wingwhere = index($line, $wing);
 121
 122                 # the interpro XML is not fully formed, so we need to convert the
 123                 # extra double quotes and ampersands into appropriate XML chracter codes
 124                 if($where > 0){
 125                         my @linearray = split /$zinc/, $line;
 126                         $finishedline = join "&quot;zincins&quot;", $linearray[0], $linearray[2];
 127                 }
 128                 if(index($line, "&") > 0){
 129                         my @linearray = split /&/, $line;
 130                         $finishedline = join "&amp;", $linearray[0], $linearray[1];
 131                 }
 132                 if($wingwhere > 0){
 133                         my @linearray = split /$wing/, $line;
 134                         $finishedline = join "&quot;Winged helix&quot;", $linearray[0], $linearray[1];
 135                 }
 136
 137                 $xml_fragment .= $finishedline;
 138                 last if $finishedline =~ m!</protein>!;
 139         }
 140         # Match <protein> but not other similar elements like <protein-matches>
 141         return unless $xml_fragment =~ /<protein[\s>]/;
 142
 143         $self->_parse_xml($xml_fragment);
 144
 145         my $dom = $self->_dom;
 146
 147         my ($protein_node) = $dom->findnodes('/protein');
 148         my @interproNodes = $protein_node->findnodes('/protein/interpro');
 149         my @DBNodes = $protein_node->findnodes('/protein/interpro/match');
 150         for(my $interpn=0; $interpn<scalar(@interproNodes); $interpn++){
 151                 my $ipnlevel = join "", "/protein/interpro[", $interpn + 1, "]";
 152                 my @matchNodes = $protein_node->findnodes($ipnlevel);
 153                 for(my $match=0; $match<scalar(@matchNodes); $match++){
 154                         my $matlevel = join "", "/protein/interpro[", $interpn+1, "]/match[",
 155                           $match+1, "]/location";
 156                         my @locNodes = $protein_node->findnodes($matlevel);
 157                         my $class_level = join "", "/protein/interpro[",$interpn+1, "]/classification";
 158                         my @goNodes = $protein_node->findnodes($class_level);
 159                         my @seqFeatures = map { Bio::SeqFeature::Generic->new(
 160                   -start => $_->getAttribute('start'),
 161                                                 -end => $_->getAttribute('end'),
 162                   -score => $_->getAttribute('score'),
 163                   -source_tag => 'IPRscan',
 164                   -primary_tag => 'region',
 165                   -display_name => $interproNodes[$interpn]->getAttribute('name'),
 166                   -seq_id => $protein_node->getAttribute('id') ),
 167                                         } @locNodes;
 168                         foreach my $seqFeature (@seqFeatures){
 169                                 $bioSeq->add_SeqFeature($seqFeature);
 170
 171                                 my $annotation1 = Bio::Annotation::DBLink->new;
 172                                 $annotation1->database($matchNodes[$match]->getAttribute('dbname'));
 173                                 $annotation1->primary_id($matchNodes[$match]->getAttribute('id'));
 174                                 $annotation1->comment($matchNodes[$match]->getAttribute('name'));
 175                                 $seqFeature->annotation->add_Annotation('dblink',$annotation1);
 176
 177                                 my $annotation2 = Bio::Annotation::DBLink->new;
 178                                 $annotation2->database('INTERPRO');
 179                                 $annotation2->primary_id($interproNodes[$interpn]->getAttribute('id'));
 180                                 $annotation2->comment($interproNodes[$interpn]->getAttribute('name'));
 181                                 $seqFeature->annotation->add_Annotation('dblink',$annotation2);
 182
 183                                 # Bug 1908 (enhancement)
 184                                 my $annotation3  = Bio::Annotation::DBLink->new;
 185                                 $annotation3->database($DBNodes[$interpn]->getAttribute('dbname'));
 186                                 $annotation3->primary_id($DBNodes[$interpn]->getAttribute('id'));
 187                                 $annotation3->comment($DBNodes[$interpn]->getAttribute('name'));
 188                                 $seqFeature->annotation->add_Annotation('dblink',$annotation3);
 189                                 # need to put in the go annotation here!
 190                                  foreach my $g (@goNodes)
 191                                  {
 192                                      my $goid = $g->getAttribute('id');
 193                                      my $go_annotation = Bio::Annotation::DBLink->new;
 194                                      $go_annotation->database('GO');
 195                                      $go_annotation->primary_id($goid);
 196                                      $go_annotation->comment($goid);
 197                                      $seqFeature->annotation->add_Annotation('dblink', $go_annotation);
 198                                  }
 199                         }
 200                 }
 201         }
 202         my $accession = $protein_node->getAttribute('id');
 203         my $displayname = $protein_node->getAttribute('name');
 204         $bioSeq->accession($accession);
 205         $bioSeq->display_name($displayname);
 206         return $bioSeq;
 207 }
 208
 209 =head2 _initialize
 210
 211  Title   : _initialize
 212  Usage   :
 213  Function:
 214  Returns :
 215  Args    :
 216
 217 =cut
 218
 219 sub _initialize {
 220   my($self,@args) = @_;
 221
 222   $self->SUPER::_initialize(@args);
 223   # hash for functions for decoding keys.
 224   $self->{'_func_ftunit_hash'} = {};
 225
 226   my %param = @args;  # From SeqIO.pm
 227   @param{ map { lc $_ } keys %param } = values %param; # lowercase keys
 228
 229   my $line = undef;
 230   # fast forward to first <protein/> record.
 231   while($line = $self->_readline()){
 232     # Match <protein> but not other similar elements like <protein-matches>
 233     if($line =~ /<protein[\s>]/){
 234       $self->_pushback($line);
 235       last;
 236     }
 237   }
 238
 239   $self->_xml_parser( XML::DOM::Parser->new() );
 240
 241   $self->_sequence_factory( Bio::Seq::SeqFactory->new
 242                            ( -verbose => $self->verbose(),
 243                              -type => 'Bio::Seq::RichSeq'))
 244     if ( ! defined $self->sequence_factory );
 245 }
 246
 247 =head2 _sequence_factory
 248
 249  Title   : _sequence_factory
 250  Usage   :
 251  Function:
 252  Returns :
 253  Args    :
 254
 255 =cut
 256
 257 sub _sequence_factory {
 258   my $self = shift;
 259   my $val = shift;
 260
 261   $self->{'sequence_factory'} = $val if defined($val);
 262   return $self->{'sequence_factory'};
 263 }
 264
 265 =head2 _xml_parser
 266
 267  Title   : _xml_parser
 268  Usage   :
 269  Function:
 270  Returns :
 271  Args    :
 272
 273 =cut
 274
 275 sub _xml_parser {
 276   my $self = shift;
 277   my $val = shift;
 278
 279   $self->{'xml_parser'} = $val if defined($val);
 280   return $self->{'xml_parser'};
 281 }
 282
 283 =head2 _parse_xml
 284
 285  Title   : _parse_xml
 286  Usage   :
 287  Function:
 288  Returns :
 289  Args    :
 290
 291 =cut
 292
 293 sub _parse_xml {
 294   my ($self,$xml) = @_;
 295   $self->_dom( $self->_xml_parser->parse($xml) );
 296   return 1;
 297 }
 298
 299 =head2 _dom
 300
 301  Title   : _dom
 302  Usage   :
 303  Function:
 304  Returns :
 305  Args    :
 306
 307 =cut
 308
 309 sub _dom {
 310   my $self = shift;
 311   my $val = shift;
 312
 313   $self->{'dom'} = $val if defined($val);
 314   return $self->{'dom'};
 315 }
 316
 317 1;
 318
 319 __END__