#!/usr/bin/perl -w
use strict;
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%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==3 or die""."Error in syntax!\n"."    ./filterClusterLength.pl <MinSameLengthCopyNum> <input.MITE.txt> <output.MITE.txt>\n"."Note:\n"."    Only those clusters with at least <MinSameLengthCopyNum> copies of the same lengths\n"."will be kept for further analysis.\n";
my($egMinSameLengthCopyNum,$egFileI,$egFileO)=@ARGV;
my%egCluster2ID=();
my@egClusterID=();
my@egLength=();
my@egTIRLength=();
my@egDRLength=();
my$egNum=0;
print"Loading the file ... ";
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=~/^#/){}elsif(@tlist>=23){$tid=$egCluster2ID{"$tlist[2]"};
if(defined$tid){$egLength[$tid].=";$tlist[6]";
$egTIRLength[$tid].=";$tlist[11]";
$egDRLength[$tid].=";$tlist[7]";}else{$tid=$egNum;
$egCluster2ID{"$tlist[2]"}=$tid;
$egClusterID[$tid]=$tlist[2];
$egLength[$tid].="$tlist[6]";
$egTIRLength[$tid].="$tlist[11]";
$egDRLength[$tid].="$tlist[7]";
$egNum++;}}}close(efI);
print" [done] [Loading clusters:$egNum]\n";
print"Screening for candidate MITEs ... ";
%tidx=();
my@pL=();
my@pTIR=();
my@pDR=();
my$pN=0;
my@pID=();
my$pCluster=0;
my@pNewID=();
for($tempi=0;$tempi<$egNum;$tempi++){@pL=();@pTIR=();@pDR=();$pN=0;
@pID=();
$pCluster=$egClusterID[$tempi];
@pL=split(/;/,$egLength[$tempi]);
@pTIR=split(/;/,$egTIRLength[$tempi]);
@pDR=split(/;/,$egDRLength[$tempi]);
$pN=scalar@pL;
for($tempj=0;$tempj<$pN;$tempj++){$pID[$tempj]=$tempj;}@pNewID=sort{$pL[$a]<=>$pL[$b]}@pID;
my$pFlag=1;
for($tempj=0;$tempj<$pN-$egMinSameLengthCopyNum;$tempj++){$pFlag=1;
for($tempk=0;$tempk<$egMinSameLengthCopyNum;$tempk++){if(($pL[$pNewID[$tempj]]==$pL[$pNewID[$tempj+$tempk]])and($pTIR[$pNewID[$tempj]]==$pTIR[$pNewID[$tempj+$tempk]])and($pDR[$pNewID[$tempj]]==$pDR[$pNewID[$tempj+$tempk]])){;}else{$pFlag=0;
last;}}if($pFlag==1){$tidx{"$pCluster"}=1;
last;}}}print" [done] [Candidate MITEs:".(scalar keys%tidx)."]\n";
print"Screening and saving ... ";
open(efI,"$egFileI")or die"Q!\n";
open(efO,">$egFileO")or die"E!\n";
$tempi=0;
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){print efO"$tline\n";}elsif(@tlist>=23){if(defined$tidx{"$tlist[2]"}){print efO"$tline\n";
$tempi++;}}}close(efO);
close(efI);
print" [done] [Saved MITE copies:$tempi]\n";
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;}}
