t/AlignIO/AlignIO.t: fix number of tests in plan (fixup c523e6bed866)
[bioperl-live.git] / Bio / Tools / Analysis / Protein / Sopma.pm
blob526554b66f1d4061f0d4eb00c729086018993cea
1 # $Id: Sopma.pm,v 1.0 2003/07/ 11
3 # BioPerl module for Bio::Tools::Analysis::Protein::Sopma
5 # Copyright Richard Adams
7 # You may distribute this module under the same terms as perl itself
9 # POD documentation - main docs before the code
12 =head1 NAME
14 Bio::Tools::Analysis::Protein::Sopma - a wrapper around the
15 Sopma protein secondary structure prediction server
17 =head1 SYNOPSIS
19 use Bio::Tools::Analysis::Protein::Sopma;
20 #get a Bio::Seq or Bio::PrimarySeq
21 my $seq;
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
28 =head1 DESCRIPTION
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.,
35 my $analysis_object = Bio::Tools::SimpleAnalysis::Protein::Sopma->new
36 ( -seq => $seq,
37 -states => 4,
38 -window_width => 15,
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.,
45 $analysis_object->states(4);
46 $analysis_object->run;
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:
52 =over 4
54 =item 1
56 The raw text of the program output.
58 my $rawdata = $analysis_object->result;
60 =item 2
62 A reference to an array of hashes of scores for each state and the
63 assigned state.
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";
69 Hash keys are 'helix', 'struc', 'sheet', 'coil', 'turn'.
71 =item 3
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.
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";
83 =item 4
85 A Bio::Seq::Meta::Array implementing sequence.
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.,
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";
95 Meta sequence names are : Sopma_helix, Sopma_sheet, Sopma_turn,
96 Sopma_coil, Sopma_struc, representing the scores for each residue.
98 Many methods common to all analyses are inherited from
99 Bio::Tools::Analysis::SimpleAnalysisBase.
101 =back
103 =head1 SEE ALSO
105 L<Bio::SimpleAnalysisI>,
106 L<Bio::Tools::Analysis::SimpleAnalysisBase>
107 L<Bio::Seq::Meta::Array>,
108 L<Bio::WebAgent>
110 =head1 FEEDBACK
112 =head2 Mailing Lists
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.
118 bioperl-l@bioperl.org - General discussion
119 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
121 =head2 Support
123 Please direct usage questions or support issues to the mailing list:
125 I<bioperl-l@bioperl.org>
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.
132 =head2 Reporting Bugs
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:
138 https://github.com/bioperl/bioperl-live/issues
140 =head1 AUTHORS
142 Richard Adams, Richard.Adams@ed.ac.uk,
144 =head1 APPENDIX
146 =cut
148 use strict;
150 package Bio::Tools::Analysis::Protein::Sopma;
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;
159 use base qw(Bio::Tools::Analysis::SimpleAnalysisBase);
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',
170 {mandatory =>'false',
171 type => 'integer',
172 name => 'similarity_threshold',
173 default => 8,
175 {mandatory =>'false',
176 type => 'integer',
177 name => 'window_width',
178 default => 17,
180 {mandatory =>'false',
181 type => 'integer',
182 name => 'states',
183 default => 4,
186 my $RESULT_SPEC =
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',
193 use constant MIN_STRUC_LEN => 3;
196 =head2 similarity_threshold
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 .
203 This method gets/sets the similarity threshold for the prediction.
205 =cut
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;
214 $self->{'_similarity_threshold'} ||= $self->input_spec->[1]{'default'};
215 return $self->{'_similarity_threshold'};
218 =head2 window_width
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.
225 This method gets/sets the window width for the prediction, . If
226 attempted to set longer than the sequence, warns of error.
228 =cut
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;
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'};
243 =head2 states
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.
250 This method gets/sets the number of states for the prediction, either
251 3 or 4 (includes turns).
253 =cut
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;
262 $self->{'_states'} ||= $self->input_spec->[3]{'default'};
263 return $self->{'_states'};
266 =head2 result
268 Usage : $job->result (...)
269 Returns : a result created by running an analysis
270 Args : various
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.
276 This implementation returns differently processed data depending on
277 argument:
279 =over 3
281 =item undef
283 Returns the raw ASCII data stream but without HTML tags
285 =item 'Bio::SeqFeatureI'
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)
292 =item 'parsed'
294 Array of hash references of scores/structure assignations
295 { helix =E<gt> , sheet =E<gt> , coil =E<gt> , struc=E<gt>}.
297 =item 'all'
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.
304 =back
307 =cut
309 sub result {
310 my ($self,$value, $run_id) = @_;
312 my @score;
313 my @fts;
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,
326 #include turn if 4states are requested
327 $score[$#score]{'turn'} = $4 if $self->states == 4;
328 #can optimize by duplicating code here
330 $self->{'_parsed'} = \@score;
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
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,
343 $self->_add_params_to_result($tag_hash);
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,
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
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};
367 push @{$type_scores{'turn'}}, $aa->{'turn'} if exists $aa->{'turn'};
370 ## convert to meta sequence array ##
371 if (!$self->seq->isa("Bio::Seq::Meta::Array")) {
372 bless ($self->seq, "Bio::Seq::Meta::Array");
374 $self->seq->isa("Bio::Seq::MetaI")
375 || $self->throw("$self is not a Bio::Seq::MetaI");
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";
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;
389 $self->seq->named_meta($meta_name,\@meta );
391 # return seq array object implementing meta sequence #
392 return $self->seq;
395 ## else return parsed data if $value is defined
396 else {
397 return $self->{'_parsed'};
400 } #endif ($value)
401 #return raw result if no return format stated
402 return $self->{'_result'};
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;
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++) {
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;
444 $self->{'_parsed_coords'} = \%Result; #temp assignment
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,
463 my $text = $self->request($request)->content;
464 return $self unless $text;
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;
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()};
485 for my $p (keys %$hash) {
486 if (!ref($self->$p) && $self->$p ne $hash->{$p}{'default'}) {
487 $tag_hash->{$p} = $self->$p;