#!/usr/local/bin/perl -w

# launch MZEF on all fasta sequence present in the directory given as parameter

# MZEF server URL
$MZEF = "http://sciclio.cshl.org/genefinder/ARAB/mzef.cgi";
$POST = "/usr/local/bin/POST";
$HTML2TXT = "/usr/local/bin/html2txt";

if (!-x $POST)
{
    die "need program $POST";
}
if (!-x $HTML2TXT)
{
    die "need program $HTML2TXT";
}

#------------------------------------------------------------

sub complement
{
    if ($_[0] eq "A") {return("T");}
    if ($_[0] eq "T") {return("A");}
    if ($_[0] eq "C") {return("G");}
    if ($_[0] eq "G") {return("C");}
    return($_[0]);
}

#------------------------------------------------------------


if ($#ARGV < 1)
{
    die "usage: param1 = directory where FASTA sequence files are; param2 = file containing a list of sequence name and length, param=value... where param is overlap or prior_proba";
}

if (!-d $ARGV[0])
{
    die "$ARGV[0] is not a directory";
}

$fasta_dir = $ARGV[0];

if ($fasta_dir !~ /\/$/)
{
    $fasta_dir .= "/";
}

# retrieve fasta files
@fasta_files = glob "${fasta_dir}*.tfa";

if ($#fasta_files == -1)
{
    die "not FASTA files (*.tfa) in $fasta_dir";
}

if (!-f $ARGV[1])
{
    die "$ARGV[1] is not a file";
}

# read sequence name to be computed
open (F_IN, "$ARGV[1]") ||
    die "can't open $ARGV[1]";
while (<F_IN>)
{
    ($seq_id) = ($_ =~ /seq(\d+)\s+.+/);
    $seq_id = int $seq_id;
    $good_seq{$seq_id} = 1;
}
close F_IN;

# get PID for temp file
$temp_file = "$$.form_out";

# loop on FASTA files
for ($seq_no=0;$seq_no<=$#fasta_files;$seq_no++)
{
    # retrieve sequence number
    ($seq_id) = ($fasta_files[$seq_no] =~ /seq(\d+)/);
    $seq_id = int $seq_id;

    # check if this sequence has to be computed
    if (!exists $good_seq{$seq_id})
    {
	next;
    }

    print STDOUT "compute seq $seq_id\n";

    # read the sequence file
    if (!-r $fasta_files[$seq_no])
    {
	die "can't read file $fasta_files[$seq_no]";
    }

    open (F_IN, "$fasta_files[$seq_no]") ||
	die "can't open $fasta_files[$seq_no]";

    $seq = "";
    $first=0;
    while (<F_IN>)
    {
	if ($first == 0)
	{
	    $first=1;
	    next;
	}
	$_ =~ s/[\s\n\r\f]//g;
	$seq .= uc($_);
    }
    close F_IN;

    $seq_length = length($seq);

    # set default values for overlap and prior_proba
    $param{'overlap'} = 0;
    $param{'prior_proba'} = '.03';

    # set user's values for overlap and prior_proba
    for ($i=2; $i <= $#ARGV; $i++)
      {
	if (($p,$v) = ($ARGV[$i] =~ /([^=]+)=(.*)/))
	  {
	    $param{$p} = $v;
	  }
      }

    # launch program for each strand
    for ($strand=1; $strand<3; $strand++)
    {
	print STDOUT "\tstrand $strand\n";

	# generate form output file
	open (F_OUT, ">$temp_file") ||
	    die "can't create temporary file $temp_file";

	print F_OUT "parm1=$strand";
	printf F_OUT "&parm2=%s", $param{'overlap'};
	printf F_OUT "&parm3=%s", $param{'prior_proba'};
	print F_OUT "&address=";
	print F_OUT "&parm4=";
	# add a sequence name otherwize the first 60 bases are skipped
	print F_OUT "%3Eseq$seq_id%0D%0A";
	# split sequence into lines of 60 bases, otherwise MZEF returns nothing
	for ($i=0;$i<$seq_length;$i+=60)
	{
	    print F_OUT substr($seq, $i, 60);
	    print F_OUT "%0D%0A";
	}
	close F_OUT;

	# launch the program
	$cmd = sprintf "cat $temp_file|$POST -e $MZEF | $HTML2TXT > seq%02dmzef%d.txt", $seq_id, $strand;
	system "$cmd";
    }
}

unlink "$temp_file";
