Bio::Tools::Signalp::ExtendedSignalp - enhanced parser for Signalp output

 use Bio::Tools::Signalp::ExtendedSignalp;
 my $params = [qw(maxC maxY maxS meanS D)];
 my $parser = new Bio::Tools::Signalp::ExtendedSignalp(
                                                       -fh      => $filehandle
                                                       -factors => $params
                                                      );

 $parser->factors($params);
 while( my $sp_feat = $parser->next_feature ) {
       #do something
       #eg
       push @sp_feat, $sp_feat;
 }

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  https://redmine.open-bio.org/projects/bioperl/

 Based on the Bio::Tools::Signalp module
 Emmanuel Quevillon <emmanuel.quevillon@versailles.inra.fr>

 The rest of the documentation details each of the object methods.
 Internal methods are usually preceded with a _

 Title   : new
 Usage   : my $obj = new Bio::Tools::Signalp::ExtendedSignalp();
 Function: Builds a new Bio::Tools::Signalp::ExtendedSignalp object
 Returns : Bio::Tools::Signalp::ExtendedSignalp
 Args    : -fh/-file => $val, # for initing input, see Bio::Root::IO

 Title   : next_feature
 Usage   : my $feat = $signalp->next_feature
 Function: Get the next result feature from parser data
 Returns : Bio::SeqFeature::Generic
 Args    : none

 Title   : _filterok
 Usage   : my $feat = $signalp->_filterok
 Function: Check if the factors required by the user are all ok.
 Returns : 1/0
 Args    : hash reference

 Title   : factors
 Usage   : my $feat = $signalp->factors
 Function: Get/Set the filters required from the user
 Returns : hash
 Args    : array reference

 Title   : _parsed
 Usage   : obj->_parsed()
 Function: Get/Set if the result is parsed or not
 Returns : 1/0 scalar
 Args    : On set 1

 Title   : _parse
 Usage   : obj->_parse
 Function: Parse the SignalP result
 Returns :
 Args    :

 Title   : _parse_summary_format
 Usage   : $self->_parse_summary_format
 Function: Method to parse summary/full format from signalp output
           It automatically fills filtered features.
 Returns :
 Args    :

 Title   : _parse_nn_result
 Usage   : obj->_parse_nn_result
 Function: Parses the Neuronal Network (NN) part of the result
 Returns : Hash reference
 Args    :

 Title   : _parse_hmm_result
 Usage   : obj->_parse_hmm_result
 Function: Parses the Hiden Markov Model (HMM) part of the result
 Returns : Hash reference
 Args    :

 Title   : _parse_short_format
 Usage   : $self->_parse_short_format
 Function: Method to parse short format from signalp output
           It automatically fills filtered features.
 Returns :
 Args    :

 Title   : create_feature
 Usage   : obj->create_feature(\%feature)
 Function: Internal(not to be used directly)
 Returns :
 Args    :

 Title   : seqname
 Usage   : obj->seqname($name)
 Function: Internal(not to be used directly)
 Returns :
 Args    :

#
# BioPerl module for Bio::Tools::Signalp::ExtendedSignalp
#
# Please direct questions and support issues to <bioperl-l@bioperl.org> 
#
# Cared for by Emmanuel Quevillon <emmanuel.quevillon@versailles.inra.fr>
#
# Copyright Emmanuel Quevillon
#
# You may distribute this module under the same terms as perl itself
#
# POD documentation - main docs before the code

package Bio::Tools::Signalp::ExtendedSignalp;

use strict;
use Data::Dumper;
use Bio::SeqFeature::Generic;
# don't need Bio::Root::Root/IO (already in inheritance tree)
use base qw(Bio::Tools::Signalp Bio::Tools::AnalysisResult);

#Supported arguments
my $FACTS = {
	     'maxC'  => 1,
	     'maxS'  => 1,
	     'maxY'  => 1,
	     'meanS' => 1,
	     'D'     => 1,
	    };

sub new {
      my($class,@args) = @_;

      my $self = $class->SUPER::new(@args);
      $self->_initialize_io(@args);

      my $factors = $self->_rearrange([qw(FACTORS)], @args);
      #To behave like the parent module (Bio::Tools::Signalp) we default factors to these two factors
      if($factors && scalar(@$factors)){
	  $factors = $factors;
      }
      else{
	  $factors = [qw(maxY meanS)];
      }
      $factors && $self->factors($factors);
						
      return $self;
}

