#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
my$tline="";
my@tlist=();
my$tempi=0;
my$tempj=0;
my$tempk=0;
my$pline="";
my@plist=();
my$tid=0;
my%tidx=();
my$egDirScript="/home/xynlab/zhouff/work/script_app";
my$egDirGenomes="/db/ncbi-genomes-bacteria/Bacteria092008";
my$egDate=`date \"+%F %T\"`;$egDate=~s/[\r\n]//g;
my$egDirNow=$0;$egDirNow=~s/[\r\n]//g;
@tlist=split(/\//,$egDirNow);
$egDirNow=~s/\/$tlist[@tlist-1]$//g;
my$egK=10;
my%egCI=("A"=>0,"T"=>1,"C"=>2,"G"=>3,);
my@egC=("A","T","C","G",);
my$egV=4;
my$egMaxID=16*16*16*16*16;
my$egDirDB=".";
my$egFileNodes="nodes.dmp";
my$egFileNames="names.dmp";
my%egName2ID=();
my%egTaxID2ID=();
my@egTaxID=();
my@egPTaxID=();
my@egRank=();
my@egName=();
my$egNum=0;
my%egTaxRank=();
my@egTaxLevel=("species","genus","family","order","class","subphylum","phylum","superkingdom");
my%egTaxLevelIdx=();
for($tempi=0;$tempi<@egTaxLevel;$tempi++){$egTaxLevelIdx{"$egTaxLevel[$tempi]"}=1;}my%egCharIdx=();
my%egCharPairScore=();
my%egCharToAnyScore=();
my$egVoidValue=-1;
my@lVA=();my@lVB=();
my$egGSeq="";
@ARGV==5 or die""."Error in syntax!\n"."./refineBoundary.pl <TemporaryDirectory> <MutationRate> <Genome.fasta> <input.MITE-filter.txt> <output.exp-MITE.txt>\n"."Note:\n"."This program will expand the MITE copy to match the boundaries of another copy\n"."with the maximum \"DR+TIR\" values, which is assumed to be the \"optimal\" MITE.\n";
my($egDirTemp,$egCutOffMutation,$egFileSeq,$egFileI,$egFileO)=@ARGV;
print"---[./refineBoundary.pl ".join(" ",@ARGV)."]---\n";
if(!(-d$egDirTemp)){if(system("mkdir -p \"$egDirTemp\"")){die"-Error!\n";}}my%mGID2ID=();
my%mMITE2ID=();
my@mGID=();
my@mBestCopyID=();
my@mID=();
my@mMITE=();
my@mEScore=();
my@mDR=();
my@mTIR=();
my@mSeq=();
my@mSeqL=();
my@mSeqR=();
my@mLenL=();
my@mLenR=();
my@mS=();
my@mE=();
my@mStrand=();
my@mIdentDR=();
my@mIdentTIR=();
my@mLength=();
my@mMITESeq=();
my@mOffsetS=();
my@mOffsetE=();
my$mN=0;
my@mLine=();
print"Loading the file ... ";
my@egHeadLines=();
open(efI,"$egFileI")or die"Error when reading the input file!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){push(@egHeadLines,$tline);}elsif(@tlist>=23){$mGID[$mN]=$tlist[0];
$mID[$mN]=$mN;
$mMITE[$mN]=$tlist[2];
$tid=$mMITE2ID{"$mMITE[$mN]"};
if(defined$tid){$tid.=";$mN";}else{$tid="$mN";}$mMITE2ID{"$mMITE[$mN]"}=$tid;
$tid=$mGID2ID{"$mGID[$mN]"};
if(defined$tid){$tid.=";$mN";}else{$tid="$mN";}$mGID2ID{"$mGID[$mN]"}=$tid;
$mDR[$mN]=$tlist[7];
$mTIR[$mN]=$tlist[11];
$mIdentDR[$mN]=$tlist[8];
$mIdentTIR[$mN]=$tlist[12];
$mEScore[$mN]=&efEScore($mDR[$mN],$mTIR[$mN],$mIdentDR[$mN],$mIdentTIR[$mN]);
my$tFlanking=int(length($tlist[20])/2)-$mDR[$mN];
$mMITESeq[$mN]=$tlist[21];
$mLength[$mN]=length($mMITESeq[$mN]);
$mSeqL[$mN]=substr($tlist[20],0,$tFlanking+$mDR[$mN]);
$mSeqR[$mN]=substr($tlist[20],$tFlanking+$mDR[$mN],$tFlanking+$mDR[$mN]);
$mSeq[$mN]=$mSeqL[$mN].$mMITESeq[$mN].$mSeqR[$mN];
$mLenL[$mN]=$tFlanking;
$mLenR[$mN]=$tFlanking;
$mOffsetS[$mN]=0;
$mOffsetE[$mN]=0;
$mS[$mN]=$tlist[3];
$mE[$mN]=$tlist[4];
$mStrand[$mN]=$tlist[5];
$mLine[$mN]=$tline;
$mN++;}}close(efI);
print" [done] [MITE Copies:$mN]\n";
my%egMITE2BestID=();
print"Screening for the best MITE copies ... ";
for($tempi=0;$tempi<$mN;$tempi++){$tid=$egMITE2BestID{"$mMITE[$tempi]"};
if(defined$tid){if(($mEScore[$tid]<$mEScore[$tempi])or(($mEScore[$tid]==$mEScore[$tempi])and($mDR[$tid]<$mDR[$tempi]))){$tid=$tempi;}}else{$tid=$tempi;}$egMITE2BestID{"$mMITE[$tempi]"}=$tid;}print" [done]\n";
print"Saving the sequences for alignments ... ";
my$tFileQuery="$egDirTemp/temp-refineBoundary.fasta";
open(efO,">$tFileQuery")or die"E!\n";
for($tempi=0;$tempi<$mN;$tempi++){print efO">$tempi\n$mMITESeq[$tempi]\n";}close(efO);
print" [Copies:$mN]";
my$tFileData="$egDirTemp/temp-refineBoundary-DB.fasta";
open(efO,">$tFileData")or die"E:DB!\n";
@tlist=keys%egMITE2BestID;
for($tempi=0;$tempi<@tlist;$tempi++){$tid=$egMITE2BestID{"$tlist[$tempi]"};
print efO">$tid\n$mMITESeq[$tid]\n";
$pline=$mMITE2ID{"$tlist[$tempi]"};
@plist=split(/;/,$pline);
for($tempj=0;$tempj<@plist;$tempj++){$mBestCopyID[$plist[$tempj]]=$tid;}}close(efO);
print" [MITE:".(scalar@tlist)."]";
print" [done]\n";
my$egEValue=1e-20;
my$tFileLog="$egDirTemp/temp-refineBoundary-LOG.txt";
my$tFileMap="$egDirTemp/temp-refineBoundary-Map.txt";
print"Mapping ... ";
if(system("formatdb -p F -i \"$tFileData\" 2>&1 1>>\"$tFileLog\"")or system("megablast -d \"$tFileData\" -i \"$tFileQuery\" -o \"$tFileMap\" -e $egEValue -m 9 -F F 2>&1 1>>\"$tFileLog\"")){die"-Error when mapping!\n";}print" [done]\n";
print"Screening for the boundary information ... ";
open(efI,"$tFileMap")or die"E:Map!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){}elsif(@tlist>=12){my($qID,$dbID)=($tlist[0],$tlist[1]);
if(($mBestCopyID[$qID]==$dbID)and($qID!=$dbID)){my($qL,$dbL)=($mLength[$qID],$mLength[$dbID]);
my($dbTIR,$dbDR)=($mTIR[$dbID],$mDR[$dbID]);
my$tIdent=$tlist[2];
my($qS,$qE,$dbS,$dbE)=($tlist[6],$tlist[7],$tlist[8],$tlist[9]);
my($offsetS,$offsetE)=(0,0);
if($dbS<=$dbE){$offsetS=$qS-$dbS;
$offsetE=($qL-$qE)-($dbL-$dbE);}else{$offsetS=$qS-($dbL-$dbS)-1;
$offsetE=($qL-$qE)-$dbE+1;}my($qSeqL,$qMITESeq,$qSeqR)=($mSeqL[$qID],$mMITESeq[$qID],$mSeqR[$qID]);
my($qLenL,$qLenM,$qLenR)=(length($qSeqL),length($qMITESeq),length($qSeqR));
my$qSeqAll=$qSeqL.$qMITESeq.$qSeqR;
my$qtMITESeq=substr($qSeqAll,$qLenL+$offsetS,$qLenM-$offsetS-$offsetE);
my$qtDRL=substr($qSeqAll,($qLenL+$offsetS)-$dbDR,$dbDR);
my$qtDRR=substr($qSeqAll,$qLenL+$qLenM-$offsetE,$dbDR);
my$qtTIRL=substr($qtMITESeq,0,$dbTIR);
my$qtTIRR=substr($qtMITESeq,length($qtMITESeq)-$dbTIR,$dbTIR);
my$qtIdentTIR=&efCheckIdentity($qtTIRL,&efDNA_RC($qtTIRR));
my$qtIdentDR=&efCheckIdentity($qtDRL,$qtDRR);
if(($qtIdentDR>=$egCutOffMutation)and($qtIdentTIR>=$egCutOffMutation)){$mOffsetS[$qID]=$offsetS;
$mOffsetE[$qID]=$offsetE;
if($qtIdentDR*length($qtDRL)>$mIdentDR[$qID]*$mDR[$qID]){$mDR[$qID]=length($qtDRL);}$mDR[$qID]=&efGetMaxDR($qSeqL,$qSeqR,$egCutOffMutation);
if($qtIdentTIR*length($qtTIRL)>$mIdentTIR[$qID]*$mTIR[$qID]){$mTIR[$qID]=length($qtTIRL);}}}}}close(efI);
print" [done]\n";
print"Re-generating the summary data ... ";
for($tempi=0;$tempi<$mN;$tempi++){$mLine[$tempi]="";}my$egSeqObj=new Bio::SeqIO(-file=>"$egFileSeq",-format=>"fasta");
while(my$egSeqO=$egSeqObj->next_seq){my($egID,$egSeq)=($egSeqO->id,$egSeqO->seq);
$tline=$mGID2ID{"$egID"};
@tlist=split(/;/,$tline);
for($tempi=0;$tempi<@tlist;$tempi++){$tid=$tlist[$tempi];
{my($tTIR,$tDR)=($mTIR[$tid],$mDR[$tid]);
$mLine[$tid]="$egID	$mID[$tid]	$mMITE[$tid]	".&efGenMITECopy($egSeq,$mS[$tid]+$mOffsetS[$tid],$mE[$tid]-$mOffsetE[$tid],$mStrand[$tid],$tTIR,$tDR,$mLenL[$tid],$mLenR[$tid])."\n";}}}print" [done]\n";
print"Saving the final data ... ";
open(efO,">$egFileO")or die"E:Output!\n";
my($cTotal,$cUpdate)=(0,0);
print efO"".join("\n",@egHeadLines)."\n";
for($tempi=0;$tempi<$mN;$tempi++){print efO"$mLine[$tempi]\n";
$cTotal++;
if(($mOffsetS[$tempi]!=0)or($mOffsetE[$tempi]!=0)){$cUpdate++;}}close(efO);
print" [done] [Total:$cTotal] [Updated:$cUpdate]\n";
sub efGetMaxDR{my($aL,$aR,$aMutation)=@_;
my@lL=split(//,"\U$aL\E");my$lenL=scalar@lL;
my@lR=split(//,"\U$aR\E");my$lenR=scalar@lR;
my$maxDR=30;
my($realDR,$realIdent,$nowIdent)=(0,0,0);
my($li,$lj)=(0,0);
for($li=1;$li<=$maxDR;$li++){$nowIdent=0;
for($lj=0;$lj<$li;$lj++){if($lL[$lenL-$li+$lj]eq$lR[$lj]){$nowIdent++;}}if(($nowIdent>$realIdent)and($nowIdent>=$li*$aMutation)){$realDR=$li;
$realIdent=$nowIdent;}}return$realDR;}sub efGenMITECopy{my($pSeq,$pS,$pE,$pStrand,$pTIR,$pDR,$pLenL,$pLenR)=@_;
my($pDRL,$pDRR)=(substr($pSeq,$pS-$pDR-1,$pDR),substr($pSeq,$pE,$pDR));
my$pMITE=substr($pSeq,$pS-1,abs($pE-$pS)+1);
my$pLength=length($pMITE);
my($pTIRL,$pTIRR)=(substr($pMITE,0,$pTIR),substr($pMITE,$pLength-$pTIR,$pTIR));
my($pFlankingL,$pFlankingR)=(substr($pSeq,$pS-$pDR-$pLenL,$pLenL),substr($pSeq,$pE+$pDR,$pLenR));
my$pIdentDR=&efCheckIdentity($pDRL,$pDRR);
my$pIdentTIR=&efCheckIdentity($pTIRL,&efDNA_RC($pTIRR));
my$pRL="$pS	$pE	$pStrand	$pLength"."	$pDR	$pIdentDR	$pDRL	$pDRR"."	$pTIR	$pIdentTIR	$pTIRL	$pTIRR"."	".&efCheckAT($pDRL.$pDRR).""."	".&efCheckAT($pTIRL.$pTIRR).""."	".&efCheckAT($pMITE).""."	".&efCheckAT($pFlankingL.$pFlankingR).""."	".$pFlankingL.$pDRL.$pFlankingR.""."	".$pFlankingL.$pDRL.$pDRR.$pFlankingR.""."	$pMITE"."	".&efEScore($pDR,$pTIR,$pIdentDR,$pIdentTIR).""."\n";}sub efEScore{my($adDR,$adTIR,$adIdentDR,$adIdentTIR)=@_;
my$adWeightIdent=20;
my$adWeightDR=2;
my$adEScore=$adIdentDR*$adDR*$adWeightDR+$adIdentTIR*$adTIR;
return$adEScore;}sub efCheckAT{my($zzz)=@_;
my@zl=split(//,"\U$zzz\E");
my($zTotal,$zAT)=(0,0);
for(my$zi=0;$zi<@zl;$zi++){if(($zl[$zi]eq"A")or($zl[$zi]eq"T")){$zAT++;}$zTotal++;}my$zATP=0;
if($zTotal>0){$zATP=$zAT/$zTotal;}return$zATP;}sub efCheckIdentity{my($aaa,$bbb)=@_;
my@la=split(//,"\U$aaa\E");
my@lb=split(//,"\U$bbb\E");
my$ii=0;
my($ni,$nt)=(0,0);
for($ii=0;($ii<@la)and($ii<@lb);$ii++){if($la[$ii]eq$lb[$ii]){$ni++;}$nt++;}my$np=0;
if($nt>0){$np=$ni/$nt;}return$np;}sub efGetReadNum{my($aFile)=@_;
my$aFormat="fasta";
if(($aFile=~/\.fastq/)or($aFile=~/\.fq/)){$aFormat="fastq";}my$aSeqObj;
if($aFile=~/\.gz$/){$aSeqObj=new Bio::SeqIO(-file=>"gunzip -c \"$aFile\" | ",-format=>"$aFormat");}elsif($aFile=~/\.bz2$/){$aSeqObj=new Bio::SeqIO(-file=>"bunzip2 -c \"$aFile\" | ",-format=>"$aFormat");}else{$aSeqObj=new Bio::SeqIO(-file=>"$aFile",-format=>"$aFormat");}my$aN=0;
while(my$aSeqO=$aSeqObj->next_seq){$aN++;}return$aN;}sub efGetSeq{my($aSeq,$aLength,$aS,$aE,$aStrand)=@_;
my($as,$ae)=($aS,$aE);
if($as<0){$as=0;}if($ae>=$aLength){$ae=$aLength-1;}my$aa="";
if($aStrand eq"+"){$aa=substr($aSeq,$as,abs($ae-$as)+1);}else{$aa=&efDNA_RC(substr($aSeq,$as,abs($ae-$as)+1));}$aa="-"x($as-$aS).$aa."-"x($aE-$ae);
return$aa;}sub efLoadBLOSUM{my($tFile)=@_;
my$tLine="";my@tList=();
my($tI,$tJ)=(0,0);
my@tCharList=();
if(-e$tFile){open(bIn,"$tFile")or die"Error when loading the BLOSUM matrix!\n";
while(<bIn>){$tLine=$_;$tLine=~s/[\r\n]//g;
@tList=split(/\s+/,$tLine);
if($tLine=~/^#/){}elsif($tLine=~/^ /){$tLine=~s/^\s+//g;
@tCharList=split(/\s+/,$tLine);
for($tI=0;$tI<@tCharList;$tI++){{$egCharIdx{"$tCharList[$tI]"}=1;}}}elsif($tLine=~/^(.)\s+(.*)$/){my($tA,$tL)=($1,$2);
if($tA ne"*"){@tList=split(/\s+/,$tL);
for($tI=0;$tI<@tList;$tI++){if(defined$tCharList[$tI]){if((!(defined$egCharPairScore{"$tA;$tCharList[$tI]"}))and(!(defined$egCharPairScore{"$tCharList[$tI];$tA"}))){$egCharPairScore{"$tA;$tCharList[$tI]"}=$tList[$tI];
$egCharPairScore{"$tCharList[$tI];$tA"}=$tList[$tI];
if($tCharList[$tI]eq"*"){$egCharToAnyScore{"$tCharList[$tI]"}=$tList[$tI];}}}}}}}close(bIn);}}sub efBLOSUMScore{my($tP1,$tP2,@tOmitPos)=@_;
my@tPA1=split(//,$tP1);
my@tPA2=split(//,$tP2);
my$tScore=0;
my%tOP=();
for(my$tI=0;$tI<@tOmitPos;$tI++){$tOP{"$tOmitPos[$tI]"}=1;}for(my$tI=0;$tI<@tPA1;$tI++){if(!(defined$tOP{"$tI"})){for(my$tJ=0;$tJ<@tPA2;$tJ++){if(!(defined$tOP{"$tI"})){my$tS=$egCharPairScore{"$tPA1[$tI];$tPA2[$tJ]"};
if(defined$tS){}else{my$tS1=$egCharToAnyScore{"$tPA1[$tI]"};
my$tS2=$egCharToAnyScore{"$tPA2[$tI]"};
if(defined$tS1){$tS=$tS1;}elsif(defined$tS2){$tS=$tS2;}else{$tS=0;}}$tScore=$tScore+$tS;}}}}return$tScore;}sub efStatList{my(@sList)=@_;
my$sAvg=0;
my$sStd=0;
my$si=0;
my@ssList=@sList;
@sList=sort{$a<=>$b}@ssList;
for($si=0;$si<@sList;$si++){$sAvg=$sAvg+$sList[$si];}if($si==0){$sAvg=0;}else{$sAvg=$sAvg/$si;}for($si=0;$si<@sList;$si++){$sStd=$sStd+($sList[$si]-$sAvg)*($sList[$si]-$sAvg);}if($si==0){$sStd=0;}else{$sStd=sqrt($sStd/$si);}my$sMin=$sList[0];
my$sMax=$sList[@sList-1];
return($sMin,$sMax,$sAvg,$sStd);}sub efFormatSeq{my($fSeq,$tBlock,$tNumber)=@_;
my$tPos=0;
my$tRSeq="";
while(1){my$tStart=$tBlock*$tPos;
my$tEnd=$tBlock*($tPos+1)-1;
if($tStart>=length($fSeq)){last;}if($tEnd>=length($fSeq)){$tEnd=length($fSeq)-1;}$tRSeq=$tRSeq.substr($fSeq,$tStart,abs($tEnd-$tStart)+1)." ";
if(($tPos>0)and($tPos%$tNumber==0)){$tRSeq=$tRSeq."\n     ";}$tPos++;}$tRSeq=~s/^ +//g;
$tRSeq=~s/ +$//g;
$tRSeq="     ".$tRSeq;
return$tRSeq;}sub efCalculateDistance{my($lMode)=@_;
my$lRV=$egVoidValue;
if($lMode==1){$lRV=&efEuclideanDistance();}elsif($lMode==2){$lRV=&efCosinDistance();}else{$lRV=$egVoidValue;}return$lRV;}sub efEuclideanDistance{my$eRV=0;
my$eMaxNumber=scalar@lVA;
if($eMaxNumber>scalar@lVB){$eMaxNumber=scalar@lVB;}for(my$ei=0;$ei<$eMaxNumber;$ei++){my$eDiff=($lVA[$ei]-$lVB[$ei]);
$eRV=$eRV+$eDiff*$eDiff;}return sqrt($eRV);}sub efCosinDistance{my$eRV=0;my$ei=0;
my$eAvgA=0;my$eAvgB=0;
my$eMaxNumber=scalar@lVA;
if($eMaxNumber>scalar@lVB){$eMaxNumber=scalar@lVB;}for($ei=0;$ei<$eMaxNumber;$ei++){$eAvgA=$eAvgA+$lVA[$ei];}if($ei==0){$eAvgA=0;}else{$eAvgA=$eAvgA/$ei;}for($ei=0;$ei<$eMaxNumber;$ei++){$eAvgB=$eAvgB+$lVB[$ei];}if($ei==0){$eAvgB=0;}else{$eAvgB=$eAvgB/$ei;}my$eFengZi=0;
my$eFengMuA=0;
my$eFengMuB=0;
for($ei=0;$ei<$eMaxNumber;$ei++){$eFengZi=$eFengZi+($lVA[$ei]-$eAvgA)*($lVB[$ei]-$eAvgB);
$eFengMuA=$eFengMuA+($lVA[$ei]-$eAvgA)*($lVA[$ei]-$eAvgA);
$eFengMuB=$eFengMuB+($lVB[$ei]-$eAvgB)*($lVB[$ei]-$eAvgB);}if(($eFengMuA<=0)or($eFengMuB<=0)){$eRV=0;}else{$eRV=$eFengZi/sqrt($eFengMuA*$eFengMuB);}return abs($eRV);}sub efBinList{my($bBinNum,@bDataA)=@_;
my$bMax=$bDataA[0];my$bMin=$bDataA[0];my$bSpacer=0;
my$bi=0;
for($bi=0;$bi<@bDataA;$bi++){if($bDataA[$bi]>$bMax){$bMax=$bDataA[$bi];}if($bDataA[$bi]<$bMin){$bMin=$bDataA[$bi];}}$bSpacer=(abs($bMax-$bMin))/$bBinNum;
my@bValue=();
for($bi=0;$bi<$bBinNum;$bi++){$bValue[$bi]=0;}for($bi=0;$bi<@bDataA;$bi++){my$bid=int(($bDataA[$bi]-$bMin)/$bSpacer);
if($bid<0){$bid=0;}if($bid>=$bBinNum){$bid=$bBinNum-1;}$bValue[$bid]++;}return($bMin,$bMax,$bSpacer,@bValue);}sub efRemoveHTMLTags{my($rLine)=@_;
$rLine=~s/<[a-zA-Z\s\d\#\"\:\-\=\/\.\_\?\%\&\,\'\"\?\+\[\]]+>//g;
$rLine=~s/\&nbsp;/ /g;
return$rLine;}sub efWGet{my($tURL,$tFile)=@_;
my$tCmdLine="wget --output-document=\"$tFile\" \"$tURL\" -q";
if(system("$tCmdLine")){die"-Error when downloading the web page!\n";}}sub efGetTaxByID{my($gTaxID)=@_;
my@gTax=();my$gNum=0;
%egTaxRank=();
while(1){my$gid=$egTaxID2ID{"$gTaxID"};
if((defined$gid)and($egTaxID[$gid]!=$egPTaxID[$gid])and($egTaxID[$gid]!=1)and(defined$egTaxLevelIdx{"$egRank[$gid]"})){$gTax[$gNum]="$egRank[$gid]\=\>$egName[$gid]";
$gTaxID=$egPTaxID[$gid];
$egTaxRank{"$egRank[$gid]"}=$egName[$gid];
$gNum++;}else{if(defined$gid){$gTaxID=$egPTaxID[$gid];}if((!(defined$gid))or(($egTaxID[$gid]==$egPTaxID[$gid])and($egTaxID[$gid]==1))){last;}}}return@gTax;}sub efLoadTaxDB{my$egFileNodes="nodes.dmp";
my$egFileNames="names.dmp";
open(efIn,"$egFileNodes")or die"Error when loading the NodeFile!\n";
while(<efIn>){$tline=$_;$tline=~s/[\r\n]//g;
@tlist=split(/\t\|\t/,$tline);
if($_=~/^#/){}elsif((@tlist>=13)and($tlist[12]=~/\t\|$/)){$egTaxID2ID{"$tlist[0]"}=$egNum;
$egPTaxID[$egNum]=$tlist[1];
$egTaxID[$egNum]=$tlist[0];
$egRank[$egNum]=$tlist[2];
$egName[$egNum]="";
$egNum++;}}close(efIn);
open(efIn,"$egFileNames")or die"Error when loading the NameFile!\n";
while(<efIn>){$tline=$_;$tline=~s/[\r\n]//g;
@tlist=split(/\t\|\t/,$tline);
if(defined$tlist[3]){$tlist[3]=~s/\t\|$//g;}if($tline=~/^#/){}elsif(@tlist>=4){$tid=$egTaxID2ID{"$tlist[0]"};
if((defined$tid)and($tlist[3]eq"scientific name")and($egName[$tid]eq"")){$egName[$tid]=$tlist[1];
$egName2ID{"$tlist[1]"}=$tid;}}}close(efIn);}sub efGetMITE{my($seq,$s,$e,$strand,$tir,$dr)=@_;
my($tir5,$tir3,$dr5,$dr3,$mseq)=("","","","","");
if($strand eq"+"){$tir5=substr($seq,$s,$tir);
$dr5=substr($seq,$s-$dr,$dr);
$tir3=substr($seq,$e-$tir+1,$tir);
$dr3=substr($seq,$e+1,$dr);
$mseq=substr($seq,$s,abs($e-$s)+1);}else{$tir3=&efDNA_RC(substr($seq,$s,$tir));
$dr3=&efDNA_RC(substr($seq,$s-$dr,$dr));
$tir5=&efDNA_RC(substr($seq,$e-$tir+1,$tir));
$dr5=&efDNA_RC(substr($seq,$e+1,$dr));
$mseq=&efDNA_RC(substr($seq,$s,abs($e-$s)+1));}return($tir5,$tir3,$dr5,$dr3,$mseq);}sub efDNA_RC{my($qLine)=@_;
my@qList=split(//,$qLine);
my$qi=0;
my@qL=();
my$qLength=scalar@qList;
if($qLength>=1){for($qi=$qLength-1;$qi>=0;$qi--){if($qList[$qi]eq"A"){$qL[$qLength-$qi-1]="T";}elsif($qList[$qi]eq"T"){$qL[$qLength-$qi-1]="A";}elsif($qList[$qi]eq"G"){$qL[$qLength-$qi-1]="C";}elsif($qList[$qi]eq"C"){$qL[$qLength-$qi-1]="G";}else{$qL[$qLength-$qi-1]="N";}}$qLine=join("",@qL);}else{$qLine="";}return$qLine;}sub efIsSameStrand{my($iSeq1,$iSeq2)=@_;
if((length($iSeq1)>0)and(length($iSeq2)>0)){my$iF1="seq1.fasta";open(iO,">$iF1")or die"Error when saving sequence 1\n";print iO">Query\n$iSeq1\n";close(iO);
my$iF2="seq2.fasta";open(iO,">$iF2")or die"Error when saving sequence 2\n";print iO">DB\n$iSeq2\n";close(iO);
my$iFO="seq.bl2seq.dat";
my$iCmdLine="bl2seq -p blastn -i \"$iF1\" -j \"$iF2\" -o $iFO";
if(system($iCmdLine)){print"-Error when running bl2seq!\n";
return-1;}else{my$iSearchIO=new Bio::SearchIO(-file=>"$iFO",-format=>"blast");
my$iFlag=1;
while(my$iResult=$iSearchIO->next_result){my$iNameQ=$iResult->query_name();
while(my$iHit=$iResult->next_hit){my$iNameH=$iHit->name();
while(my$iHSP=$iHit->next_hsp){my$iStrandQ=$iHSP->strand('query');
my$iStrandH=$iHSP->strand('hit');
if(($iStrandQ eq$iStrandH)and($iStrandQ ne"")){$iFlag=1;}else{$iFlag=0;}goto lbGotIt;}}}lbGotIt:return$iFlag;}}else{return-1;}}sub efID2Seq{my($fInID)=@_;
my@fList=();
for(my$fi=0;$fi<$egK;$fi++){$fList[$egK-1-$fi]=0;
my$fj=$fInID%$egV;
$fList[$egK-1-$fi]=$fj;
$fInID=int(($fInID-$fj)/$egV);}my$fLine="";
for(my$fi=0;$fi<$egK;$fi++){my$fc=$egC[$fList[$fi]];
if(defined$fc){$fLine=$fLine.$fc;}else{die"Error! C Index: $fList[$fi]\n";}}return$fLine;}sub efSeq2ID{my(@flist)=@_;
my$fInID=0;
if(@flist==$egK){$fInID=0;
for(my$fi=0;$fi<$egK;$fi++){my$fid=$egCI{"$flist[$fi]"};
if(defined$fid){$fInID=($fInID+$fid)*$egV;}else{die"Error! \$fid => $flist[$fi]\n";}if($fi==$egK-1){$fInID=int($fInID/$egV);}}}else{die"Error! Input line: $egGSeq\n";}return$fInID;}sub efRunPiler{my($tFileIn)=@_;
if(system("pals -self \"$tFileIn\" -out \"$tFileIn.gff\" 2>/dev/null 1>/dev/null")){die"-Error when running pals!\n";}else{print" [pals]";}if(system("piler -trs \"$tFileIn.gff\" -out \"$tFileIn.trs.gff\" 2>/dev/null 1>/dev/null")){die"-Error when running piler!\n";}else{print" [piler]";}if((system("mkdir -p fams"))or(system("piler -trs2fasta \"$tFileIn.trs.gff\" -seq \"$tFileIn\" -path fams 2>/dev/null 1>/dev/null"))){die"-Error when running piler -trs2fasta!\n";}else{print" [piler-trs2fasta]";}my@tFileList=`ls fams`;
my$ti=0;
for($ti=0;$ti<@tFileList;$ti++){$tFileList[$ti]=~s/[\r\n]//g;}if(system("mkdir -p aligned_fams")){die"-Error when creating directory!\n";}for($ti=0;$ti<@tFileList;$ti++){if(system("muscle -in \"fams\/$tFileList[$ti]\" -out \"aligned_fams\/$tFileList[$ti]\" -maxiters 1 -diags 1 2>/dev/null 1>/dev/null")){die"-Error when running muscle for: fams\/$tFileList[$ti]\n";}}print" [muscle]";
if(system("mkdir -p cons")){die"-Error when creating directory!\n";}for($ti=0;$ti<@tFileList;$ti++){if(system("piler -cons \"fams\/$tFileList[$ti]\" -out \"cons\/$tFileList[$ti]\" -label \"$tFileList[$ti]\" 2>/dev/null 1>/dev/null")){die"-Error when running piler-cons\n";}}if(system("cat cons\/* > piler_library.fasta")or system("cp -r cons \"$egDirNow\/poplar_piler\"")or system("cp -r piler_library.fasta \"$egDirNow\/poplar_piler\"")){die"-Error!\n";}else{print" [DONE:piler]";}}sub efgdtrc{my($ta,$tb,$tc)=@_;
my$trv=0;
eval{$trv=gdtrc($ta,$tb,$tc);};
if($@){$trv=0;}else{}return$trv;}sub efIsOverlap{my($s1,$e1,$s2,$e2)=@_;
my$t=0;
if($s1>$e1){$t=$s1;$s1=$e1;$e1=$t;}if($s2>$e2){$t=$s2;$s2=$e2;$e2=$t;}if((($s2>=$s1)and($s2<=$e1))or(($s1>=$s2)and($s1<=$e2))){return 1;}else{return 0;}}
