#!/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;
@ARGV==3 or die"Error in syntax!\n"."   ./filterSSR.pl <cutoff> <input.fasta> <output.valid_id.txt>\n";
my$egFileI=$ARGV[1];
my$egFileO=$ARGV[2];
my$egCutOff=$ARGV[0];
my$egScriptTRF="./MUST/trf400.linux.exe";my$egParamTRF="2 7 7 75 20 50 500 -h -d";my$egParamTRFTag="2.7.7.75.20.50.500";
my$egScriptSPU="./MUST/sputnik.program";my$egParamSPU="";my$egParamSPUTag="";
my%egID2Pos=();my@egID=();my@egSeq=();my@egStatus=();my$egNumber=0;
print"Loading the sequences ... ";
my$egSeqObj=new Bio::SeqIO(-file=>"$egFileI",-format=>"fasta");
while(my$egSeqO=$egSeqObj->next_seq){my$tID=$egSeqO->id;my$tDesc=$egSeqO->desc;
my$tSeq=$egSeqO->seq;
if((defined$tID)and(defined$tSeq)and($tSeq ne"")){$egID2Pos{"$tID"}=$egNumber;$egID[$egNumber]=$tID;$egSeq[$egNumber]=$tSeq;$egStatus[$egNumber]="";$egNumber++;}}print" [Number:$egNumber] [done]\n";
my$egCmdLine="";
my$egDirData="./tempData";
if(-d$egDirData){}elsif(system("mkdir -p \"$egDirData\"")){die"-Error when initizting!\n";}print"Scanning for the SSR using TRF ... ";
chdir($egDirData);
$egCmdLine="\"../$egScriptTRF\" \"../$egFileI\" $egParamTRF";
if(system("$egCmdLine 1>/dev/null 2>/dev/null")){}else{die"-Erro when TRF!\n";}chdir("..");
my$tFile="$egFileI.$egParamTRFTag.dat";
$tempi=0;
$tid="";
if(-e$tFile){open(efI,"$tFile")or die"Error!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\s+/,$tline);
if($tline=~/^Sequence: (.*)$/){$tid=$egID2Pos{"$1"};
if(defined$tid){;}else{$tid="";}}elsif((@tlist>=15)and(defined$tid)and($tid ne"")){$egStatus[$tid]=$egStatus[$tid].";$tlist[0]-$tlist[1]|$tlist[13]|TRF";
$tempi++;}}close(efI);}print" [done] [TRF:$tempi]\n";
if(1==0){print"Scanning for the SSR using sputnik ... ";
chdir($egDirData);
$egCmdLine="\"../$egScriptSPU\" \"../$egFileI\" $egParamSPU";
if(system($egCmdLine)){die"-Erro when SPU!\n";}chdir("..");
my$tFile="$egFileI.$egParamSPUTag.dat";
$tempi=0;
$tid="";
if(-e$tFile){open(efI,"$tFile")or die"Error!\n";
while(<efI>){$tline=$_;$tline=~s/[\r\n]//g;@tlist=split(/\s+/,$tline);
if($tline=~/^>(.*)$/){$tid=$egID2Pos{"$1"};
if(defined$tid){;}else{$tid="";}}elsif(($tline=~/^[a-zA-Z\d]+ (\d+) : (\d+) --/)and(defined$tid)and($tid ne"")){my($pS,$pE)=($1,$2);
$egStatus[$tid]=$egStatus[$tid].";$pS-$pE|sputnik";
$tempi++;}}close(efI);}print" [done] [SPU:$tempi]\n";}print"Saving the final ID list ... ";
$tempk=0;
my$tLength=0;my$tSSRLength=0;
my($tS,$tE,$tUnit,$tProg)=(0,0,"","");
open(efO,">$egFileO")or die"Error when saving!\n";
for($tempi=0;$tempi<$egNumber;$tempi++){@tlist=split(/;/,$egStatus[$tempi]);
$tLength=length($egSeq[$tempi]);
$tSSRLength=0;
for($tempj=0;$tempj<@tlist;$tempj++){($tS,$tE,$tUnit,$tProg)=split(/[-\|]/,$tlist[$tempj]);
if((defined$tS)and(defined$tE)and(defined$tUnit)and($tUnit ne"")){$tSSRLength=$tSSRLength+abs($tE-$tS)+1;}}if($tSSRLength>=$egCutOff*$tLength){$tempk++;}else{print efO"$egID[$tempi]\n";}}close(efO);
print" [done] [RemovedNumber:$tempk]\n";
