Bio::DB::Universal: move into its own distribution

[bioperl-live.git] / Bio / Tools / Signalp / ExtendedSignalp.pm

#
# 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

=head1 NAME

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

=head1 SYNOPSIS

 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;
 }

=head1 DESCRIPTION

# Please direct questions and support issues to I<bioperl-l@bioperl.org> 

Parser module for Signalp.

Based on the EnsEMBL module Bio::EnsEMBL::Pipeline::Runnable::Protein::Signalp
originally written by Marc Sohrmann (ms2 a sanger.ac.uk) Written in BioPipe by
Balamurugan Kumarasamy (savikalpa a fugu-sg.org) Cared for by the Fugu
Informatics team (fuguteam@fugu-sg.org)

You may distribute this module under the same terms as perl itself

Compared to the original SignalP, this method allow the user to filter results
out based on maxC maxY maxS meanS and D factor cutoff for the Neural Network (NN)
method only. The HMM method does not give any filters with 'YES' or 'NO' as result.

The user must be aware that the filters can only by applied on NN method.
Also, to ensure the compatibility with original Signalp parsing module, the user
must know that by default, if filters are empty, max Y and mean S filters are
automatically used to filter results.

If the used gives a list, then the parser will only report protein having 'YES'
for each factor.

This module supports parsing for full, summary and short output form signalp.
Actually, full and summary are equivalent in terms of filtering results.

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list.  Your participation is much appreciated.

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

=head2 Support 

Please direct usage questions or support issues to the mailing list:

I<bioperl-l@bioperl.org>

rather than to the module maintainer directly. Many experienced and 
reponsive experts will be able look at the problem and quickly 
address it. Please include a thorough description of the problem 
with code and data examples if at all possible.

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via
the web:

  https://github.com/bioperl/bioperl-live/issues

=head1 AUTHOR

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

=head1 APPENDIX

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

=cut

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,
            };

=head2 new

 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


=cut

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;
}

=head2 next_feature

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


=cut

sub next_feature {

    my ($self) = @_;

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

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

}

=head2 _filterok

 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


=cut

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;

}

=head2 factors

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


=cut

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};

}

=head2 _parsed

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


=cut

sub _parsed {

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

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

    return $self->{_parsed};

}

=head2 _parse

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


=cut

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;
}

=head2 _parse_summary_format

 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    :

=cut

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;
}


=head2 _parse_nn_result

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


=cut

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;
}


=head2 _parse_hmm_result

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

=cut

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;
}

=head2 _parse_short_format

 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    :

=cut

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;
}

=head2 create_feature

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


=cut

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;

}

=head2 seqname

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


=cut

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

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

    return $self->{'seqname'};

}


1;