#! /usr/bin/perl -w

######################################################################
# Predict beta pairing map for a singel sequence with user defined 
# secondary structure.
# Input: predict_ssa.sh, beta predictor, 3-line input file(>name, seq, ss), output file
# Output foramt (name, seq, ss, bp1(all 0), bp2(all 0), matrix)

#copy from predict_fasta_beta.pl (to use user defined ss instead of prediction)
# Author: Jianlin Cheng, 3/20/2004
#######################################################################

if (@ARGV != 4)
{
	die "need 4 parameters: ssa predictor, beta predictor, input file (>name, seq, ss), output file.\n"; 
}
$ssa_predictor = shift @ARGV;
$beta_predictor = shift @ARGV; 
$fasta_file = shift @ARGV;
$output_file = shift @ARGV; 

if (! -f $ssa_predictor)
{
	die "can't find the ss, sa predictor.\n"; 
}
if (! -f $beta_predictor)
{
	die "can't find beta predictor.\n"; 
}
if (! -f $fasta_file)
{
	die "can't find the fasta file.\n"; 
}
open(FASTA, "$fasta_file") || die "can't open fasta file.\n"; 
$target_name = <FASTA>; 
$user_seq = <FASTA>;
$user_ss = <FASTA>;
$user_ss =~ s/\s//g;
close FASTA;
#bug fix. Must remove SS before use blast, otherwise, blast treat it as AA.
open(FASTA, ">$fasta_file") || die "can't overwite fasta file.\n";
print FASTA "$target_name$user_seq";
close FASTA;
$target_name = substr($target_name,1); 

#construct a temporay file for prediction of ss and sa
$pos = rindex($fasta_file, "/"); 
if ($pos < 0)
{
	$ssa_file = $fasta_file; 
}
else
{
	$ssa_file = substr($fasta_file, $pos+1); 
}
$ssa_file =~ s/\./_/g; 
$ssa_file .= "ssa"; 
$align_file = "${ssa_file}align"; 
`$ssa_predictor $fasta_file $ssa_file`; 
#notice: two files are generated from ssa predictor: one is ssa output, one is alignment file. 

#create tmp file for beta prediction
open(TMP, ">$output_file.bak") || die "can't create temporary file.\n"; 
open(INPUT, "$ssa_file") || die "can't open the ssa file.\n"; 
<INPUT>;
$seq = <INPUT>; 
chomp($seq); 
$ss = <INPUT>; 
chomp($ss); 
$sa = <INPUT>; 
chomp($sa); 
close INPUT; 
$length = length($seq); 
if ($length != length($ss) || $length != length($sa))
{
	die "sequence length doesn't match.\n"; 
}

print TMP "$target_name$seq\n$user_ss\n$sa\n";

#predict beta-mapping
system("$beta_predictor $output_file.bak $align_file $output_file"); 

#remove temporary files
`rm $output_file.bak`; 
`rm $ssa_file`; 
`rm $align_file`;