#!/usr/local/bin/perl5.6.0 -w

=head1 NAME

signalp2embl - Generate EMBL annotation of signal sequences predicted
by SignalP2.0

=head1 SYNOPSIS

Examples:

 signalp2embl -t 11 -y gram- embl_file > signalp.tab
 signalp2embl -t 1 -y euk -c 0 embl_file > signalp.tab

The EMBL feature table is written to STDOUT.

=head1 DESCRIPTION

This program requires an EMBL file containing annotated CDS
features. The features will be translated using the specified
translation table (default is table 11; bacterial), except for those
with /pseudo, /partial or /transl_except qualifiers.

The output is parsed by the script to produce an EMBL table on
STDOUT.

Setting the -y switch (to euk, gram+ or gram- is important, as is
choosing the correct translation table with the -t switch. The -c
switch sets the colour of the features (using the Artemis colour
codes). Set the -p switch to the threshold probability (0-1.0) below
which you do not want to report predictions.

This script requires the Bio::PSU Perl modules see:

http://www.sanger.ac.uk/Users/kdj/software.html

Don't forget to cite the authors if you use this:

"Identification of prokaryotic and eukaryotic signal peptides and
prediction of their cleavage sites", Henrik Nielsen, Jacob
Engelbrecht, Soeren Brunak and Gunnar von Heijne Protein Engineering,
vol. 10, pages 1-6, 1997.

"Prediction of signal peptides and signal anchors by a hidden Markov
model", Henrik Nielsen and Anders Krogh In: Proceedings of the Sixth
International Conference on Intelligent Sys- tems for Molecular
Biology (ISMB 6), AAAI Press, Menlo Park, California, pages 122-130,
1998.

=head1 METHODS

None

=head1 AUTHOR

Keith James (kdj@sanger.ac.uk)

=head1 COPYRIGHT

Copyright (C) 2001 Genome Research Limited. All Rights Reserved.

decodeanhmm, version 1.1e Copyright (C) 1998-1999 by Anders Krogh. All
rights reserved.

=head1 DISCLAIMER

This software is provided "as is" without warranty of any kind. It may
be used, redistributed and/or modified under the same conditions as
Perl itself.

=cut

use strict;
use Getopt::Std;
use Bio::PSU::SeqFactory;
use Bio::PSU::FeatureFactory;

# Configuration of SignalP variables
my $bindir       = "/nfs/disk222/yeastpub/bio-soft/signalp-2.0";
my $signalp_bin  = "$bindir/signalp";
my $signalp_args = "-trunc 70 -f summary -m hmm";

my (%opts, $colour, $prob, $table, $type);
my $usage = "Usage: signalp2embl [-c <colour>] [-p prob] [-t <table>] -y <type> file <signalp options>

\t-c Colour code of features  (default  0)
\t-p Cutoff probability       (default 0.6)
\t-t Translation table number (default 11)
\t-y Organism type            (gram+, gram-, euk)

\tTables available are:

\t 1) Standard
\t 2) Vertebrate Mitochondrial
\t 3) Yeast Mitochondrial
\t11) Bacterial
";

# Process options
getopts('c:hp:t:y:', \%opts);

defined $opts{h} and die $usage;
defined $opts{c} and $colour  = $opts{c};
defined $opts{p} and $prob    = $opts{p};
defined $opts{t} and $table   = $opts{t};
defined $opts{y} and $type    = $opts{y};

# Default colour is white
$colour = 0 unless defined $colour;

# Default cutoff is 0.6
$prob = 0.6 unless defined $prob;

# Default translation table is bacterial (11)
$table ||= '11';

unless ($type =~ /^(gram\+$|gram-|euk)/) { die $usage }
$signalp_args .= " -t $type";

my $file = shift;
unless (defined $file) { die $usage }

# Read EMBL stream and get entry
my $es = Bio::PSU::SeqFactory->make(-file => "<$file");
my $e  = $es->next_seq;

# Handler for printing features as they appear, rather than printing a
# big feature table at the end
my $fth = Bio::PSU::IO::FTHandler->new;

my $anf = Bio::PSU::FeatureFactory->new;

foreach my $f ($e->features)
{
    # Discard non-CDS features
    next unless $f->key eq 'CDS';

    # Discard any CDS with /pseudo or /partial qualifier
    next if $f->qexists('pseudo|partial|transl_except');

	# Find and fix the id
	my $id = $f->gene;

	$id ||= "Unknown_CDS";
	$id   =~ s/[,;:]$//;

    # Translate with specified translation table
    my $aa = $f->translate($table);

    # Fasta output stream
    my $time        = time;
    my $signalp_tmp = "signalp.$time";

    open(FA, ">$signalp_tmp");
    my $fo = Bio::PSU::SeqFactory->make(-fh     => \*FA,
										-format => 'fasta');
    $fo->write_seq($aa);
    close(FA);

    open(SIGP, "$signalp_bin $signalp_args $signalp_tmp |") or
	die "Unable to open pipe from $signalp_bin: $!\n";

    my ($pred, $sprob, $aprob, $cprob, $cpos1, $cpos2);
    while (<SIGP>)
    {
		next if /^>/;
		chomp;

		if (/^Prediction:\s+(.*)/)
		{
			$pred = $1;
			next;
		}

		if (/^Signal peptide probability:\s+(\S+)/)
		{
			$sprob = $1;
			next;
		}

		if (/^Signal anchor probability:\s+(\S+)/)
		{
			$aprob = $1;
			next;
		}

		if (/^Max cleavage site probability:\s+(\S+)\s+between pos\.\s+(\S+)\s+and\s+(\S+)/)
		{
			($cprob, $cpos1, $cpos2) = ($1, $2, $3);
			next;
		}
    }

    close(SIGP);
    unlink($signalp_tmp);

    unless (defined $pred  &&
			defined $sprob &&
			defined $cprob &&
			defined $cpos1 &&
			defined $cpos2)
    {
		die "Unable to parse SignalP result for CDS:\n" .
			$fth->make_feature_block($f);
    }

	next if $pred eq 'Non-secretory protein';

	if ($type eq 'euk')
	{
		die "Unable to parse SignalP result (anchor probability) )for CDS:\n" .
			$fth->make_feature_block($f) unless defined $aprob;
	}

	my $note = "";

	if ($type eq 'euk')
	{
		next if (($sprob < $prob) && ($aprob < $prob));

		$note = sprintf("%s predicted for %s by SignalP 2.0 HMM (Signal peptide probabilty %1.3f, signal anchor probability %1.3f) with cleavage site probability %1.3f between residues %d and %d",
						$pred,
						$id,
						$sprob,
						$aprob,
						$cprob,
						$cpos1,
						$cpos2);
	}
	else
	{
		next if $sprob < $prob;

		$note = sprintf("%s predicted for %s by SignalP 2.0 HMM (Signal peptide probabilty %1.3f) with cleavage site probability %1.3f between residues %d and %d",
						$pred,
						$id,
						$sprob,
						$cprob,
						$cpos1,
						$cpos2);
	}

    my $sigpf = $anf->annotate_translation($f, 1, $cpos1);
    $sigpf->qadd("note", $note);
    $sigpf->qadd("colour", $colour);

    # Print features as they are created using the handler
    print $fth->make_feature_block($sigpf);
}
