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