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

# PIR, March 2003,
# Modified Oct, 2004 
# Major changes:

# 1. Take off domain match.
# 2. Add sub-families
# 3. Replace cut-off by the following conditions:
#    L(matched)/Length >= 0.8
#     and
#    Score - mean(S) >= -Std(S)  or  (Score - mean(S)  >= -2.5Std(S)) and (Score >= min(S)) 
#     and
#    |Length-mean(L)| < 3.5Std(L)  or  |Length-mean(L)| < 50
#
#    where Length is the query seq length and Score is the query HMM score.
#          mean(L) and std(L) are for the length of the family
#          mean(S), std(S) and min(S) are for the HMM score of the family 


# Modified April 2005
# Major changes:

# 1. condition: Score - mean(S) >= -Std(S)  or  (Score - mean(S)  >= -2.5Std(S)) and (Score >= min(S))
#     becomes  Score >= min(S)
# 2. if a sequence can be placed into a subfamily then it will be placed its parent PIRSF
# 3. add BLAST search for most PIRSF:  predicted sequence is checked by BLAST for the best and "majority" hit to the PIRSF.


########################################################################
# 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";
#BLAST program
$blastall="/drive1/WWW/bin/blastall";
#HMM models
$sf_hmm_subf="sf_hmm_subf";
$sf_hmm="sf_hmm";
#sf.tb file
$sftb="sf.tb";
#sf.seq file
$sfseq="sf.seq";

#Read in pirsf.dat 
open(IN,$pirsf_dat);
while($line=<IN>)
{ chop $line;
  if ($line=~/^\>/) 
   {  $x=(split / /,$line)[0];
      $x=~s/\>//;
      if (index($line,"child")>0) { $child{$x}=(split /: /,$line)[1]; }
      $line=<IN>; chop $line;
      $name{$x}=$line;
      $line=<IN>; chop $line;
      ($lm{$x},$lsd{$x},$smin{$x},$sm{$x},$ssd{$x})=split (/ /,$line);
      # lm=mean(L); lsd=Std(L); smin=min(S); sm=mean(S); ssd=Std(S);
      $line=<IN>; chop $line;
      if ($line=~/Yes/) { $blast{$x}=1; }
   }
}
close IN;

#get seq length
$infile=$ARGV[0];
open(IN,$infile);
while($line=<IN>)
{ chop $line;
  if ($line=~/^\>/){ $leng=0;}
  else { $line=~s/\s+//g; 
         $leng= $leng + length($line);} 
}
close IN;

#Run HMM search for full-length models and get information
@sf_out=` $hmmpfam -E 0.01 --acc $sf_hmm $infile`;

foreach $x (@sf_out)
{ if ($x=~/^PIRSF/ && index ($x,"\/")>0 )  
     { chop $x; 
       ($a,$b,$c,$s,$e)=(split /\s+/,$x)[0,2,3,8,9];
       $ovl=($c-$b+1)/$leng;
       $ld=abs($leng-$lm{$a});
       $sf{$a}=1;
       $sf_out{$a}="$x\n";
       if($ovl>=0.8 && ($s>=$smin{$a}) && ($ld<3.5*$lsd{$a}||$ld< 50 ))
        { 
          $e{$a}=$e;
        }
     }
  elsif ($x=~/^Query/)
    { $id=(split / /,$x)[2]; chop $id;}
  elsif (index($x,"seq-f seq-t")>0) { $sf_line=$x; }
}

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


#sub-family search
foreach $z (keys %sf) 
{if ($child{$z})
 { foreach $x ((split / /,$child{$z}))
     { $cf{$x}=1;  $subf=1;
       $parent{$x}=$z; 
     }
 }
}
if ($subf)
 { @sf_out_subf=` $hmmpfam -E 0.01 --acc $sf_hmm_subf $infile`;
   foreach $x (@sf_out_subf)
       { if ($x=~/^PIRSF/ && index ($x,"\/")>0 )
          { chop $x;
            ($a,$s,$e)=(split /\s+/,$x)[0,8,9];
            if($cf{$a} && $s>=$smin{$a})
             { 
               $ce{$a}=$e;
               $sf_out_subf{$a}="$x\n";
             }
          }
        elsif (index($x,"seq-f seq-t")>0) { $sf_line_subf=$x; }
       }
   #find the min. mean e-value
    $min=100;
    foreach $x (sort keys %ce)
     { if ($ce{$x} < $min )
        { $min=$ce{$x}; $SF_subf=$x;}
     }
 }
if ($SF_subf && !$SF) {$SF=$parent{$SF_subf};}


#BLAST search
if ($SF && $blast{$SF})
{
  @blast_tmp=`blastall -p blastp -F F -e 0.0001 -b 10000 -d $sfseq -i $infile -m 8`; 
  open (IN,$sftb);
  while($line=<IN>)
   { chop $line;
     ($a,$b)=(split / /,$line)[0,1];
     $SFn{$a}=$b;
   }
  close IN;

  %n=();
  foreach $x (@blast_tmp)
  { 
  $y=(split /\s+/,$x)[1]; 
  $n{$y}++;
  if ($n{$y}<2) 
   {$sf=(split /-/,$y)[1];
    if (!$sf1){$sf1=$sf;}
    if ($sf eq $sf1 ){ $hit++; }
   }
  }

 foreach $x ( keys %n)
  { if ($hit>9 || ($SFn{$sf1} && $hit/$SFn{$sf1}>0.33334))
    { $bl_sf="PIR$sf1";}
  }
if ($bl_sf ne $SF) {$SF="";$SF_subf="";}
}


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

if ($SF_subf)
{print "and matches Sub-Family $SF_subf: $name{$SF_subf}\n
$sf_line_subf\n$sf_out_subf{$SF_subf}";}