sub next_feature {

    my ($self) = @_;

    if(!$self->_parsed()){
	$self->_parse();
    }

    return shift @{$self->{_features}} || undef;

}

sub _filterok {

    my($self, $hash) = @_;

    #We hope everything will be fine ;)
    my $bool = 1;

    #If the user did not give any filter, we keep eveything
    return $bool unless keys %{$self->{_factors}};

    #If only one of the factors parsed is equal to NO based on the user factors cutoff
    #Then the filter is not ok.
    foreach my $fact (keys %{$self->factors()}){
	if(exists($hash->{$fact}) && $hash->{$fact} =~ /^N/){
	    $bool = 0;
	}
    }

    return $bool;

}

sub factors {

    my($self, $array) = @_;

    if($array){
	$self->{_factors} = { };
	foreach my $f (@$array){
	    if(exists($FACTS->{$f})){
		$self->{_factors}->{$f} = 1;
	    }
	    else{
		$self->throw("[$f] incorrect factor. Supported:\n- ".join("\n- ", keys %$FACTS)."\n");
	    }
	}
    }

    return $self->{_factors};

}

sub _parsed {

    my($self, $parsed) = @_;

    if(defined($parsed)){
	$self->{_parsed} = $parsed;
    }

    return $self->{_parsed};

}

