#!/usr/bin/perl
use Cwd;
####################################################
# This script is used to demonstrate how to use "system" call to run other tools, so that you can automate data processing when you need to run multiple tools

#Author: Jack Min  Jan. 2009
#Program:  secPipelineW.pl
# this version to use whole data set for every step
# This pipeline runs serveral programs to screen sereted proteins 
# The input is fasta sequences in a directory, then all the output will be in the same 
# directory.  The protocol:
# (1) to split multiple fasta sequence file into individual files
# (2) run signalp, then tmhmm, phobius, targetp, ps-scan
################################################
if ($#ARGV !=0)
{
	print "usage: secPipelineW.pl input\n";
	print "			input: ref seq in mult-fas format\n";
	exit;
}

$input = $ARGV[0];
open (INPUT, $input) || die ("Cannot open input");
@inputA = <INPUT>;
close (INPUT);
@outputA = ();
$count = 0;
$flag = 0;
$flag1 = 0;
foreach $line (@inputA)
{
	chomp ($line);

	if ($line =~ /^>/)
	{
		if ($count != 0)
		{

			push (@outputA, ">".$gi."\n".$seq."\n");
			push (@seqA, $gi."\t".$seq."\n");
		$seq = '';
		}
		$query = $line;
		@tmp = split(/\|/, $line);
		$gi = $tmp[1];
		$acc = $tmp[3];
		@tmp2 = split (/\./, $acc);
		$acc = $tmp2[0];
		$fun = $tmp[4];
		@tmp3 = split (/\[/, $fun);
		$fun = $tmp3[0];
		@spA = split(/[\[\]]/, $line);
		$species = $spA[1];
		$species =~ s/ /_/g;
		push (@infA, $gi."\t".$acc."\t".$fun."\t".$species."\n");
		$count++;
		$flag = 1;
	}
	elsif ($line =~ /\w/)
	{
		$seq =$seq.$line;
	}
}
push (@outputA, ">".$gi."\n".$seq."\n");
push (@seqA, $gi."\t".$seq."\n");
$seq = '';
$outF = $input.".oneLine";
open (OUTPUT, ">$outF");
print OUTPUT (@outputA);
close OUTPUT;

$outF2 = $species.".info";
open (OUTPUT, ">$outF2");
print OUTPUT (@infA);
close OUTPUT;
$outF3 = $species.".seq_db";
open (OUTPUT, ">$outF3");
print OUTPUT (@seqA);
close OUTPUT;
###################### split #########
$input = $outF;
$line = '';
print "splitting file ...in progress ...\n";
open (IN, $input) || die "Cannot open file $input\n";
$flag = 0;
$count = 0;
$seqID = '';
$f = 'f';
@idMap = ();
system ("mkdir tmp");
chmod 0777, tmp;
while (<IN>)
{
	$line = $_;
	if ($line =~ /^>/)
	{
		$count++;
		if ($count > 1)
		{
		close (OUT);
		}
		$file = "tmp/".$f.$count.".fas";
		open (OUT, ">$file");
		print OUT ($line);
		chomp ($line);
		$seqID = substr ($line, 1);
		push (@idMap, $f.$count."\t".$seqID."\n");
	}
	else
	{
	 	$subLine = substr($line, 0, 70);
		print OUT ($subLine);
	}

}
close (IN);
close (OUT);
print "splitting file is completed..\n";
print "run signalP ...\n";
#########
$sigalp = "/home/jack/Tools/Others/signalp3/signalp";
$dir = getcwd();
$input = $dir."/tmp/f";
system ("mkdir sp_out");
chmod 0777, sp_out;
$sp_out = ">".$dir."/sp_out/f";
$int = 0;
$ext = ".sp";
$ext2 = ".fas";
$total = $count; #14086;
for ($int = 1; $int <= $total; $int++)
{
	system ("$sigalp -t euk -f summary $input$int$ext2 $sp_out$int$ext");
	# /home/jack/Tools/Others/signalp3/signalp -t euk -f summary test.fas >test.out
	system ("rm $input$int$ext2");
}
system ("rm -r tmp");
#### sp is done, now parsing sp output #####
$sigalParse = "perl /home/jack/Tools/YSU_tools/signalParse_strict.pl";
$input = $dir."/sp_out/f";
system ("mkdir sp_parse");
chmod 0777, sp_parse;
$output = $dir."/sp_parse/f";
$int = 0; 
$ext1 = ".sp";
$ext2 = ".parse";
print "sp parsing ...\n";
for ($int = 1; $int <= $total; $int++)
{
	system ("$sigalParse $input$int$ext1 $output$int$ext2");
}
system ("rm -r sp_out");
######## cat file #########
$sp_dir = $dir."/sp_parse";
opendir (DIR, $sp_dir);
while ($file = readdir DIR)
{
	if ($file =~ /parse/)
	{
		push (@fileA, $file."\n");
	}
}
close DIR;

foreach $afile (@fileA)
{
	chomp $afile;
	$file = $sp_dir."/".$afile;
	open (INPUT, $file) || die ("Cannot open file $file\n\n");
	while (<INPUT>)
	{
        	$line = $_; #$_, default input variable
        	if ($file =~ /.id/)
		{
			push (@idA, $line);
		}
		elsif ($file =~ /.db/)
		{
			push(@dbA, $line);
		}
		else
		{
			push (@sumA, $line);
		}
	}
	close INPUT;
	system ("rm $file");
}
system ("rm -r sp_parse");
$out1 = "sp.summary";
$out2 = "sp.id_info";
$out3 = "sp.db";
open (OUTPUT, ">$out1") || die ("Cannot open for writing file $out1 !\n");
print OUTPUT (@sumA);
close OUTPUT;
open (OUTPUT, ">$out2") || die ("Cannot open for writing file $out2 !\n");
print OUTPUT (@idA);
close OUTPUT;
open (OUTPUT, ">$out3") || die ("Cannot open for writing file $out3 !\n");
print OUTPUT (@dbA);
close OUTPUT;
#############
system "/home/jack/Tools/Others/TMHMM2/bin/tmhmm $outF -short -noplot >hmm.out";
system "perl /home/jack/Tools/YSU_tools/hmmParse.pl hmm.out hmm.parse";
print "hmm is done\n";
#######
system "/home/jack/Tools/Others/phobius/phobius.pl $outF -short >phobius.out";
system "perl /home/jack/Tools/YSU_tools/phobiusParse.pl phobius.out phobius+.id";
print "phobius is done\n";
# note: targetP generates wrong results if the number of sequences is too large - the exact number is not known, but 5992 seq input is not handled properly 
#system "/home/jack/Tools/Others/targetp/targetp -c -N $outF >targetp.out";
# split file
open (IN, $outF) || die "Cannot open file $outF\n";
$flag = 0;
$count = 0;
$seqID = '';
$f = 'f';
system ("mkdir tmp");
chmod 0777, tmp;
while (<IN>)
{
	$line = $_;
	if ($line =~ /^>/)
	{
		$count++;
		if ($count > 1)
		{
		close (OUT);
		}
		$file = "tmp/".$f.$count.".fas";
		open (OUT, ">$file");
		print OUT ($line);
		chomp ($line);
		$seqID = substr ($line, 1);
		push (@idMap, $f.$count."\t".$seqID."\n");
	}
	else
	{
	 	print OUT ($line);
	}
}
close (IN);
close (OUT);
##########
print "run TargetP ...\n";
#########
$targetp = "/home/jack/Tools/Others/targetp/targetp";
$dir = getcwd();
$input = $dir."/tmp/f";
system ("mkdir tp_out");
chmod 0777, tp_out;
$tp_out = ">".$dir."/tp_out/f";
$int = 0;
$ext = ".tp";
$ext2 = ".fas";
$total = $count; #14086;
for ($int = 1; $int <= $total; $int++)
{
	system ("$targetp -c -N $input$int$ext2 >$tp_out$int$ext");
	# /home/jack/Tools/Others/signalp3/signalp -t euk -f summary test.fas >test.out
	system ("rm $input$int$ext2");
}
system ("rm -r tmp");
#### sp is done, now parsing sp output #####
$tp_dir = $dir."/tp_out";
opendir (DIR, $tp_dir);
@fileA = ();
while ($file = readdir DIR)
{
	if ($file =~ /tp/)
	{
		push (@fileA, $file."\n");
	}
}
close DIR;

foreach $afile (@fileA)
{
	chomp $afile;
	$file = $tp_dir."/".$afile;
	open (INPUT, $file) || die ("Cannot open file $file\n\n");
	while (<INPUT>)
	{
        	$line = $_; #$_, default input variable
        	push (@tpA, $line);
	}
	close INPUT;
	system ("rm $file");
}
system ("rm -r tp_out");
$out1 = "tp.out";
open (OUTPUT, ">$out1") || die ("Cannot open for writing file $out1 !\n");
print OUTPUT (@tpA);
close OUTPUT;
###############

system "perl /home/jack/Tools/YSU_tools/targetpParse.pl tp.out tp.parse";
system "perl /home/jack/Tools/Others/ps_scan/ps_scan.pl $outF -p PS00014 -o scan -d /home/jack/Tools/Others/ps_scan/prosite.dat >ps.out";
system "perl /home/jack/Tools/YSU_tools/psParse.pl ps.out ps.parse";
print "ps-scan is done\n";
print "The work is done!!\n";


exit;
