1#! /usr/bin/perl
2
3$top_builddir  = shift;
4$top_srcdir    = shift;
5$wrkdir        = shift;
6$tblfile       = shift;
7$msafile       = shift;
8$fafile        = shift;
9$outfile       = shift;
10
11$phmmer     = "$top_builddir/src/phmmer";
12$searchopts = "-E 200 --cpu 1";
13
14if (! -d $top_builddir)  { die "didn't find build directory $top_builddir"; }
15if (! -d $top_srcdir)    { die "didn't find source directory $top_srcdir"; }
16if (! -x $phmmer)        { die "didn't find executable $jackhmmer"; }
17if (! -e $wrkdir)        { die "$wrkdir doesn't exist"; }
18
19open(OUTFILE,">$outfile") || die "failed to open $outfile";
20open(TABLE, "$tblfile")   || die "failed to open $tblfile";
21MSA:
22while (<TABLE>)
23{
24    ($msaname) = split;
25
26    $cmd  = "esl-afetch -o $wrkdir/$msaname.sto $msafile $msaname";                                      $output = `$cmd`;     if ($?) { print "FAILED: $cmd\n"; next MSA; }   # Fetch the query MSA from the benchmark; tmp .sto file here
27    $cmd  = "esl-seqstat --amino -a $wrkdir/$msaname.sto | grep '^=' | awk '{print \$2}'";               $output = `$cmd`;     if ($?) { print "FAILED: $cmd\n", next MSA; }   # Extract list of indiv seq names. --amino for robustness, some msa's v. small
28    @qnames = split(/^/,$output);
29    chop (@qnames);
30    $qname = $qnames[0];
31    $cmd = "esl-sfetch -o $wrkdir/$msaname.query $wrkdir/$msaname.sto $qname > /dev/null";                         `$cmd`;     if ($?) { print "FAILED: $cmd\n"; next MSA; }   # Pick a single seq (first one) to tmp file; tmp .query file here
32
33    $cmd = "$phmmer $searchopts --tblout $wrkdir/$msaname.tmp $wrkdir/$msaname.query $fafile > /dev/null";         `$cmd`;     if ($?) { print "FAILED: $cmd\n"; next MSA; }   # phmmer against benchmark db; tmp .tmp output file here
34
35    if (! open(OUTPUT, "$wrkdir/$msaname.tmp")) { print "FAILED: to open $wrkdir/$msaname.tmp"; next MSA; }
36    while (<OUTPUT>)
37    {
38	if (/^\#/) { next; }
39	@fields   = split(' ', $_, 7);
40	$target   = $fields[0];
41	$pval     = $fields[4];
42	$bitscore = $fields[5];
43	printf OUTFILE "%g %.1f %s %s\n", $pval, $bitscore, $target, $msaname;
44    }
45
46    close OUTPUT;
47    unlink "$wrkdir/$msaname.tmp";
48    unlink "$wrkdir/$msaname.query";
49    unlink "$wrkdir/$msaname.sto";
50}
51close TABLE;
52close OUTFILE;
53