#!/usr/bin/perl -w
use strict;
use Bio::SearchIO;
my$egDirTemp="/tmp/isa_temp";
my$egDirOutput="isa_output";
@ARGV>=8 or die""."Error in syntax!\n"."   ./blast_profile.pl <-d[n/p] DataBase.fasta> <-q[n/p] QuerySeq.fasta> <-o Output.TabDelimited.dat> <-e EValueCutOff> [-dirtemp DirTemp] [-diroutput DirOutput] [-blast blast_program] [-displaynumber 10000]\n"."Note:\n"."   Directory information in the output file name will be discarded! And the output file will be placed\n"."in subdirectory \"$egDirOutput\" of the current directory.\n"."   And you could specify the blast program as \"blastn\".\n";
my$egDirScript="/hphome/www/html/ffzhou/MUST/MUST";
my$egBlastType="";
my$egMaxDisplayNumber=10000;
my$tBlastP="";
my$tline="";
my@tlist=();
my$tempi=0;
my$tempj=0;
my$tempk=0;
my$pline="";
my@plist=();
my$tid=0;
my%tidx=();
my$egFileD="";
my$egTypeD="";
my$egFileQ="";
my$egTypeQ="";
my$egFileO="";
my$egEValueCutOff="";
my$egFileMatrix="BLOSUM62";
if(&efSyntax()==1){}else{die"-Error in syntax!!----\n";}my$egCmdLine="";
my$egOldTime=`date +%s`;$egOldTime=~s/[\r\n]//g;
my$egNowTime=`date +%s`;$egNowTime=~s/[\r\n]//g;
print"Blasting the query sequences against the database ... ";
&efBlast();
$egNowTime=`date +%s`;$egNowTime=~s/[\r\n]//g;
print" [done] [".($egNowTime-$egOldTime)." seconds]\n";
$egOldTime=$egNowTime;
print"Analyzing the BlastAll results [$egDirTemp\/$egFileO\.$egBlastType] ... ";
&efAnalyze();
if(("$egDirTemp/$egFileO" ne"$egDirOutput/$egFileO")and system("cp \"$egDirTemp\/$egFileO\" \"$egDirOutput\/\"")){die"-Error when copying the result file back!\n";}$egNowTime=`date +%s`;$egNowTime=~s/[\r\n]//g;
print" [Output:$egDirOutput\/$egFileO] [done] [".($egNowTime-$egOldTime)." seconds]\n";
$egOldTime=$egNowTime;
sub efAnalyze{my%tStrand=(-1=>"-",0=>"NA",1=>"+");
open(efOut,">$egDirTemp\/$egFileO")or die"Error when saving data!\n";
print efOut"#NameQ	LengthQ	TypeQ	EValue	Score	Identical%	Conserved%	NameH	LengthH	TypeH	StartAQ	EndAQ	LengthAQ	StrandAQ	StartAH	EndAH	LengthAH	StrandAH	NumIdentical	NumConserved\n";
my$aSearchIO=new Bio::SearchIO(-file=>"$egDirTemp\/$egFileO\.$egBlastType",-format=>"blast");
while(my$aResult=$aSearchIO->next_result){my$aNameQ=$aResult->query_name();
my$aLengthQ=$aResult->query_length();
while(my$aHit=$aResult->next_hit){my$aNameH=$aHit->name();
my$aLengthH=$aHit->length();
while(my$aHSP=$aHit->next_hsp){my$aLengthAQ=$aHSP->length('query');
my$aStartAQ=$aHSP->start('query');
my$aEndAQ=$aHSP->end('query');
my$aStrandAQ=$aHSP->strand('query');
my$aLengthAH=$aHSP->length('hit');
my$aStartAH=$aHSP->start('hit');
my$aEndAH=$aHSP->end('hit');
my$aStrandAH=$aHSP->strand('hit');
my$aNumIdentical=$aHSP->num_identical;
my$aNumConserved=$aHSP->num_conserved;
my$aEValue=$aHSP->evalue;
my$aScore=$aHSP->score;
print efOut"$aNameQ	$aLengthQ	$egTypeQ	$aEValue	$aScore	".($aNumIdentical/$aLengthQ)."	".($aNumConserved/$aLengthQ)."	$aNameH	$aLengthH	$egTypeD	"."$aStartAQ	$aEndAQ	$aLengthAQ	".$tStrand{$aStrandAQ}."	$aStartAH	$aEndAH	$aLengthAH	".$tStrand{$aStrandAH}."	$aNumIdentical	$aNumConserved\n";}}}close(efOut);}sub efBlast{if(system("mkdir -p \"$egDirTemp\"")or system("mkdir -p \"$egDirOutput\"")){die"Error when creating the data directories!\n";}print" [Dir:Created]";
$egCmdLine="";
if((-e"$egFileD.nsq")or(-e"$egFileD.psq")or(-e"$egFileD.00.nsq")or(-e"$egFileD.00.psq")){print" [DatabaseAlreadyFormatted]";}elsif($egTypeD eq"N"){$egCmdLine="\"$egDirScript/formatdb\" -p F -i \"$egFileD\"";}else{$egCmdLine="\"$egDirScript/formatdb\" -p T -i \"$egFileD\"";}if((length($egCmdLine)>0)and(system($egCmdLine))){die"-Error when formatdb database!--\n";}print" [FormatDB]";
if(($egTypeD eq"N")and($egTypeQ eq"N")){$egBlastType="tblastx";}elsif(($egTypeD eq"N")and($egTypeQ eq"P")){$egBlastType="tblastn";}elsif(($egTypeD eq"P")and($egTypeQ eq"N")){$egBlastType="blastx";}elsif(($egTypeD eq"P")and($egTypeQ eq"P")){$egBlastType="blastp";}else{die"--Error in sequence types in database or query file!-\n";}if(length($tBlastP)>0){$egBlastType=$tBlastP;}print" [Type:$egBlastType]";
if(system("\"$egDirScript/blastall\" -W 8 -v $egMaxDisplayNumber -b $egMaxDisplayNumber -p $egBlastType -F F -e $egEValueCutOff -d \"$egFileD\" -i \"$egFileQ\" -o \"$egDirTemp\/$egFileO\.$egBlastType\"")){die"--Error when running BlastAll [$egBlastType]!\n";}print" [BlastAll]";}sub efSyntax{$tBlastP="";
for($tempi=0;$tempi<int(@ARGV/2);$tempi++){if($ARGV[$tempi*2]=~/^\-d(.)$/i){$tid="\U$1\E";
if(($tid eq"N")or($tid eq"P")){$egTypeD=$tid;
$egFileD=$ARGV[$tempi*2+1];}}elsif($ARGV[$tempi*2]=~/^\-q(.)$/i){$tid="\U$1\E";
if(($tid eq"N")or($tid eq"P")){$egTypeQ=$tid;
$egFileQ=$ARGV[$tempi*2+1];}}elsif($ARGV[$tempi*2]=~/^\-e$/i){$egEValueCutOff=$ARGV[$tempi*2+1];}elsif($ARGV[$tempi*2]=~/^\-o$/i){$tline=$ARGV[$tempi*2+1];
@tlist=split(/\//,$tline);
$egFileO=$tlist[@tlist-1];}elsif($ARGV[$tempi*2]=~/^\-dirtemp$/i){$egDirTemp=$ARGV[$tempi*2+1];}elsif($ARGV[$tempi*2]=~/^\-diroutput$/i){$egDirOutput=$ARGV[$tempi*2+1];}elsif($ARGV[$tempi*2]=~/^\-blast$/i){$tBlastP=$ARGV[$tempi*2+1];}elsif($ARGV[$tempi*2]=~/^\-displaynumber$/i){$egMaxDisplayNumber=$ARGV[$tempi*2+1];}}if($egMaxDisplayNumber=~/^\d+$/){}else{$egMaxDisplayNumber=10000;}if((&efValidString($egFileD)==1)and(&efValidString($egTypeD)==1)and((-e$egFileD)or(-e"$egFileD.nsq")or(-e"$egFileD.psq")or(-e"$egFileD.00.nsq")or(-e"$egFileD.00.psq"))and(&efValidString($egFileQ)==1)and(&efValidString($egTypeQ)==1)and(-e$egFileQ)and(&efValidString($egFileO)==1)and(-e$egFileMatrix)){return 1;}else{return 0;}}sub efValidString{my($vA)=@_;
if((!(defined$vA))or(length($vA)==0)){return 0;}else{return 1;}}
