#!/usr/bin/perl -w
use strict;
use warnings;
use MIME::Base64;
use Bio::SeqIO;
use Bio::SearchIO;
use threads;
use threads::shared;
my$tline="";my@tlist=();my$tempi=0;my$tempj=0;my$tempk=0;my$pline="";my@plist=();my$tid=0;my%tidx=();
my$egDirTemp="tempData";my$egScriptMUST="MUST";my$egScriptPWD="PairWiseDistance";my$egScriptFilter="filterRedundancy.new.pl";my$egLogFilter="log.filterRedundancy.new.log";
if(system("echo \"\" > \"$egLogFilter\"")){die"-Error when clearing the log file!\n";}my$egDirScript=$0;
@tlist=split(/\//,$egDirScript);
if(@tlist==1){$egDirScript=".";}else{$egDirScript=~s/\/$tlist[@tlist-1]$//g;}$|=1;
my$gMinTIR=8;
my$gMaxTIR=50;
my$gMinDR=2;
my$gMaxDR=30;
my$gMinMITE=100;
my$gMaxMITE=600;
my$gFixedFlanking=50;
my$gMutationRate=0.80;
my$egFileSeq="";
my$egFileOut="";
my$egPTNum=10;
if((@ARGV==3)or(@ARGV==12)){}else{print"Error in syntax!\n";
&efSyntax();
die"-----------------------------------------------------------------------------\n";}$egFileSeq=$ARGV[0];$egFileOut=$ARGV[1];$egDirTemp=$ARGV[2];
my$egFileSize=-s$egFileSeq;
my$egSizeOne=int($egFileSize/$egPTNum+1.0);
if(-d$egDirTemp){}else{if(system("mkdir -p \"$egDirTemp\"")){die"-Error when creating the temporary directory!\n";}}if(@ARGV==12){$egPTNum=$ARGV[3];$gMinTIR=$ARGV[4];$gMaxTIR=$ARGV[5];$gMinDR=$ARGV[6];$gMaxDR=$ARGV[7];$gMinMITE=$ARGV[8];$gMaxMITE=$ARGV[9];$gFixedFlanking=$ARGV[10];$gMutationRate=$ARGV[11];}my$egCutOffIdent=$gMutationRate;
my%mMITE2ID=();
my@mID=();
my@mGID=();
my@mStart=();
my@mEnd=();
my@mStrand=();
my@mLength=();
my@mDR=();
my@mDRIdent=();
my@mDRLeft=();
my@mDRRight=();
my@mTIR=();
my@mTIRIdent=();
my@mTIRLeft=();
my@mTIRRight=();
my@mDRAT=();
my@mTIRAT=();
my@mMITEAT=();
my@mFDRAT=();
my@mPreSite=();
my@mPostSite=();
my@mMITE=();
my@mCopyNum=();
my@mCopies=();
my@mTotalScore=();
my@mFlag=();
my$mNum=0;
my%mData=();
my%egExistCopy=();
my$egCutOffCopyScore=10;
my$egLineLength=60;
my$egTimeOld=`date +%s`;
my$egTimeNew=`date +%s`;
my$egSeq="";
my$egID="";
my$lFile="";
my$lFSeq="";
my$lCmdLine="";
my$lFileDat="";
my$lFileTIR="";
my$lFileEdge="";
$egTimeOld=`date +%s`;
my$egTID=1;
my$egFileSub="$egDirTemp/sub-genome.$egTID.fasta";
my$egSizeNow=$egSizeOne;
print"Spliting the genome into $egPTNum sub-genomes ... ";
open(efO,">$egFileSub")or die"Error:$egTID!\n";
my$egSeqObj=new Bio::SeqIO(-file=>"$egFileSeq",-format=>"fasta");
while(my$egSeqO=$egSeqObj->next_seq){$tline=$egSeqO->seq;$egSeq="\U$tline\E";$egID=$egSeqO->id;
my$tStartPos=0;
my$tLength=length($egSeq);
while(1){if($egTID>=$egPTNum){print efO">$tStartPos;;$egID\n".substr($egSeq,$tStartPos,$tLength-$tStartPos)."\n";
$egSizeNow=0;$tStartPos=0;
last;}elsif($egSizeNow>=$tLength-$tStartPos){print efO">$tStartPos;;$egID\n".substr($egSeq,$tStartPos,$tLength-$tStartPos)."\n";
$egSizeNow-=$tLength-$tStartPos;
$tStartPos=0;
last;}else{print efO">$tStartPos;;$egID\n".substr($egSeq,$tStartPos,$egSizeNow)."\n";
$tStartPos+=$egSizeNow-2*($gFixedFlanking+$gMaxDR+$gMaxMITE+$gMaxDR+$gFixedFlanking);
close(efO);
$egTID++;
open(efO,">$egDirTemp/sub-genome.$egTID.fasta")or die"EEE:$egTID!\n";
$egSizeNow=$egSizeOne;}}}$egTimeNew=`date +%s`;
close(efO);
print" [SubFiles:$egTID] [done] [Elapsed time: ".($egTimeNew-$egTimeOld)." seconds]\n";
print"Scanning the nucleotide sequences for potential MITEs ... ";
$egTimeOld=`date +%s`;
my@egPTID=();
my$egSleepTime=5;
for($tempi=1;$tempi<=$egPTNum;$tempi++){my$tPID=threads->new(\&efPTMUST,$tempi,"$egDirTemp/sub-genome.$tempi.fasta");
$egPTID[$tempi]=$tPID;}$egTimeNew=`date +%s`;
%tidx=();
while(1){for($tempi=1;$tempi<=$egPTNum;$tempi++){$tid=$egPTID[$tempi];
if(($tid->is_joinable())and(!(defined$tidx{"$tempi"}))){$tid->join();$tidx{"$tempi"}=1;}}if((scalar keys%tidx)==$egPTNum){last;}sleep($egSleepTime);}$egTimeNew=`date +%s`;
print"Total calculation time elapsed: ".($egTimeNew-$egTimeOld)." seconds.\n";
print"Loading the annotation data ... ";
$egTimeOld=`date +%s`;
$tempj=0;$tempk=0;
for($tempi=1;$tempi<=$egPTNum;$tempi++){$tempj=&efLoadMUST("$egDirTemp/$tempi.temp_seq.dat");
$tempk+=$tempj;}$egTimeNew=`date +%s`;
print" [RawMITEs:$tempk] [done] [Time:".($egTimeNew-$egTimeOld)." seconds]\n";
my$egCmdLine="";
print"Removing redundancy in the predicted MITEs ... ";
my@mSortID=();
my@mNewID=();
for($tempi=0;$tempi<$mNum;$tempi++){$mSortID[$tempi]=$tempi;}@mNewID=sort{$mGID[$a]cmp$mGID[$b]or$mStart[$a]<=>$mStart[$b]or$mEnd[$a]<=>$mEnd[$b]}@mSortID;
$tempk=0;
for($tempi=0;$tempi<@mNewID-1;$tempi++){my($newA,$newB)=($mNewID[$tempi],$mNewID[$tempi+1]);
$tempj=1;
while(($mFlag[$newA]==1)and($mEnd[$newA]>=$mStart[$newB])){if($mFlag[$newB]==0){$tempj++;$newB=$mNewID[$tempi+$tempj];}elsif($mTotalScore[$newA]>=$mTotalScore[$newB]){$mFlag[$newB]=0;
$tempk++;
$tempj++;
$newB=$mNewID[$tempi+$tempj];}else{$mFlag[$newA]=0;
$tempk++;}}}print" [RuleBasedRemoving:$tempk]";
$tempj=0;
for($tempi=0;$tempi<$mNum;$tempi++){if($mFlag[$tempi]==1){$tempj++;}}$egTimeNew=`date +%s`;
print" [MITE:$tempj] [done] [Time:".($egTimeNew-$egTimeOld)." seconds]\n";
$egTimeOld=$egTimeNew;
my$egCutOffCopyNum=3;
my$egCutOffEValue=1e-20;
print"Clustering ... ";
my@mClusterID=();
my%mC2IDs=();
for($tempi=0;$tempi<$mNum;$tempi++){$mClusterID[$tempi]=-1;}my$mCluster=0;
my$mMITEWithCluster=0;
my%tInCluster=();
my$tFirstMITE=-1;
my@tInClusterList=();
my@tNeighborList=();
my$egClusterMode=1;
my$egFileMITESeq="$egDirTemp/temp-mite-seq.fasta";
open(efO,">$egFileMITESeq")or die"E!\n";
for($tempi=0;$tempi<$mNum;$tempi++){print efO">$tempi\n$mMITE[$tempi]\n";}close(efO);
system("formatdb -p F -i \"$egFileMITESeq\"");
system("megablast -d \"$egFileMITESeq\" -i \"$egFileMITESeq\" -F F -e $egCutOffEValue -m 9 -o \"$egDirTemp/temp-mite-seq.all-vs-all.megablast.txt\" 2>&1 1>>\"$egDirTemp/temp-mite-seq.all-vs-all.megablast.log\" ");
print" [BLAST]";
my%idxIdent=();
my%idxSameStrand=();
my@mCC=();
my@mFF=();
open(efI,"$egDirTemp/temp-mite-seq.all-vs-all.megablast.txt")or die"EE!\n---[$egDirTemp/temp-mite-seq.all-vs-all.megablast.txt]---\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){}elsif(@tlist>=12){if($tlist[0]<$tlist[1]){$idxIdent{$tlist[0].$tlist[1]}=$tlist[2]*$tlist[3]/100;
if($tlist[8]<=$tlist[9]){$idxSameStrand{$tlist[0].$tlist[1]}=1;}}elsif($tlist[0]>$tlist[1]){$idxIdent{$tlist[1].$tlist[0]}=$tlist[2]*$tlist[3]/100;
if($tlist[8]<=$tlist[9]){$idxSameStrand{$tlist[1].$tlist[0]}=1;}}}}close(efI);
for($tempi=0;$tempi<$mNum;$tempi++){$mCC[$tempi]="$tempi";
$mFF[$tempi]=0;
$mClusterID[$tempi]=0;}for($tempi=0;$tempi<$mNum;$tempi++){for($tempj=$tempi+1;$tempj<$mNum;$tempj++){$tid=$idxIdent{$tempi.$tempj};
if(defined$tid){if(($tid>=$egCutOffIdent*$mLength[$tempi])and($tid>=$egCutOffIdent*$mLength[$tempj])){$mCC[$tempi].=" $tempj";
$mCC[$tempj].=" $tempi";}}}}my@mCFlag=();
my$egClusterNumber=1;
for($tempi=0;$tempi<$mNum;$tempi++){my$kCID=$egClusterNumber;
$mCC[$tempi]=~s/^ +//g;
@tlist=split(/ /,$mCC[$tempi]);
if((scalar@tlist)<$egCutOffCopyNum){$mCFlag[$tempi]=0;}else{$mCFlag[$tempi]=1;
for($tempj=0;$tempj<@tlist;$tempj++){if(($mClusterID[$tlist[$tempj]]>0)and($mClusterID[$tlist[$tempj]]<$kCID)){$kCID=$mClusterID[$tempj];}}for($tempj=0;$tempj<@tlist;$tempj++){if($mClusterID[$tlist[$tempj]]==0){$mClusterID[$tlist[$tempj]]=$kCID;}}if($mClusterID[$tempi]==0){$mClusterID[$tempi]=$kCID;}if($kCID==$egClusterNumber){$egClusterNumber++;}}}$mCluster=$egClusterNumber-1;
for($tempi=0;$tempi<$mNum;$tempi++){if($mCFlag[$tempi]==1){$tid=$mC2IDs{"$mClusterID[$tempi]"};
if(defined$tid){$tid.=";$tempi";}else{$tid="$tempi";}$mC2IDs{"$mClusterID[$tempi]"}=$tid;}}$egTimeNew=`date +%s`;
print" [done] [ClusterNum:$mCluster] [Time:".($egTimeNew-$egTimeOld)." seconds]\n";
$egTimeOld=$egTimeNew;
print"Deciding the strand information of each valid MITE ... ";
for($tempi=0;$tempi<$mNum;$tempi++){if(($mFlag[$tempi]==1)and($mCFlag[$tempi]==1)and($mStrand[$tempi]eq"")){$mStrand[$tempi]="+";
for($tempj=$tempi+1;$tempj<$mNum;$tempj++){if(($mFlag[$tempj]==1)and($mCFlag[$tempj]==1)and($mClusterID[$tempj]==$mClusterID[$tempi])){if(defined$idxSameStrand{$tempi.$tempj}){$mStrand[$tempj]="+";}else{$mStrand[$tempj]="-";}}}}}$egTimeNew=`date +%s`;
print" [done] [Time:".($egTimeNew-$egTimeOld)." seconds]\n";
$egTimeOld=$egTimeNew;
print"Saving the data ... ";
$tempj=0;
open(efOut,">$egFileOut\.before_filter")or die"Error when saving the predicted MITEs!\n";
print efOut"#GID	ID	Cluster	Start	End	Strand	Length	DR	DRIdent	DRLeft	DRRight	TIR	TIRIdent	TIRLeft	TIRRight	DRAT	TIRAT	MITEAT	FDRAT	PreSite	PostSite	MITE	TotalScore\n";
for($tempi=0;$tempi<$mNum;$tempi++){if(($mCFlag[$tempi]==1)and($mClusterID[$tempi]>0)and($mFlag[$tempi]==1)){$mCopies[$tempi]=~s/\s/,/g;
print efOut"$mGID[$tempi]	$mID[$tempi]	$mClusterID[$tempi]	$mStart[$tempi]	$mEnd[$tempi]	$mStrand[$tempi]	$mLength[$tempi]	$mDR[$tempi]	$mDRIdent[$tempi]	$mDRLeft[$tempi]	$mDRRight[$tempi]	$mTIR[$tempi]	$mTIRIdent[$tempi]	$mTIRLeft[$tempi]	$mTIRRight[$tempi]	$mDRAT[$tempi]	$mTIRAT[$tempi]	$mMITEAT[$tempi]	$mFDRAT[$tempi]	$mPreSite[$tempi]	$mPostSite[$tempi]	$mMITE[$tempi]	$mTotalScore[$tempi]\n";
$tempj++;}}close(efOut);
if(system("$egDirScript/filterCopyNum.pl \"$egFileOut.before_filter\" \"$egFileOut.filter-1\" $egCutOffCopyNum >/dev/null")or system("$egDirScript/filterSegDup-2.pl \"$egFileOut.filter-1\" \"$egFileOut\" 4 3 >/dev/null")){die"-Error when filtering!\n";}$egTimeNew=`date +%s`;
print" [done] [ValidMITENum:$tempj] [Time:".($egTimeNew-$egTimeOld)." seconds]\n";
$egTimeOld=$egTimeNew;
print"Screening for all the MITE copies ... ";
open(efI,"$egFileOut")or die"Q!\n";
open(efO,">$egFileOut.before_expansion")or die"!!\n";
$tempi=0;
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){print efO"$tline\n";}elsif(@tlist>=23){$tlist[1]=$tempi;
$tlist[22]=&efEScore($tlist[7],$tlist[11],$tlist[8],$tlist[12]);
print efO"".join("\t",@tlist)."\n";
$tempi++;}}close(efO);
close(efI);
print" [BackupMITEs:$tempi]";
my$tFileTemp="$egDirTemp/temp-templates.fasta";
my@tGID=();
my@tCluster=();
my@tMITEID=();
my@tS=();
my@tE=();
my@tStrand=();
my@tDR=();
my@tTIR=();
my@tDRIdent=();
my@tTIRIdent=();
my@tEScore=();
my@tMITESeq=();
my$tN=0;
my$tMITECopy=0;
$tempk=0;
open(efI,"$egFileOut.before_expansion")or die"E:TemplateGeneration!\n";
open(efO,">$tFileTemp")or die"E:SaveTemplate!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){;}elsif(@tlist>=23){print efO">MITE;;$tlist[1];;$tlist[2];;$tlist[7];;$tlist[8];;$tlist[11];;$tlist[12]\n$tlist[21]\n";
$tGID[$tN]=$tlist[0];$tMITEID[$tN]=$tlist[1];$tCluster[$tN]=$tlist[2];$tS[$tN]=$tlist[3];$tE[$tN]=$tlist[4];$tStrand[$tN]=$tlist[5];
$tDR[$tN]=$tlist[7];$tDRIdent[$tN]=$tlist[8];$tTIR[$tN]=$tlist[11];$tTIRIdent[$tN]=$tlist[12];
$tEScore[$tN]=&efEScore($tDR[$tN],$tTIR[$tN],$tDRIdent[$tN],$tTIRIdent[$tN]);$tMITESeq[$tN]=$tlist[21];
$tempk++;$tN++;}}close(efO);
close(efI);
print" [Templates:$tempk]";
$tMITECopy=$tN;
my$tFileSeq="$egDirTemp/temp-genome-seq.fasta";
my$tFileBLAT="$egDirTemp/temp-genome-seq.blat-psl.txt";
if(system("cp \"$egFileSeq\" \"$tFileSeq\"")or system("formatdb -p F -i \"$tFileSeq\"")or system("blat \"$tFileSeq\" \"$tFileTemp\" -minIdentity=95 -minScore=100 \"$tFileBLAT\" 2>&1 1>>\"$egDirTemp/temp-genome-seq.blat.log\" ")){die"-Error when mapping!\n";}print" [BLAT]";
my@pGID=();
my@pS=();
my@pE=();
my@pStrand=();
my@pMID=();
my@pMCluster=();
my@pMDR=();
my@pMDRIdent=();
my@pMTIR=();
my@pMTIRIdent=();
my@pMBitScore=();
my$pN=0;
open(efI,"$tFileBLAT")or die"Error!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){}elsif(($tline=~/^\d+/)and(@tlist>=12)){my$tstrand="+";
$tstrand=$tlist[8];
my($tmid,$tcluster,$tmdr,$tmdrident,$ttir,$ttirident)=(-1,-1,0,0,0,0);
if($tlist[9]=~/^MITE;;(\d+);;(\d+);;(\d+);;([\d\.]+);;(\d+);;([\d\.]+)/){($tmid,$tcluster,$tmdr,$tmdrident,$ttir,$ttirident)=($1,$2,$3,$4,$5,$6);}$tid=-1;
for($tempi=0;$tempi<$tN;$tempi++){if(($tGID[$tempi]eq$tlist[13])and(&efIsOverlap($tS[$tempi],$tE[$tempi],$tlist[15],$tlist[16])==1)){$tid=-2;last;}}if($tid==-1){for($tempi=0;$tempi<$pN;$tempi++){if(($pGID[$tempi]eq$tlist[13])and(&efIsOverlap($pS[$tempi],$pE[$tempi],$tlist[15],$tlist[16])==1)){$tid=$tempi;last;}}if($tid!=-1){if($pMBitScore[$tid]<$tlist[0]){$pGID[$tid]=$tlist[13];$pS[$tid]=$tlist[15];$pE[$tid]=$tlist[16];$pStrand[$tid]=$tstrand;
$pMID[$tid]=$tmid;$pMCluster[$tid]=$tcluster;$pMDR[$tid]=$tmdr;$pMDRIdent[$tid]=$tmdrident;$pMTIR[$tid]=$ttir;$pMTIRIdent[$tid]=$ttirident;
$pMBitScore[$tid]=$tlist[0];}}else{$tid=$pN;
$pGID[$tid]=$tlist[13];$pS[$tid]=$tlist[15];$pE[$tid]=$tlist[16];$pStrand[$tid]=$tstrand;
$pMID[$tid]=$tmid;$pMCluster[$tid]=$tcluster;$pMDR[$tid]=$tmdr;$pMDRIdent[$tid]=$tmdrident;$pMTIR[$tid]=$ttir;$pMTIRIdent[$tid]=$ttirident;
$pMBitScore[$tid]=$tlist[0];
$pN++;}}}}close(efI);
print" [OtherCopy:$pN]";
my@pSeqL=();
my@pSeqM=();
my@pSeqR=();
my%pGID2ID=();
if($pN==0){goto labelNoExpanding;}for($tempi=0;$tempi<$pN;$tempi++){$tid=$pGID2ID{"$pGID[$tempi]"};
if(defined$tid){$tid.=";$tempi";}else{$tid="$tempi";}$pGID2ID{"$pGID[$tempi]"}=$tid;}$egSeqObj=new Bio::SeqIO(-file=>"$egFileSeq",-format=>"fasta");
while(my$egSeqO=$egSeqObj->next_seq){my($ppID,$ppSeq)=($egSeqO->id,$egSeqO->seq);
$tline=$pGID2ID{"$ppID"};
if(defined$tline){@tlist=split(/;/,$tline);
for($tempj=0;$tempj<@tlist;$tempj++){$tid=$tlist[$tempj];
$pSeqL[$tid]=substr($ppSeq,$pS[$tid]-1-$gFixedFlanking,$gFixedFlanking);
$pSeqM[$tid]=substr($ppSeq,$pS[$tid],abs($pE[$tid]-$pS[$tid])+1);
if($pStrand[$tid]ne"+"){$pSeqM[$tid]=&efDNA_RC($pSeqM[$tid]);}$pSeqR[$tid]=substr($ppSeq,$pE[$tid]-1+1,$gFixedFlanking);}}}if(system("cp \"$egFileOut.before_expansion\" \"$egFileOut.be-filter1\"")){die"-Error!\n";}open(efO,">>$egFileOut.be-filter1")or die"E:Expansion!\n";
$tempk=$tN;
for($tempi=0;$tempi<$pN;$tempi++){$tempk++;
print efO"$pGID[$tempi]	$tempk	$pMCluster[$tempi]	$pS[$tempi]	$pE[$tempi]	$pStrand[$tempi]	".(abs($pE[$tempi]-$pS[$tempi])+1)."	$mDR[$tempi]	0			$pMTIR[$tempi]	0			0	0	0	0	".($pSeqL[$tempi].$pSeqR[$tempi])."	".($pSeqL[$tempi].$pSeqR[$tempi])."	".$pSeqM[$tempi]."	0\n";}close(efO);
print" [AllCopySaved:$tempk]";
if(system("\"$egDirScript/refineBoundary.pl\" \"$egDirTemp\" $gMutationRate \"$tFileSeq\" \"$egFileOut.be-filter1\" \"$egFileOut.be-filter2\" 2>/dev/null 1>/dev/null")){print"\"$egDirScript/refineBoundary.pl\" \"$egDirTemp\" $gMutationRate \"$tFileSeq\" \"$egFileOut.be-filter1\" \"$egFileOut.be-filter2\"\n";
die"-Error when refineBoundary!\n";}print" [Boundary Refined]";
labelNoExpanding:
my($tF,$tP)=(0,0);
if($pN==0){open(efI,"$egFileOut.before_expansion")or die"E!\n";}else{open(efI,"$egFileOut.be-filter2")or die"E!\n";}open(efO,">$egFileOut")or die"Q!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){print efO"$tline	Copy status\n";}elsif(@tlist>=23){$tlist[22]=&efEScore($tlist[7],$tlist[11],$tlist[8],$tlist[12]);
$tlist[23]="";
if(($tlist[8]>=$gMutationRate)and($tlist[12]>=$gMutationRate)){$tlist[23]="Full Copy";$tF++;}else{$tlist[23]="Partial Copy";$tlist[7]=0;$tlist[8]=0;$tlist[9]="";$tlist[10]="";$tlist[11]=0;$tlist[12]=0;$tlist[13]="";$tlist[14]="";$tlist[15]=0;$tlist[16]=0;$tP++;}print efO"".join("\t",@tlist)."\n";}}close(efO);
close(efI);
print" [Full:$tF] [Partial:$tP]";
print" [done]\n";
if(system("rm -rf \"$egDirTemp\"")or system("rm -rf \"$egFileOut.before_filter\"")or system("rm -rf \"$egFileOut.filter-1\"")or system("rm -rf *.log")or system("rm -rf \"$egFileOut.before_expansion\"")or system("rm -rf \"$egFileOut.be-filter1\"")or system("rm -rf \"$egFileOut.be-filter2\"")){print"---Error when removing the temporary files!\n";}my$qMITENum=0;
my$qClusterNum=0;
%tidx=();
open(efI,"$egFileOut")or die"E!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\t/,$tline);
if($tline=~/^#/){}elsif(@tlist>=5){$qMITENum++;
$tidx{"$tlist[2]"}=1;}}$qClusterNum=scalar keys%tidx;
close(efI);
print"There are $qMITENum MITEs detected in $qClusterNum clusters.\n";
$|=0;
sub efEScore{my($adDR,$adTIR,$adIdentDR,$adIdentTIR)=@_;
my$adWeightIdent=20;
my$adWeightDR=2;
my$adEScore=$adIdentDR*$adDR*$adWeightDR+$adIdentTIR*$adTIR;
return$adEScore;}sub efLoadMUST{my($tFileDat)=@_;
%mData=();
my$tStartID=$mNum;
my$tNowNum=0;
open(tIn,"$tFileDat")or die"Error when loading the Dat!\n";
while(<tIn>){$tline=$_;
$tline=~s/[\r\n]//g;
my($ppp,$qqq)=("","");
if($tline=~/^\/\/$/){$egID=$mData{"GID"};
$mData{"SEQID"}=$egID;
my@ttlist=();
my$tNowID=&efValidString($mData{"ID"});
if($tNowID=~/^\d+$/){$mMITE2ID{"$egID $tNowID"}=$mNum;
$mID[$mNum]=&efValidString($tNowID);
$mGID[$mNum]=&efValidString($egID);
($ppp,$qqq)=split(/;;/,$mGID[$mNum],2);
if(!(defined$qqq)){$qqq=$ppp;$ppp=0;}$mGID[$mNum]=$qqq;
@ttlist=split(/-/,&efValidString($mData{"POSITION"}),2);
$mStart[$mNum]=$ppp+&efValidString($ttlist[0]);
$mEnd[$mNum]=$ppp+&efValidString($ttlist[1]);
if(!(defined$egExistCopy{"$mGID[$mNum]|$mStart[$mNum]-$mEnd[$mNum]"})){$mStrand[$mNum]="";
$mLength[$mNum]=abs($mEnd[$mNum]-$mStart[$mNum])+1;
$mDR[$mNum]=&efValidString($mData{"DR"});
$mDRIdent[$mNum]=&efValidString($mData{"DRIDENT"});
$mDRLeft[$mNum]=&efValidString($mData{"DRLEFT"});
$mDRRight[$mNum]=&efValidString($mData{"DRRIGHT"});
$mTIR[$mNum]=&efValidString($mData{"TIR"});
$mTIRIdent[$mNum]=&efValidString($mData{"TIRIDENT"});
$mTIRLeft[$mNum]=&efValidString($mData{"TIRLEFT"});
$mTIRRight[$mNum]=&efValidString($mData{"TIRRIGHT"});
$mDRAT[$mNum]=&efValidString($mData{"DRAT"});
$mTIRAT[$mNum]=&efValidString($mData{"TIRAT"});
$mMITEAT[$mNum]=&efValidString($mData{"MITEAT"});
$mFDRAT[$mNum]=&efValidString($mData{"FDRAT"});
$mPreSite[$mNum]=&efValidString($mData{"PRESITE"});
$mPostSite[$mNum]=&efValidString($mData{"POSTSITE"});
$mMITE[$mNum]=&efValidString($mData{"MITE"});
$mTotalScore[$mNum]=&efEScore($mDR[$mNum],$mTIR[$mNum],$mDRIdent[$mNum],$mTIRIdent[$mNum]);
$mCopyNum[$mNum]=1;
$mCopies[$mNum]="$mNum";
$mFlag[$mNum]=1;
$egExistCopy{"$mGID[$mNum]|$mStart[$mNum]-$mEnd[$mNum]"}=1;
$tNowNum++;
$mNum++;}}}elsif($tline=~/^([a-zA-Z]+)\t(.*)$/){my$tID=&efValidString($1);
my$tDat=&efValidString($2);
if(($tID ne"")and($tDat ne"")){$mData{"$tID"}=$tDat;}}}close(tIn);
return$tNowNum;}sub efPTMUST{my($tmPID,$tmFileSeq)=@_;
my$tmgTimeOld=`date +%s`;
my$tmSeqObj=new Bio::SeqIO(-file=>"$tmFileSeq",-format=>"fasta");
while(my$tmSeqO=$tmSeqObj->next_seq){my$tmTimeOld=`date +%s`;
my$ttmSeq=$tmSeqO->seq;my$tmSeq="\U$ttmSeq\E";my$tmID=$tmSeqO->id;
my$tmFile="";my$tmFSeq="";my$tmCmdLine="";my$tmFileDat="";my$tmFileTIR="";
$tmFSeq=&efFormatSeq($tmSeq,$egLineLength,1);
$tmFSeq=~s/ +//g;
$tmFile="$egDirTemp/$tmPID.temp_seq.raw";
open(tmOut,">$tmFile")or die"Error when saving the temporary sequence!\n";
print tmOut"$tmFSeq\n";
close(tmOut);
$tmFileDat="$egDirTemp/$tmPID.temp_seq.dat";$tmFileTIR="$egDirTemp/$tmPID.temp_tir.dat";
$tmCmdLine="\"$egDirScript\/$egScriptMUST\" \"$tmFile\" \"$tmID\" \"$tmFileDat.before_filter\" \"$tmFileTIR\" $gMinTIR $gMaxTIR $gMinDR $gMaxDR $gMinMITE $gMaxMITE $gFixedFlanking $gMutationRate";
if(system("$tmCmdLine >/dev/null")){die"-Error occurred when running the script \"$egDirScript\/$egScriptMUST\"!\n";}$tmCmdLine="\"$egDirScript/$egScriptFilter\" \"$tmFileDat.before_filter\" \"$tmFileDat\"";
if(system("$tmCmdLine 2>/dev/null 1>/dev/null")){die"-Error occurred when running the script \"$egDirScript\/$egScriptFilter\"!\n";}}my$tmgTimeNew=`date +%s`;}sub efValidString{my($tStr)=@_;
if((defined$tStr)and($tStr ne"")){return$tStr;}else{return"";}}sub efSyntax{my$egSyntax=decode_base64("TUlURSBVbmNvdmVyaW5nIFN5c1RlbSAgICB2ZXJzaW9uIDIuNC4wMDEKLS0tLS0tLS0tLS0tLS0t
LS0tLS0tKGMpIEZlbmdmZW5nIFpob3UgKEZlbmdmZW5nWmhvdUBnbWFpbC5jb20pIDIwMTMtMTIt
MjgKLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0t
LS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCi4vTVVTVF9QaXBlLnBsIDxpbnB1dC5mYXN0YT4g
PG91dHB1dC5NSVRFLmRhdD4gPERpclRlbXA+IFtvcHRpb25zXQpOb3RlOgpJZiBhbnkgb3B0aW9u
cyBhcmUgcHJvdmlkZWQsIHRoZSB1c2VyIGhhdmUgdG8gaW5wdXQgdGhlIHZhbHVlcyBmb3IgYWxs
IHRoZSBvcHRpb25zLgpPcHRpb24gbGlzdDoKLTxUaHJlYWRfTnVtPiAgICAgICAgICAgICBbMTBd
Ci08TUlOX1RJUl9sZW5ndGg+ICAgICAgICAgWzhdCi08TUFYX1RJUl9sZW5ndGg+ICAgICAgICAg
WzUwXQotPE1JTl9EUl9sZW5ndGg+ICAgICAgICAgIFsyXQotPE1BWF9EUl9sZW5ndGg+ICAgICAg
ICAgIFszMF0KLTxNSU5fTUlURV9sZW5ndGg+ICAgICAgICBbMTAwXQotPE1BWF9NSVRFX2xlbmd0
aD4gICAgICAgIFs2MDBdCi08RklYRURfRkxBTktJTkdfbGVuZ3RoPiAgWzUwXQotPE11dGF0aW9u
X1JhdGU+ICAgICAgICAgIFswLjgwXQotLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0t
LS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQo="
);
print$egSyntax;}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 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;}}sub efIsSameStrand{my($iSeq1,$iSeq2)=@_;
if((length($iSeq1)>0)and(length($iSeq2)>0)){my$iF1="$egDirTemp/seq1.fasta";
open(iO,">$iF1")or die"Error when saving sequence 1\n";
print iO">Query\n$iSeq1\n";
close(iO);
my$iF2="$egDirTemp/seq2.fasta";
open(iO,">$iF2")or die"Error when saving sequence 2\n";
print iO">DB\n$iSeq2\n";
close(iO);
my$iFO="$egDirTemp/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 efPrintData{my(%pIdx)=@_;
my@pKeys=keys%pIdx;
for(my$pi=0;$pi<@pKeys;$pi++){print"--[$pKeys[$pi]] => [".$pIdx{"$pKeys[$pi]"}."]---\n";}print"----------------------------------------------------------------------\n";}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;}
