1#! /usr/bin/perl -w
2
3BEGIN {
4    $top_builddir  = shift;
5    $top_srcdir    = shift;
6    $wrkdir        = shift;
7    $tblfile       = shift;
8    $msafile       = shift;
9    $fafile        = shift;
10    $outfile       = shift;
11}
12use lib "${top_srcdir}/easel/demotic";
13use demotic_fasta;
14
15$fasta      = "${top_builddir}/bin/ssearch36";
16$opts       = "-q -E 200 -d 0 -H";
17# -q     = quiet (batch) mode
18# -E 200 = report top hits deeper into noise, down to E=200 (default was 10)
19# -d 0   = suppresses alignment output; we only need the hit list
20# -H     = suppresses histogram output.
21
22if (! -d $top_builddir)                                 { die "didn't find FASTA/SSEARCH build directory $top_builddir"; }
23if (! -d $top_srcdir)                                   { die "didn't find H3 source directory $top_srcdir"; }
24if (! -x $fasta)                                        { die "didn't find executable $fasta"; }
25if (! -e $wrkdir)                                       { die "$wrkdir doesn't exist"; }
26
27open(OUTFILE,">$outfile") || die "failed to open $outfile";
28open(TABLE, "$tblfile")   || die "failed to open $tblfile";
29MSA:
30while (<TABLE>)
31{
32    ($msaname) = split;
33
34    $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
35    $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
36    @qnames = split(/^/,$output);
37    chop (@qnames);
38    $qname = $qnames[0];
39    $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
40
41    # search it against the benchmark db
42    $cmd = "$fasta $opts $wrkdir/$msaname.query $fafile |";                           if (! open(FASTA, "$cmd"))            { print "FAILED: $cmd\n", next MSA; }
43
44    if (! demotic_fasta::parse(\*FASTA)) { print "FAILED: demotic fasta parser on $msaname\n"; next MSA; }
45    for ($i = 0; $i < $demotic_fasta::nhits; $i++)
46    {
47	$target   = $demotic_fasta::hit_target[$i];
48	$pval     = $demotic_fasta::hit_Eval[$i];
49	$bitscore = $demotic_fasta::hit_bitscore[$i];
50	printf OUTFILE "%g %.1f %s %s\n", $pval, $bitscore, $target, $msaname;
51    }
52    close FASTA;
53    unlink "$wrkdir/$msaname.query";
54    unlink "$wrkdir/$msaname.sto";
55}
56close TABLE;
57close OUTFILE;
58