#!/bin/sh -- perl
eval 'exec /usr/bin/perl -w -S $0 ${1+"$@"}'
    if 0;

# PIR, March 2003,
# Modified May, 2003 

# File name: pirsf.pl
# Usage: perl pirsf.pl seq_file > out_file
# PIR SF scan program for one sequence in fasta format

#settings
#PIRSF data file
$pirsf_dat="pirsf.dat";
#HMMER program
$hmmpfam="/drive1/WWW/bin/hmmpfam";
#HMM models
$dm_hmm="dm_hmm";
$sf_hmm="sf_hmm";

#Read in pirsf.dat 
open(IN,$pirsf_dat);
while($line=<IN>)
{ chop $line;
  if ($line=~/^\>/) 
   {($a,$b)=(split /\s+/,$line)[0,1]; 
     $a=~s/\>//; 
     $b1="";
     if ($b)
      {foreach $x ((split /:/,$b))
        { if (index($x,"(")>0)
           { $x=~s/\)/\}/;$x=~s/\(/\)\{/; $x=~s/\-/\,/;
             $b1.="(".$x;
            }
          else {$b1.=$x;}
        }
      }
     $doma{$a}=$b1;
   }
  else {$name{$a}=$line;} 
}
close IN;

#Run HMM search for domain models and get information
@dm_out=` $hmmpfam --cut_tc $dm_hmm $ARGV[0]`;
foreach $x (@dm_out)
{ if ($x=~/^PIRSF/ && index ($x,"\/")>0 ) 
   { $x=~s/-/ /;chop $x; 
     ($a,$b,$e)=(split / /,$x)[0,1,-1]; 
     $domain{$a}.="$b";
     $domain_e{$a}+=$e;
     $domain_n{$a}++;
     $domain_out{$a}.="$x\n";
   }
  elsif (index($x,"seq-f seq-t")>0) { $dm_line=$x; }
}
%temp=%domain;
foreach $u (keys %temp)
   { foreach $v (keys %temp)
      {if (($temp{$u} ne $temp{$v}) && index($temp{$u},$temp{$v})>=0) {delete $domain{$v};}}
   }

#Run HMM search for full-length models and get information
@sf_out=` $hmmpfam --cut_tc --acc $sf_hmm $ARGV[0]`;
foreach $x (@sf_out)
{ if ($x=~/^PIRSF/ && index ($x,"\/")>0 )  
     { chop $x; 
       ($a,$e)=(split / /,$x)[0,-1];
       $sf{$a}=1;
       $sf_e{$a}=$e;
       $sf_out{$a}="$x\n"; 
     }
  elsif ($x=~/^Query/)
    { $id=(split / /,$x)[2]; chop $id;}
  elsif (index($x,"seq-f seq-t")>0) { $sf_line=$x; }
}

#Check domains and calculate mean e-value
foreach $x (sort keys %sf) 
{ #if ($domain{$x}) { chop $domain{$x};}
  if (!$doma{$x} && !$domain{$x}) { $e{$x}=$sf_e{$x};}
  elsif ($doma{$x} && $domain{$x} && $domain{$x}=~/^$doma{$x}$/)   
    { $e{$x}=($sf_e{$x}+$domain_e{$x}/$domain_n{$x})/2;}
} 

#find the min. mean e-value
$min=100;
foreach $x (sort keys %e)
{ if ($e{$x} < $min )
   { $min=$e{$x}; $SF=$x;}
}

#Output
if ($SF) 
{$out=" matches $SF: $name{$SF}\n
Full-length match:\n$sf_line\n$sf_out{$SF}\nDomain match:\n$dm_line\n$domain_out{$SF}";}
else { $out="No Match";}
print "Query sequence: $id $out\n";