Bio/Tools/Analysis/Protein/Sopma.pm

   1 # $Id: Sopma.pm,v 1.0 2003/07/ 11
   2 #
   3 # BioPerl module for Bio::Tools::Analysis::Protein::Sopma
   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
  12 =head1 NAME
  13
  14 Bio::Tools::Analysis::Protein::Sopma - a wrapper around the
  15 Sopma protein secondary structure prediction server
  16
  17 =head1  SYNOPSIS
  18
  19   use Bio::Tools::Analysis::Protein::Sopma;
  20   #get a Bio::Seq or Bio::PrimarySeq
  21   my $seq;
  22
  23   my $sopma = Bio::Tools::Analysis::Protein::Sopma->new
  24       (-seq=>$seq, states=>4);
  25   $sopma->run;
  26   print $sopma->result;# #raw text to standard error
  27
  28 =head1  DESCRIPTION
  29
  30 A module to remotely retrieve predictions of protein secondary
  31 structure.  Each residue in the protein receives a score representing
  32 the likelihood of existing in each of four different states (helix,
  33 coil, turn or sheet), e.g.,
  34
  35   my $analysis_object = Bio::Tools::SimpleAnalysis::Protein::Sopma->new
  36       ( -seq          => $seq,
  37         -states       => 4,
  38         -window_width => 15,
  39       );
  40
  41 creates a new object.  Compulsory argument -seq.  Optional arguments
  42 -states, -window_width,-similarity_threshold. These arguments can also be
  43 set by direct methods , e.g.,
  44
  45   $analysis_object->states(4);
  46   $analysis_object->run;
  47
  48 submits the query to the server and obtains raw text output. Given an
  49 amino acid sequence the results can be obtained in 4 formats,
  50 determined by the argument to the result method:
  51
  52 =over 4
  53
  54 =item 1
  55
  56 The raw text of the program output.
  57
  58   my $rawdata = $analysis_object->result;
  59
  60 =item 2
  61
  62 A reference to an array of hashes of scores for each state and the
  63 assigned state.
  64
  65   my $data_ref = $analysis_object->result('parsed');
  66   print "score for helix at residue 2 is $data_ref->[1]{'helix'}\n";
  67   print "predicted struc  at residue 2 is $data_ref->[1]{'struc}\n";
  68
  69 Hash keys are 'helix', 'struc', 'sheet', 'coil', 'turn'.
  70
  71 =item 3
  72
  73 An array of Bio::SeqFeature::Generic objects where each feature is a
  74 predicted unit of secondary structure. Only stretches of helix/sheet
  75 predictions for longer than 4 residues are defined as helices/sheets.
  76
  77   my @fts = $analysis_object->result(Bio::SeqFeatureI);
  78   for my $ft (@fts) {
  79       print " From ",  $ft->start, " to  ",$ft->end, " struc: " ,
  80              ($ft->each_tag_value('type'))[0]  ,"\n";
  81   }
  82
  83 =item 4
  84
  85 A Bio::Seq::Meta::Array implementing sequence.
  86
  87 This is a Bio::Seq object that can also hold data about each residue
  88 in the sequence.  In this case, the sequence can be associated with a
  89 arrays of Sopma prediction scores.  e.g.,
  90
  91   my $meta_sequence = $analysis_object->result('meta');
  92   print "scores from residues 10 -20 are ",
  93       $meta_sequence->named_submeta_text("Sopma_helix",10,20), "\n";
  94
  95 Meta sequence names are : Sopma_helix, Sopma_sheet, Sopma_turn,
  96 Sopma_coil, Sopma_struc, representing the scores for each residue.
  97
  98 Many methods common to all analyses are inherited from
  99 Bio::Tools::Analysis::SimpleAnalysisBase.
 100
 101 =back
 102
 103 =head1 SEE ALSO
 104
 105 L<Bio::SimpleAnalysisI>,
 106 L<Bio::Tools::Analysis::SimpleAnalysisBase>
 107 L<Bio::Seq::Meta::Array>,
 108 L<Bio::WebAgent>
 109
 110 =head1 FEEDBACK
 111
 112 =head2 Mailing Lists
 113
 114 User feedback is an integral part of the evolution of this and other
 115 Bioperl modules. Send your comments and suggestions preferably to one
 116 of the Bioperl mailing lists.  Your participation is much appreciated.
 117
 118   bioperl-l@bioperl.org                  - General discussion
 119   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 120
 121 =head2 Support
 122
 123 Please direct usage questions or support issues to the mailing list:
 124
 125 I<bioperl-l@bioperl.org>
 126
 127 rather than to the module maintainer directly. Many experienced and
 128 reponsive experts will be able look at the problem and quickly
 129 address it. Please include a thorough description of the problem
 130 with code and data examples if at all possible.
 131
 132 =head2 Reporting Bugs
 133
 134 Report bugs to the Bioperl bug tracking system to help us keep track
 135 the bugs and their resolution.  Bug reports can be submitted via the
 136 web:
 137
 138   https://github.com/bioperl/bioperl-live/issues
 139
 140 =head1 AUTHORS
 141
 142 Richard Adams, Richard.Adams@ed.ac.uk,
 143
 144 =head1 APPENDIX
 145
 146 =cut
 147
 148 use strict;
 149
 150 package Bio::Tools::Analysis::Protein::Sopma;
 151
 152 use IO::String;
 153 use Bio::SeqIO;
 154 use HTTP::Request::Common qw (POST);
 155 use Bio::SeqFeature::Generic;
 156 use Bio::Seq::Meta::Array;
 157 $ENV{PERL_LWP_SSL_VERIFY_HOSTNAME} = 0;
 158
 159 use base qw(Bio::Tools::Analysis::SimpleAnalysisBase);
 160
 161 #extends array for 2struc.
 162 my $URL = 'https://npsa-prabi.ibcp.fr/cgi-bin/secpred_sopma.pl';
 163 my $ANALYSIS_NAME= 'Sopma';
 164 my $ANALYSIS_SPEC= {name => 'Sopma', type => 'Protein'};
 165 my $INPUT_SPEC = [
 166                   {mandatory=>'true',
 167                    type     => 'Bio::PrimarySeqI',
 168                    'name'   => 'seq',
 169                   },
 170                   {mandatory =>'false',
 171                    type      => 'integer',
 172                    name      => 'similarity_threshold',
 173                    default   => 8,
 174                   },
 175                   {mandatory  =>'false',
 176                    type       => 'integer',
 177                    name       => 'window_width',
 178                    default    => 17,
 179                   },
 180                   {mandatory  =>'false',
 181                    type       => 'integer',
 182                    name       => 'states',
 183                    default    => 4,
 184                   },
 185                  ];
 186 my  $RESULT_SPEC =
 187     {
 188      ''   => 'bulk',              # same as undef
 189      raw  => '[{struc=>, helix=>, turn=>, coil=>, sheet=>}]',
 190      meta => 'Bio::Seq::Meta::Array object',
 191      'Bio::SeqFeatureI' => 'ARRAY of Bio::SeqFeature::Generic',
 192     };
 193 use constant MIN_STRUC_LEN => 3;
 194
 195
 196 =head2  similarity_threshold
 197
 198   Useage  : $job->similarity_threshold(...)
 199   Returns : The  similarity threshold used in the analysis
 200   Args    : None (retrieves value) or  an integer (default = 8)
 201             that sets the similarity threshold .
 202
 203 This method gets/sets the  similarity threshold for the prediction.
 204
 205 =cut
 206
 207 sub similarity_threshold {
 208     my ($self, $value) = @_;
 209     if ($value) {
 210         $self->throw ("similarity_threshold must be integer")
 211             unless $value =~ /^\d+$/;
 212         $self->{'_similarity_threshold'} = $value;
 213     }
 214     $self->{'_similarity_threshold'} ||= $self->input_spec->[1]{'default'};
 215     return $self->{'_similarity_threshold'};
 216 }
 217
 218 =head2  window_width
 219
 220   Usage    : $job->window_width(...)
 221   Returns  : The window width used in the analysis
 222   Args     : None (retrieves value) or  an integer (default = 17)
 223              that sets the window width.
 224
 225 This method gets/sets the window width for the prediction, .  If
 226 attempted to set longer than the sequence, warns of error.
 227
 228 =cut
 229
 230 sub window_width {
 231     my ($self, $value) = @_;
 232     if ($value) {
 233         $self->throw ("window_width must be integer")
 234             unless $value =~ /^\d+$/;
 235         $self->{'_window_width'} = $value;
 236     }
 237     $self->{'_window_width'} ||= $self->input_spec->[2]{'default'};
 238     $self->warn ("window width longer than sequence!")
 239         unless $self->{'_window_width'} < $self->seq->length;
 240     return $self->{'_window_width'};
 241 }
 242
 243 =head2  states
 244
 245   Usage    : $job->states(...)
 246   Returns  : The number of secondary structure prediction states
 247   Args     : None (retrieves value) or either '3' or '4' to set
 248              prior to running analysis.
 249
 250 This method gets/sets the number of states for the prediction, either
 251 3 or 4 (includes turns).
 252
 253 =cut
 254
 255 sub states {
 256     my ($self, $value) = @_;
 257     if ($value) {
 258         $self->throw ("number of states must be 3 or 4")
 259             unless $value == 3 or $value ==4;
 260         $self->{'_states'} = $value;
 261     }
 262     $self->{'_states'} ||= $self->input_spec->[3]{'default'};
 263     return $self->{'_states'};
 264 }
 265
 266 =head2 result
 267
 268   Usage   : $job->result (...)
 269   Returns : a result created by running an analysis
 270   Args    : various
 271
 272 The method returns a result of an executed job. If the job was
 273 terminated by an error the result may contain an error message instead
 274 of the real data.
 275
 276 This implementation returns differently processed data depending on
 277 argument:
 278
 279 =over 3
 280
 281 =item undef
 282
 283 Returns the raw ASCII data stream but without HTML tags
 284
 285 =item 'Bio::SeqFeatureI'
 286
 287 The argument string defines the type of bioperl objects returned in an
 288 array.  The objects are L<Bio::SeqFeature::Generic>.  Feature primary
 289 tag is "2ary".  Feature tags are "type" (which can be helix, sheet
 290 coil, or turn if 4 state prediction requested) "method" (Sopma)
 291
 292 =item 'parsed'
 293
 294 Array of hash references of scores/structure assignations
 295 { helix =E<gt> , sheet =E<gt> , coil =E<gt> , struc=E<gt>}.
 296
 297 =item 'all'
 298
 299 A Bio::Seq::Meta::Array object. Scores can be accessed using methods
 300 from this class. Meta sequence names are Sopma_helix, Sopma_sheet,
 301 Sopma_coil, Sopma_turn (if defined), and Sopma_struc.
 302
 303
 304 =back
 305
 306
 307 =cut
 308
 309 sub result {
 310     my ($self,$value, $run_id) = @_;
 311
 312     my @score;
 313     my @fts;
 314
 315     if ($value ) {
 316         if (!exists($self->{'_parsed'} )) {
 317             my $result = IO::String->new($self->{'_result'});
 318             while (my $line = <$result>) {
 319                 next unless $line =~ /^[HCET]\s/; # or for sopma/hnn  /^[A-Z]\s/
 320                 $line =~/^([A-Z])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/; # or for so
 321                 push @score, { struc => $1,
 322                                helix => $2,
 323                                sheet => $3,
 324                                coil => $5,
 325                              };
 326                 #include turn if 4states are requested
 327                 $score[$#score]{'turn'} = $4 if $self->states == 4;
 328                 #can optimize by duplicating code here
 329             }
 330             $self->{'_parsed'} = \@score;
 331         }
 332         if ($value eq 'Bio::SeqFeatureI') {
 333             $self->_get_2ary_coords();
 334             for my $type (keys %{$self->{'_parsed_coords'}} ) {
 335                 next if $type =~  /\w{2,}/; #if not H,C,E or T
 336
 337                                 ## these 2 are added to distinguish features on same
 338                ## sequence run with different params
 339                                 my $tag_hash = {
 340                                                                 type   => $type,
 341                                 method => $self->analysis_name,
 342                                                                 };
 343                                 $self->_add_params_to_result($tag_hash);
 344
 345                                 ## now make feature object
 346                 for my $loc (@{$self->{'_parsed_coords'}{$type}} ) {
 347                     push  @fts,   Bio::SeqFeature::Generic->new
 348                         (-start   => $loc->{'start'},
 349                          -end     => $loc->{'end'},
 350                          -source  => 'Sopma',
 351                          -primary => 'Domain',
 352                          -tag => $tag_hash,
 353                                  );
 354                 }               #end of array of strucs of type
 355             }                   # end of all 2nd struc elements
 356             delete $self->{'_parsed_coords'}; #remove temp data
 357             return @fts;
 358         }                       #endif BioSeqFeature
 359
 360         elsif ($value eq 'meta') {
 361             #1st of all make 3 or 4 arrays of scores for each type from column data
 362             my %type_scores;
 363             for my $aa (@{$self->{'_parsed'}}) {
 364                 for my $type (qw(struc helix sheet  coil)) {
 365                     push @{$type_scores{$type}}, $aa->{$type};
 366                 }
 367                 push @{$type_scores{'turn'}}, $aa->{'turn'} if  exists $aa->{'turn'};
 368             }
 369
 370                         ## convert to meta sequence array ##
 371                         if (!$self->seq->isa("Bio::Seq::Meta::Array")) {
 372                          bless ($self->seq, "Bio::Seq::Meta::Array");
 373                                 }
 374             $self->seq->isa("Bio::Seq::MetaI")
 375                 || $self->throw("$self is not a Bio::Seq::MetaI");
 376
 377
 378             $Bio::Seq::Meta::Array::DEFAULT_NAME = 'Sopma_struc';
 379             for my $struc_type (keys %type_scores) {
 380                 my $meta_name = "Sopma". "_" . "$struc_type";
 381                                 if ($run_id) {
 382                                         $meta_name .= "|$run_id";
 383                                 }
 384                 my @meta = map{$_->{$struc_type}} @{$self->{'_parsed'}};
 385                 if (grep{$_ eq $meta_name}$self->seq->meta_names >0) {
 386                     $self->warn ("$meta_name already exists , not overwriting!");
 387                     next;
 388                 }
 389                 $self->seq->named_meta($meta_name,\@meta );
 390             }
 391             # return  seq array object implementing meta sequence #
 392             return $self->seq;
 393
 394         }
 395                 ## else return parsed data if $value is defined
 396                  else {
 397             return $self->{'_parsed'};
 398         }
 399
 400     }                           #endif ($value)
 401     #return raw result if no return format stated
 402     return $self->{'_result'};
 403 }
 404
 405 sub _init {
 406     my $self = shift;
 407     $self->url($URL);
 408     $self->{'_ANALYSIS_SPEC'} = $ANALYSIS_SPEC;
 409     $self->{'_INPUT_SPEC'}    = $INPUT_SPEC;
 410     $self->{'_RESULT_SPEC'}   = $RESULT_SPEC;
 411     $self->{'_ANALYSIS_NAME'} = $ANALYSIS_NAME;
 412     return $self;
 413 }
 414
 415 sub _get_2ary_coords {
 416     #helper sub for result;
 417     ##extracts runs of structure > MIN_STRUC_LENresidues or less if Turn:
 418     #i.e., helical prediction for 1 residue isn't very meaningful...
 419     ## and poulates array of hashes with start/end values.
 420     ##keys of $Result are 'H' 'T' 'C' 'E'.
 421     my ($self) = @_;
 422     my @prot = @{$self->{'_parsed'}};
 423     my %Result;
 424     for (my $index = 0; $index <= $#prot; $index++) {
 425
 426         my $type        = $prot[$index]{'struc'};
 427         next unless $type && $type =~ /[HTCE]/;
 428         my $length = 1;
 429         for (my $j = $index + 1; $j <= $#prot; $j++) {
 430             my $test = $prot[$j];
 431             if ($test->{'struc'} eq $type) {
 432                 $length++;
 433             } elsif (  $length > MIN_STRUC_LEN  ||
 434                        ($length <= MIN_STRUC_LEN && $type eq 'T') ) {
 435                 push @{$Result{$type}}, {start => $index + 1 ,  end => $j};
 436                 $index += $length -1;
 437                 last;
 438             } else {
 439                 $index += $length - 1;
 440                 last;
 441             }
 442         }
 443     }
 444     $self->{'_parsed_coords'} = \%Result; #temp assignment
 445 }
 446
 447 sub  _run {
 448     my $self  = shift;
 449     $self->delay(1);
 450     # delay repeated calls by default by 3 sec, set delay() to change
 451     $self->sleep;
 452     $self->status('TERMINATED_BY_ERROR');
 453     my $request = POST 'https://npsa-prabi.ibcp.fr/cgi-bin/secpred_sopma.pl',
 454         Content_Type => 'form-data',
 455             Content  => [title     => "",
 456                          notice    => $self->seq->seq,
 457                          ali_width => 70,
 458                          states    => $self->states,
 459                          threshold => $self->similarity_threshold ,
 460                          width     => $self->window_width,
 461                         ];
 462
 463     my $text = $self->request($request)->content;
 464     return $self unless $text;
 465
 466     #### get text only version of results ##
 467     my ($next) = $text =~ /Prediction.*?=(.*?)>/;
 468     return $self unless $next;
 469     my $out    = "http://npsa-pbil.ibcp.fr/". "$next";
 470     my $req2   = HTTP::Request->new(GET=>$out);
 471     my $resp2  = $self->request ($req2);
 472     $self->{'_result'} = $resp2->content;
 473     $self->status('COMPLETED') if $resp2 ne '';
 474     return $self;
 475 }
 476
 477 sub _add_params_to_result{
 478         ## called when making Seqfeature objects
 479         my ($self, $tag_hash) = @_;
 480         my $hash;
 481         ## adds input parameter values to SeqFeatureI results where multiple
 482     ##  parameter values are possible. Only adds value if not default.
 483         map{$hash->{$_->{'name'}} = $_}@{$self->input_spec()};
 484
 485         for my $p (keys %$hash) {
 486                 if (!ref($self->$p) && $self->$p ne $hash->{$p}{'default'}) {
 487                         $tag_hash->{$p} = $self->$p;
 488                 }
 489         }
 490
 491 }
 492
 493
 494 1;