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