sub _parse {

    my($self) = @_;

    #Let's read the file...
    while (my $line = $self->_readline()) {

	chomp $line;
	#We want to be sure to catch the first non empty line to be ablte to determine
	#which format we are working with...
	next unless ($line =~ /^>(\S+)|^# SignalP-[NHM]+ \S+ predictions/);

	if($line =~ /^>(\S+)/){
	    $self->_pushback($line);
	    $self->_parse_summary_format();
	    last;
	}
	elsif($line =~ /^# SignalP-[NHM]+ \S+ predictions/){
	    $self->_pushback($line);
	    $self->_parse_short_format();
	    last;
	}
	else{
	    $self->throw("Unable to determine the format type.");
	}
    }

    return;
}

sub _parse_summary_format {

    my($self) = @_;

    my $feature = undef;
    my $ok = 0;

    while(my $line = $self->_readline()){

	if($line =~ /^SignalP-NN result:/){
	    $self->_pushback($line);
	    $feature = $self->_parse_nn_result($feature);
	}
	if($line =~ /^SignalP-HMM result:/){
	    $self->_pushback($line);
	    $feature = $self->_parse_hmm_result($feature);
	}

	if($line =~ /^---------/ && $feature){
	    my $new_feature = $self->create_feature($feature);
	    push @{$self->{_features}}, $new_feature if $new_feature;
	    $feature = undef;
	}
    }

    return;
}

sub _parse_nn_result {

    my($self, $feature) = @_;

    my $ok   = 0;
    my %facts;

    #SignalP-NN result:
    #>MGG_11635.5           length = 100
    ## Measure  Position  Value  Cutoff  signal peptide?
    #  max. C    37       0.087   0.32   NO
    #  max. Y    37       0.042   0.33   NO
    #  max. S     3       0.062   0.87   NO
    #  mean S     1-36    0.024   0.48   NO
    #       D     1-36    0.033   0.43   NO

    while(my $line = $self->_readline()){

	chomp $line;

	if($line =~ /^SignalP-NN result:/){
	    $ok = 1;
	    next;
	}

	$self->throw("Wrong line for parsing NN results.") unless $ok;

	if ($line=~/^\>(\S+)\s+length/) {
	    $self->seqname($1);
	    %facts = ();
	    next;
	}
	elsif($line =~ /max\.\s+C\s+(\S+)\s+\S+\s+\S+\s+(\S+)/) {
	    $feature->{maxCprob} = $1;
	    $facts{maxC} = $2;
	    next;
	}
	elsif ($line =~ /max\.\s+Y\s+(\S+)\s+\S+\s+\S+\s+(\S+)/) {
	    $feature->{maxYprob} = $1;
	    $facts{maxY} = $2;
	    next;
	}
	elsif($line =~ /max\.\s+S\s+(\S+)\s+\S+\s+\S+\s+(\S+)/) {
	    $feature->{maxSprob} = $1;
	    $facts{maxS} = $2;
	    next;
	}
	elsif ($line=~/mean\s+S\s+(\S+)\s+\S+\s+\S+\s+(\S+)/) {
	    $feature->{meanSprob} = $1;
	    $facts{meanS} = $2;
	    next;
	}
	elsif ($line=~/\s+D\s+(\S+)\s+\S+\s+\S+\s+(\S+)/) {
	    $feature->{Dprob} = $1;
	    $facts{D} = $2;
	    next;
	}
	#If we don't have this line it means that all the factors cutoff are equal to 'NO'
	elsif ($line =~ /Most likely cleavage site between pos\.\s+(\d+)/) {
	    #if($self->_filterok(\%facts)){
		#$feature->{name}       = $self->seqname();
		#$feature->{start}      = 1;
		$feature->{end}        = $1 + 1; #To be consistent with end given in short format
	    #}
	    #return $feature;
	}
	elsif($line =~ /^\s*$/){
	    last;
	}
    }

    if($self->_filterok(\%facts)){
	$feature->{name}  = $self->seqname();
	$feature->{start} = 1;
	$feature->{nnPrediction} = 'signal-peptide';
    }

    return $feature;
}

sub _parse_hmm_result {

    my ($self, $feature_hash) = @_;

    my $ok = 0;

    #SignalP-HMM result:
    #>MGG_11635.5
    #Prediction: Non-secretory protein
    #Signal peptide probability: 0.000
    #Signal anchor probability: 0.000
    #Max cleavage site probability: 0.000 between pos. -1 and  0

    while(my $line = $self->_readline()){

	chomp $line;
	next if $line =~ /^\s*$/o;

	if($line =~ /^SignalP-HMM result:/){
	    $ok = 1;
	    next;
	}

	$self->throw("Wrong line for parsing HMM result.") unless $ok;

	if($line =~ /^>(\S+)/){
	    #In case we already seen a name with NN results
	    $feature_hash->{name} = $1 unless $self->seqname();
	}
        elsif($line =~ /Prediction: (.+)$/){
            $feature_hash->{hmmPrediction} = $1;
        }
	elsif($line =~ /Signal peptide probability: ([0-9\.]+)/){
            $feature_hash->{peptideProb} = $1;
        }
	elsif($line =~ /Signal anchor probability: ([0-9\.]+)/){
            $feature_hash->{anchorProb} = $1;
        }
	elsif($line =~ /Max cleavage site probability: (\S+) between pos. \S+ and (\S+)/){
	    $feature_hash->{cleavageSiteProb} = $1;
	    #Strange case, if we don't have an end value in NN result (no nn method launched)
	    #We try anyway to get an end value, unless this value is lower than 1 which is
	    #the start
	    $feature_hash->{end}   = $2 if($2 > 1 && !$feature_hash->{end});
	    $feature_hash->{start} = 1 unless $feature_hash->{start};
	    last;
	}
    }

    return $feature_hash;
}

sub _parse_short_format {

                                my($self) = @_;

    my $ok = 0;
    my $method = undef;
    $self->{_oformat} = 'short';

    #Output example
    # SignalP-NN euk predictions                                   	                # SignalP-HMM euk predictions
    # name                Cmax  pos ?  Ymax  pos ?  Smax  pos ?  Smean ?  D     ? 	# name      !  Cmax  pos ?  Sprob ?
    #Q5A8M1_CANAL          0.085  27 N  0.190  35 N  0.936  27 Y  0.418 N  0.304 N	Q5A8M1_CANAL  Q  0.001  35 N  0.002 N
    #O74127_YARLI          0.121  21 N  0.284  21 N  0.953  11 Y  0.826 Y  0.555 Y	O74127_YARLI  S  0.485  23 N  0.668 Y
    #Q5VJ86_9PEZI          0.355  24 Y  0.375  24 Y  0.798  12 N  0.447 N  0.411 N	Q5VJ86_9PEZI  Q  0.180  23 N  0.339 N
    #Q5A8U5_CANAL          0.085  27 N  0.190  35 N  0.936  27 Y  0.418 N  0.304 N	Q5A8U5_CANAL  Q  0.001  35 N  0.002 N

    while(my $line = $self->_readline()){
	
	chomp $line;
	next if $line =~ /^\s*$|^# name/;

	if($line =~ /^#/){
	    $method = $line =~ /SignalP-NN .+ SignalP-HMM/ ?
	                                            'both' : $line =~ /SignalP-NN/ ?
							                      'nn' : 'hmm';
	    next;
	}

	#$self->throw("It looks like the format is not 'short' format.") unless($ok);

	my @data = split(/\s+/, $line);
	$self->seqname($data[0]);

	my $factors = { };
	my $feature = { };

	#NN results gives more fields than HMM
	if($method eq 'both' || $method eq 'nn'){

	    $feature->{maxCprob} = $data[1];
	    $factors->{maxC}     = $data[3];
	    $feature->{maxYprob} = $data[4];
	    $factors->{maxY}     = $data[6];
	    $feature->{maxSprob} = $data[7];
	    $factors->{maxS}     = $data[9];
	    $feature->{meanSprob}= $data[10];
	    $factors->{meanS}    = $data[11];
	    $feature->{Dprob}    = $data[12];
	    $factors->{D}        = $data[13];
	    #It looks like the max Y position is reported as the most likely cleavage position
	    $feature->{end}      = $data[5];
	    $feature->{nnPrediction} = 'signal-peptide';

	    if($method eq 'both'){
		$feature->{hmmPrediction}    = $data[15] eq 'Q' ? 'Non-secretory protein' : 'Signal peptide';
		$feature->{cleavageSiteProb} = $data[16];
		$feature->{peptideProb}      = $data[19];
	    }
	}
	elsif($method eq 'hmm'){
	    #In short output anchor probability is not given
	    $feature->{hmmPrediction}    = $data[1] eq 'Q' ? 'Non-secretory protein' : 'Signal peptide';
	    $feature->{cleavageSiteProb} = $data[2];
	    $feature->{peptideProb}      = $data[5];
	    #It looks like the max cleavage probability position is given by the Cmax proability
	    $feature->{end} = $data[3];
	}

	#Unfortunately, we cannot parse the filters for hmm method.
	if($self->_filterok($factors)){
	    $feature->{name}        = $self->seqname();
	    $feature->{start}       = 1;
	    $feature->{source}      = 'Signalp';
	    $feature->{primary}     = 'signal_peptide';
	    $feature->{program}     = 'Signalp';
	    $feature->{logic_name}  = 'signal_peptide';

	    my $new_feat = $self->create_feature($feature);
	    push @{$self->{_features}}, $new_feat if $new_feat;
	}
    }

    return;
}

sub create_feature {

    my ($self, $feat) = @_;

    #If we don't have neither start nor end, we return.
    unless($feat->{name} && $feat->{start} && $feat->{end}){
	return;
    }

    # create feature object
    my $feature = Bio::SeqFeature::Generic->new(
						-seq_id     => $feat->{name},
						-start      => $feat->{start},
						-end        => $feat->{end},
						-score      => defined($feat->{peptideProb}) ? $feat->{peptideProb} : '',
						-source     => 'Signalp',
						-primary    => 'signal_peptide',
						-logic_name => 'signal_peptide',
					       );

    $feature->add_tag_value('peptideProb', $feat->{peptideProb});
    $feature->add_tag_value('anchorProb', $feat->{anchorProb});
    $feature->add_tag_value('evalue',$feat->{anchorProb});
    $feature->add_tag_value('percent_id','NULL');
    $feature->add_tag_value("hid",$feat->{primary});
    $feature->add_tag_value('signalpPrediction', $feat->{hmmPrediction});
    $feature->add_tag_value('cleavageSiteProb', $feat->{cleavageSiteProb}) if($feat->{cleavageSiteProb});
    $feature->add_tag_value('nnPrediction', $feat->{nnPrediction})         if($feat->{nnPrediction});
    $feature->add_tag_value('maxCprob', $feat->{maxCprob})   if(defined($feat->{maxCprob}));
    $feature->add_tag_value('maxSprob', $feat->{maxSprob})   if(defined($feat->{maxSprob}));
    $feature->add_tag_value('maxYprob', $feat->{maxYprob})   if(defined($feat->{maxYprob}));
    $feature->add_tag_value('meanSprob', $feat->{meanSprob}) if(defined($feat->{meanSprob}));
    $feature->add_tag_value('Dprob', $feat->{Dprob})         if(defined($feat->{Dprob}));

    return $feature;

}

sub seqname{
    my ($self,$seqname)=@_;

    if (defined($seqname)){
        $self->{'seqname'} = $seqname;
    }

    return $self->{'seqname'};

}


1;