Bio/Tools/Analysis/Protein/GOR4.pm

   1 # $Id: GOR4.pm,v 1.0 2003/07/ 11
   2 #
   3 # BioPerl module for Bio::Tools::Analysis::Protein::GOR4
   4 #
   5 # Copyright Richard Adams
   6 #
   7 # You may distribute this module under the same terms as perl itself
   8
   9 # POD documentation - main docs before the code
  10
  11 =head1  NAME
  12
  13 Bio::Tools::Analysis::Protein::GOR4 - a wrapper around GOR4 protein
  14 secondary structure prediction server
  15
  16 =head1  SYNOPSIS
  17
  18   use Bio::Tools::Analysis::Protein::GOR4;
  19   #get a Bio::Seq or Bio::PrimarySeq
  20   use Bio::PrimarySeq;
  21   $seq = Bio::PrimarySeq->new
  22     (-seq=>'IKLCVHHJHJHJHJHJHJHNLAILAKAHLIELALAL',
  23      -primary_id=>'test'); # a Bio::PrimarySeqI object
  24
  25   my $gor4 = Bio::Tools::Analysis::Protein::GOR4->new (-seq=>$seq);
  26   $gor4->run;
  27   print $gor4->result;# #raw text to standard error
  28
  29 =head1  DESCRIPTION
  30
  31 A module to remotely retrieve predictions of protein secondary
  32 structure.  Each residue in the protein receives a score representing
  33 the likelihood of existing in each of three different states (helix,
  34 coil or sheet), e.g.,
  35
  36   my $analysis_object = Bio::Tools::SimpleAnalysis::Protein::GOR4->
  37       new(-seq => $seq);
  38
  39 creates a new object
  40
  41   $analysis_object->run;
  42
  43 submits the query to the server and obtains raw text output
  44
  45 Given an amino acid sequence the results can be obtained in 4 formats,
  46 determined by the argument to the result method
  47
  48 =over 4
  49
  50 =item 1
  51
  52 The raw text of the program output
  53
  54   my $rawdata = $analysis_object->result;
  55
  56 =item 2
  57
  58 An reference to an array of hashes of scores for each state and the
  59 assigned state.
  60
  61   my $data_ref = $analysis_object->result('parsed');
  62   print "score for helix at residue 2 is $data_ref->[1]{'helix'}\n";
  63   print "predicted struc  at residue 2 is $data_ref->[1]{'struc}\n";
  64
  65 =item 3
  66
  67 An array of Bio::SeqFeature::Generic objects where each feature is a
  68 predicted unit of secondary structure. Only stretches of helix/sheet
  69 predictions for longer than 4 residues are defined as helices. See
  70 Bio::Tools::Analysis::Domcut.pm for examples of how to add sequence
  71 features.
  72
  73   my @fts = $analysis_object->result(Bio::SeqFeatureI);
  74   for my $ft (@fts) {
  75       print " From ",  $ft->start, " to  ",$ft->end, " struc: " ,
  76              ($ft->each_tag_value('type'))[0]  ,"\n";
  77   }
  78
  79 =item 4
  80
  81 A Bio::Seq::Meta::Array implementing sequence.
  82
  83 This is a Bio::Seq object that can also hold data about each residue
  84 in the sequence In this case, the sequence can be associated with a
  85 single array of GOR4 prediction scores.  e.g.,
  86
  87   my $meta_sequence = $analysis_object->result('all');
  88   print "helix scores from residues 10-20 are ",
  89       $meta_sequence->named_submeta_text("GOR4_helix",10,20), "\n";
  90
  91 Meta sequence names are : GOR4_helix, GOR4_sheet, GOR4_coil,
  92 GOR4_struc, representing the scores for each residue.
  93
  94 Many methods common to all analyses are inherited from
  95 Bio::Tools::Analysis::SimpleAnalysisBase.
  96
  97 =back
  98
  99 =head1 SEE ALSO
 100
 101 L<Bio::SimpleAnalysisI>,
 102 L<Bio::Tools::Analysis::SimpleAnalysisBase>,
 103 L<Bio::Seq::Meta::Array>,
 104 L<Bio::WebAgent>
 105
 106 =head1 FEEDBACK
 107
 108 =head2 Mailing Lists
 109
 110 User feedback is an integral part of the evolution of this and other
 111 Bioperl modules. Send your comments and suggestions preferably to one
 112 of the Bioperl mailing lists.  Your participation is much appreciated.
 113
 114   bioperl-l@bioperl.org                  - General discussion
 115   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 116
 117 =head2 Support
 118
 119 Please direct usage questions or support issues to the mailing list:
 120
 121 I<bioperl-l@bioperl.org>
 122
 123 rather than to the module maintainer directly. Many experienced and
 124 reponsive experts will be able look at the problem and quickly
 125 address it. Please include a thorough description of the problem
 126 with code and data examples if at all possible.
 127
 128 =head2 Reporting Bugs
 129
 130 Report bugs to the Bioperl bug tracking system to help us keep track
 131 the bugs and their resolution.  Bug reports can be submitted via the
 132 web:
 133
 134   https://github.com/bioperl/bioperl-live/issues
 135
 136 =head1 AUTHORS
 137
 138 Richard Adams, Richard.Adams@ed.ac.uk,
 139
 140 =head1 APPENDIX
 141
 142 The rest of the documentation details each of the object
 143 methods. Internal methods are usually preceded with a _
 144
 145 =cut
 146
 147 use strict;
 148
 149 package Bio::Tools::Analysis::Protein::GOR4;
 150
 151 use IO::String;
 152 use Bio::SeqIO;
 153 use HTTP::Request::Common qw(POST);
 154 use Bio::SeqFeature::Generic;
 155 use Bio::Seq::Meta::Array;
 156 $ENV{PERL_LWP_SSL_VERIFY_HOSTNAME} = 0;
 157
 158 use base qw(Bio::Tools::Analysis::SimpleAnalysisBase);
 159
 160 use constant MIN_STRUC_LEN => 3;
 161 my $URL = 'https://npsa-prabi.ibcp.fr/cgi-bin/secpred_sopma.pl';
 162 my $ANALYSIS_NAME = 'GOR4';
 163 my $ANALYSIS_SPEC = {name => 'Gor4', type => 'Protein'};
 164 my $INPUT_SPEC    = [
 165                      {mandatory =>'true',
 166                       type      => 'Bio::PrimarySeqI',
 167                       'name'    => 'seq',
 168                   },
 169                  ];
 170 my  $RESULT_SPEC =
 171     {
 172      ''                 => 'bulk',              # same as undef
 173      'Bio::SeqFeatureI' => 'ARRAY of Bio::SeqFeature::Generic',
 174      raw                => '[ {struc =>, helix=> ,sheet=>, coil=>}]',
 175      meta                => 'Bio::Seq::Meta::Array object',
 176     };
 177
 178 =head2 result
 179
 180  Name    : result
 181  Usage   : $job->result (...)
 182  Returns : a result created by running an analysis
 183  Args    : see keys of $RESULT_SPEC
 184
 185 The method returns a result of an executed job. If the job was
 186 terminated by an error the result may contain an error message instead
 187 of the real data.
 188
 189 This implementation returns differently processed data depending on
 190 argument:
 191
 192 =over 3
 193
 194 =item undef
 195
 196 Returns the raw ASCII data stream but without HTML tags
 197
 198 =item 'Bio::SeqFeatureI'
 199
 200 The argument string defines the type of bioperl objects returned in an
 201 array.  The objects are L<Bio::SeqFeature::Generic>.  Feature primary
 202 tag is "2ary".  Feature tags are "type" (which can be helix, sheet or
 203 coil) "method" (GOR4).
 204
 205 =item 'parsed'
 206
 207 Array of hash references of { helix =E<gt>, sheet =E<gt> , coil =E<gt> , struc=E<gt>}.
 208
 209 =item 'meta'
 210
 211 A Bio::Seq::Meta::Array object. Scores can be accessed using methods
 212 from this class. Meta sequence names are GOR4_helix, GOR4_sheet,
 213 GOR4_coil, GOR4_struc.
 214
 215
 216 =back
 217
 218
 219 =cut
 220
 221
 222 sub result {
 223     my ($self,$value) = @_;
 224
 225     my @scores;
 226     my @fts;
 227
 228     if ($value ) {
 229         #parse into basic raw form, store this as well as '_result'
 230         if (!exists($self->{'_parsed'}) ) {
 231             my $result = IO::String->new($self->{'_result'});
 232             while (my $line = <$result>) {
 233                 next unless $line =~ /^\w\s/; # or for sopma/hnn  /^[A-Z]\s/
 234                 $line =~/(\w)\s+(\d+)\s+(\d+)\s+(\d+)/; # or for so
 235                 push @scores, { struc => $1,
 236                                 helix => $2,
 237                                 sheet => $3,
 238                                 coil  => $4,
 239                               };
 240             }
 241             $self->{'_parsed'} = \@scores;
 242         }
 243         if ($value eq 'Bio::SeqFeatureI') {
 244             $self->_get_2ary_coords();
 245             for my $type (keys %{$self->{'_parsed_coords'}} ) {
 246                 next if $type =~  /\w{2,}/; #if not H,C,E or T
 247                 for my $loc (@{$self->{'_parsed_coords'}{$type}} ) {
 248                     push @fts, Bio::SeqFeature::Generic->new
 249                         (-start   => $loc->{'start'},
 250                          -end     => $loc->{'end'},
 251                          -source  => 'GOR4',
 252                          -primary => 'Region',
 253                          -tag => {
 254                                   type => $type,
 255                                   method => $self->analysis_name,
 256                                  });
 257                 }               #end of array of strucs of type
 258             }                   # end of all 2nd struc elements
 259             delete $self->{'_parsed_coords'}; #remove temp data
 260             return @fts;
 261         }                       #endif BioSeqFeature
 262
 263         elsif ($value eq 'meta') {
 264             #1st of all make 3 or 4 arrays of scores for each type from column data
 265             my %type_scores;
 266             for my $aa (@{$self->{'_parsed'}}) {
 267                 push @{$type_scores{'struc'}}, $aa->{'struc'};
 268                 push @{$type_scores{'helix'}}, $aa->{'helix'};
 269                 push @{$type_scores{'sheet'}}, $aa->{'sheet'};
 270                 push @{$type_scores{'coil'}}, $aa->{'coil'};
 271             }
 272
 273                         ## bless if necessary ##
 274                         if (!$self->seq->isa("Bio::Seq::Meta::Array")){
 275                 bless ($self->seq, "Bio::Seq::Meta::Array");
 276                                 }
 277             $self->seq->isa("Bio::Seq::MetaI")
 278                 || $self->throw("$self is not a Bio::Seq::MetaI");
 279             $Bio::Seq::Meta::Array::DEFAULT_NAME = 'GOR4_struc';
 280
 281             ## now make meta_Sequence
 282             for my $struc_type (keys %type_scores) {
 283                 my $meta_name = "GOR4". "_" . "$struc_type";
 284                 my @meta = map{$_->{$struc_type}} @{$self->{'_parsed'}};
 285                 if (grep{$_ eq $meta_name}$self->seq->meta_names ) {
 286                     $self->warn ("$meta_name already exists , not overwriting!");
 287                     next;
 288                 }
 289                 $self->seq->named_meta($meta_name,\@meta );
 290             }
 291             # return  seq array object implementing meta sequence #
 292             return $self->seq;
 293
 294         }
 295                 else  {
 296            return $self->{'_parsed'};
 297          }
 298     }  #endif ($value)
 299
 300     #return raw result if no return fomrt stated
 301     return $self->{'_result'};
 302 }
 303
 304 sub _get_2ary_coords {
 305     #helper sub for result;
 306     ##extracts runs of structure > MIN_STRUC_LENresidues or less if Turn:
 307     #i.e., helical prediction for 1 residue isn't very meaningful...
 308     ## and poulates array of hashes with start/end values.
 309     ##keys of $Result are 'H' 'T' 'C' 'E'.
 310     #could be put into a secondary base class if need be
 311     my ($self) = @_;
 312     my @prot = @{$self->{'_parsed'}};
 313     my %Result;
 314     for (my $index = 0; $index <= $#prot; $index++) {
 315
 316         my $type        = $prot[$index]{'struc'};
 317         next unless $type =~ /[HTCE]/;
 318         my $length = 1;
 319         for (my $j = $index + 1; $j <= $#prot; $j++) {
 320             my $test = $prot[$j];
 321             if ($test->{'struc'} eq $type) {
 322                 $length++;
 323             } elsif (  $length > MIN_STRUC_LEN  ||
 324                        ($length <= MIN_STRUC_LEN && $type eq 'T') ) {
 325                 push @{$Result{$type}}, {start => $index + 1 ,  end => $j};
 326                 $index += $length -1;
 327                 last;
 328             } else {
 329                 $index += $length - 1;
 330                 last;
 331             }
 332         }
 333     }
 334     $self->{'_parsed_coords'} = \%Result; #temp assignment
 335 }
 336
 337 sub _init {
 338     my $self = shift;
 339     $self->url($URL);
 340     $self->{'_ANALYSIS_SPEC'} =$ANALYSIS_SPEC;
 341     $self->{'_INPUT_SPEC'} =$INPUT_SPEC;
 342     $self->{'_RESULT_SPEC'} =$RESULT_SPEC;
 343     $self->{'_ANALYSIS_NAME'} =$ANALYSIS_NAME;
 344     return $self;
 345 }
 346
 347
 348 sub  _run {
 349     my $self  = shift;
 350     $self->delay(1);
 351     # delay repeated calls by default by 3 sec, set delay() to change
 352     $self->sleep;
 353     $self->status('TERMINATED_BY_ERROR');
 354     my $request = POST $self->url,
 355         Content_Type => 'form-data',
 356             Content  => [title => "",
 357                          notice => $self->seq->seq,
 358                          ali_width => 70,
 359                         ];
 360
 361     my $content = $self->request($request);
 362     my $text = $content->content;
 363     return unless $text;
 364     my ($next) = $text =~ /Prediction.*?=(.*?)>/;
 365     return unless $next;
 366     my $out = 'http://npsa-pbil.ibcp.fr/'.$next;
 367     my $req2 = HTTP::Request->new(GET=>$out);
 368     my $resp2 = $self->request($req2);
 369     $self->status('COMPLETED') if $resp2 ne '';
 370     $self->{'_result'} = $resp2->content;
 371 }
 372
 373
 374 1;