#!/usr/bin/perl
# Description: This script runs the NetPhos 3.1b.ws0 Web Service. It reads a FASTA file from STDIN and produces predictions in a simple table.
# Author: Edita Bartaseviciute
# Email: edita@cbs.dtu.dk
# Version: 3.1b ws0
# Date: 2009-07-16
# usage: perl netphos.pl [-s] [-c float] < example.fsa


use strict;
# include standard XML::Compile helper functions (used to initiate WSDL proxys)
require "xml-compile.pl"; # downloadable from the same site as this script


#taking options from a command line
my ($soption, $cutoff);
for (my  $i = 0; $i  <= $#ARGV; $i ++) {
	if ($ARGV[$i] eq "-s") {
		$soption = defined;
	}
	elsif ($ARGV[$i] eq "-c") {
		$cutoff = $ARGV[$i+1];
	}
}


# create proxy to NetPhos Web Service
my $netphos = WSDL2proxy ( 'http://www.cbs.dtu.dk/ws/NetPhos/NetPhos_3_1b_ws0.wsdl' );

# append schema definitions
$netphos = appendSchemas ( $netphos , 
	"http://www.cbs.dtu.dk/ws/common/ws_common_1_0b.xsd" ,
	"http://www.cbs.dtu.dk/ws/NetPhos/ws_netphos_3_1_ws0.xsd"
);
# create hash of operations from proxy
my %ops = addOperations ( $netphos ) ;

# Get sequence in fasta format from STDIN
my @fasta;
my $entry = -1;

while (<STDIN>) {
	if (/^>(.*)/) {
		my ($id , $comment) = split (" ",$1);
		$entry++;
		$fasta[$entry]->{id} = $id;
		$fasta[$entry]->{comment} = $comment if defined $comment;
	} elsif (/^([A-Za-z]+)/) {
		$fasta[$entry]->{seq} .= $1;
	}
}
# Create sequence for request
my @sequence;
for ( my $i = 0 ; $i < scalar ( @fasta ) ; $i ++ ) {
	push @sequence , { id => $fasta[$i]->{id} , comment => $fasta[$i]->{comment} , seq => $fasta[$i]->{seq} };
}
# Do the request
my $response;
if ((defined $soption) and (defined $cutoff)) {
	$response = $ops{runService}->(
	 parameters => {
	  parameters => {
	   ser => 'required' ,
	   cutoff => $cutoff ,
	    sequencedata => {sequence => [@sequence]} } });
}
elsif ((defined $soption) and !(defined $cutoff)) {
	$response = $ops{runService}->(
	 parameters => {
	  parameters => {
	   ser => 'required' ,
	     sequencedata => {sequence => [@sequence]} } });
}
elsif (!(defined $soption) and (defined $cutoff)) {
	$response = $ops{runService}->(
	 parameters => {
	  parameters => {
	   cutoff => $cutoff ,
	    sequencedata => {sequence => [@sequence]} } });
}
else {
	$response = $ops{runService}->(
	 parameters => {
	  parameters => {
	    sequencedata => {sequence => [@sequence]} } });
}

# uncomment the two following lines to inspect the structure of $response
# use Data::Dumper;
# print Dumper($response);   

#get job id which can be used to get the results later
my $jobid;

if ( ! defined ( $response->{parameters}->{queueentry}) ) {
	die "error obtaining jobid\n";
} else {
	$jobid = $response->{parameters}->{queueentry}->{jobid};
	print STDERR "# waiting for job $jobid";
	my $status = "UNKNOWN";;
	# poll the queue
	while ( $status =~ /ACTIVE|RUNNING|QUEUED|WAITING|PENDING|UNKNOWN/ ) {
		my $response = $ops{pollQueue}->( job => { job  => { jobid => $jobid }   }) ;
		$status = $response->{queueentry}->{queueentry}->{status};
		print STDERR  ".";
	}
	die "\nunexpected job status '$status'\n" unless $status eq "FINISHED";
	print STDERR "\n# job has finished\n";
}
# when the job is done, fetch the result
$response = $ops{fetchResult}->(job => { jobid => $jobid });

# uncomment the two following lines to inspect the structure of $response
# use Data::Dumper;
# print Dumper($response); 

#printing the results (suitable for one sequence)
foreach my $ann (@{$response->{parameters}->{anndata}->{ann}}) {
	my $sequence = $sequence[0]->{seq};
	my $length = length ($sequence);
	print ">$ann->{sequence}->{id}\t$length amino acids\n";	
	print "#\n# netphos-3.1b prediction results\n#\n";
	print "# Sequence                 # x   Context     Score   Kinase    Answer\n",
		"# -------------------------------------------------------------------";
	my $prev_pos = "";
	foreach my $annrecord (@{$ann->{annrecords}->{annrecord}}) {
		my ($context, $answer);
		my $pos = sprintf ("%4s", $annrecord->{pos});
		my $score = sprintf ("%.3f", $annrecord->{score}[0]->{value});
		my $aa = substr ($sequence, $pos - 1, 1);
		
		#getting sequence context of the analysed residue
		my $start_pos = $annrecord->{pos} - 5;
		if ($start_pos < 0) {
			my $n = 0;
			for (my $m = $start_pos; $m < 0; $m++) {
				$context .= "-";
				$n++;
			}
			$n = 9 - $n;
			$context .= substr ($sequence, 0, $n);
		}
		elsif ($annrecord->{pos} + 4 > $length) {
			my $n = $length - $start_pos;
			$context .= substr ($sequence, $start_pos, $n);
			for (my $m = $n; $m < 9; $m++) {
				$context .= "-";
			}
		}
		else {
			$context = substr ($sequence, $start_pos, 9);
		}
		#getting the content for 'answer' column
		if (exists $annrecord->{comment}) {
			$answer = "YES";
		}
		else {
			$answer = " . ";
		}
		#printing result lines
		if ($pos ne $prev_pos) {
			print "#\n# $ann->{sequence}->{id}\t\t$pos $aa   $context   $score   $annrecord->{feature}\t$answer\n";
		}
		else {
			print "# $ann->{sequence}->{id}\t\t$pos $aa   $context   $score   $annrecord->{feature}\t$answer\n";
		}
		$prev_pos = $pos;
	}
	print "#\n";
}