#!/usr/bin/perl

# file: RBS.pl
# Author: Xingyi Guo
#use strict -w;
#use FileHandle;
use Getopt::Std;
%options=();
getopt("iIsSobgdexy",\%options);

use constant USAGE =><<END;
Usage: $0 -i <chrA_cds_fasta_seq> -I <chrB_cds_fasta_seq> -s <chrA_site> 
-S <chrB_site> -o <OUTPUT> -b <blast_path>

++++++++++++++++++++++++ RBS Version:1.0 ++++++++++++++++++++++++++++++++

 Options
 -i input chromosome A FASTA cds sequence          [*]     
 -I input chromosome B FASTA cds sequence          [*]
 -s input chromosome A FASTA cds site information  [*]     
 -S input chromosome B FASTA cds site information  [*]
 -o get output file                                [*]
 -b set absolute blast PATH                        [*]
 -g use gnuplot program to plot data,get *.ps file [default: T]
 -d get output data, you can use excel to plot     [default: F]
 -e Blast E-value                                  [default: 1e-14]   
 -x set xlable name for gnuplot                    [default: <chrA_cds_fasta_seq>]
 -y set ylable name for gnuplot                    [default: <chrB_cds_fasta_seq>]

***Attention: parameter [*] means that you must give the specific value 

Note: For a detailed description of how to use RBS_v1.0.pl,please read the README file
===============================================================================
Please address questions and bug reports to: bioinplant\@zju.edu.cn
homepage: http://ibi.zju.edu.cn/bioinplant/
last Revised JUN 6,2005

END


die USAGE if $options{'i'} eq "";
die USAGE if $options{'I'} eq "";
die USAGE if $options{'s'} eq "";
die USAGE if $options{'S'} eq "";
die USAGE if $options{'S'} eq "";
die USAGE if $options{'o'} eq "";

$ARGV[0]=$options{'i'};
$ARGV[1]=$options{'I'};
$ARGV[2]=$options{'s'};
$ARGV[3]=$options{'S'};
$ARGV[4]=$options{'o'};

$options{'d'}= "F" if($options{'d'} eq "");
$options{'g'}= "T" if($options{'g'} eq "");
$options{'e'}= "1e-14" if($options{'e'} eq "");
$options{'x'}= "$ARGV[0]" if($options{'x'} eq "");
$options{'y'}= "$ARGV[1]" if($options{'y'} eq "");

my $blastwarn="Warning! please check your blast path\n";
my $Error="Error: you should alternative -d T or -g T\n"; 

die $blastwarn if $options{'b'} eq "";
die $Error if ($options{'d'} eq "F" && $options{'g'} eq "F");



#######if blast path is not corect,please modify the blast path##################
my $blast_path="$options{'b'}";#/usr/local/bin/blast or /home/biotools/blast

##########query vs database(blastn e-14),if necessary,you can change the standard of E-value

system("$blast_path/formatdb -p F -o T -i ./$ARGV[1]");
system("$blast_path/blastall -p blastn -i ./$ARGV[0] -d ./$ARGV[1] -o ./RBS_blastn.out -v 1 -b 1 -F F -e $options{'e'}");
#query hits out
my %hash_hits1=&BLASTOUT("./RBS_blastn.out");
system("rm -rf ./RBS_blastn.out");
system("rm -rf ./$ARGV[1].*");
system("$blast_path/formatdb -p F -o T -i ./$ARGV[0]");
system("$blast_path/blastall -p blastn -i ./$ARGV[1] -d ./$ARGV[0] -o ./RBS_blastn.out -v 1 -b 1 -F F -e $options{'e'}");

#database hits out
my %hash_hits2=&BLASTOUT("./RBS_blastn.out");
#clean blastout
system("rm -rf ./RBS_blastn.out");

system("rm -rf ./$ARGV[0].*");
system("rm -rf ./$ARGV[1].*");

##############get cds site infomation#############################
my %siteA=&SITE("./$ARGV[2]");
my %siteB=&SITE("./$ARGV[3]");

################get data for Excel ###############
if($options{'e'} eq "T")
{
        open EXCEL,">./$ARGV[4].exl";
        open EXCEL_REV,">./$ARGV[4]_Rev.exl";
}


###########use gnuplot program to plot data ####################
if($options{'g'} eq "T")
{
        open OUT,">./$ARGV[4]";
        open OUT1,">./$ARGV[4]_reverse";
        print OUT "#-- FORWARD MATCHES\n0\t0\n\n";
        print OUT1 "#-- REVERSE MATCHES\n0\t0\n\n";
}
foreach my $key(sort keys %hash_hits1)
{
        
        if(exists $hash_hits2{"$hash_hits1{$key}"})
        {
                my @tmp_site1=split/\s+/,$siteA{$key};
                my @tmp_site2=split/\s+/,$siteB{$hash_hits1{$key}};
                my $tmp_max1=$tmp_site1[1]-$tmp_site1[0];
                my $tmp_max2=$tmp_site2[1]-$tmp_site2[0];
                my @newsite1=@tmp_site1;
                my @newsite2=@tmp_site2;
                if($tmp_max1 < 0)
                {
                        @newsite1="";
                        @newsite1=&MAX(@tmp_site1);
                        my $tmp_flag1=1;
                 }
                 if($tmp_max2 < 0)
                {
                        @newsite2="";
                        @newsite2=&MAX(@tmp_site2);
                        $tmp_flag2=1;
                 }
                 $A_max=$newsite1[1] if $A_max < $newsite1[1];
                 $B_max=$newsite2[1] if $B_max < $newsite2[1];
                 if(($tmp_max1 < 0 && $tmp_max2 > 0) || ($tmp_max1 > 0 && $tmp_max2 < 0))
                 {
                        print EXCEL_REV "$newsite1[0]\t$newsite2[0]\n" if ($options{'e'} eq "T");
                        print OUT1 "$newsite1[1]\t$newsite2[0]\n$newsite1[0]\t$newsite2[1]\n\n" if($options{'g'} eq "T");
                }
                 if(($tmp_max1 > 0 && $tmp_max2 > 0) || ($tmp_max1 < 0 && $tmp_max2 < 0))
                 {
                        print EXCEL "$newsite1[0]\t$newsite2[0]\n" if ($options{'e'} eq "T");
                        print OUT "$newsite1[0]\t$newsite2[0]\n$newsite1[1]\t$newsite2[1]\n\n" if($options{'g'} eq "T");
             
                 }
                
        }
#        if($hash_hits1{$key}=~/^(.+)reverse/i)
#        {
#                my $tmp_match=$1;
#                if(exists $hash_hits2{$tmp_match})
#                {
#                        my @tmp_site1=split/\s+/,$siteA{$key};
#                        my @tmp_site2=split/\s+/,$siteB{$tmp_match};
#                        my $tmp_max1=$tmp_site1[1]-$tmp_site1[0];
#                        my $tmp_max2=$tmp_site2[1]-$tmp_site2[0];
#                        my @newsite1=@tmp_site1;
#                        my @newsite2=@tmp_site2;
#                        if($tmp_max1 < 0)
#                        {
#                                 @newsite1=&MAX(@tmp_site1);
#                         }
#                         if($tmp_max2 < 0)
#                        {
#                                @newsite2="";
#                                @newsite2=&MAX(@tmp_site2);
#                         }
#                        $A_max=$newsite1[1] if $A_max < $newsite1[1];
#                        $B_max=$newsite2[1] if $B_max < $newsite2[1];
#                         print EXCEL_REV "$newsite1[0]\t$newsite2[0]\n" if ($options{'e'} eq "T");
#                        print OUT1 "$newsite1[1]\t$newsite2[0]\n$newsite1[0]\t$newsite2[1]\n\n" if($options{'g'} eq "T");
#                  
#                }
#        }
}

if($options{'g'} eq "T")
{
        print OUT "\n";
        system("rm -rf ./$ARGV[4].plot");
        system("cat ./$ARGV[4] ./$ARGV[4]_reverse >> ./$ARGV[4].plot");
        system("rm -rf ./$ARGV[4]");
        system("rm -rf ./$ARGV[4]_reverse");
        open GP,">$ARGV[4].gp";
        print GP "set terminal postscript color solid\n";
        print GP "set output \"$ARGV[4].ps\"\n";
        print GP "set nokey\nset nogrid\nset xlabel \"$options{'x'}\"\n";
        print GP "set ylabel \"$options{'y'}\"\n";
        print GP "set xrange [0:$A_max]\n";
        print GP "set yrange [0:$B_max]\n";
        print GP "plot \"$ARGV[4].plot\" index 0 title \"Fwd\" w l lt 1 lw 2, \"$ARGV[4].plot\" index 1 title \"Rev\" w l lt 2 lw 2\n";
        system("gnuplot $ARGV[4].gp");
}
###########get site information#############################
sub SITE
{
      
        my %hash_site;
        my $site_file=shift;
        open SITE_FILE,"$site_file";
        while(<SITE_FILE>)
        {
                
                ########if your site file do not meet the format,please change the follow matching 
#                if(/(.+?)\|.+\s+(\d+)/)
#                {
#                        $hash_site{$1}=$2;
#                }
                  split;
                  $hash_site{$_[0]}="$_[1] $_[2]";
                 
        }
        return %hash_site;
}   
#########Deal with the blast out############################

sub BLASTOUT
{
        
        my %hash_hits;
        my $blastout=shift;
        open BOUT_FILE,"$blastout";
        while(<BOUT_FILE>)
        {
               
                if(/query=\s+(.+?)\#/i)
                {
                        our $query=$1;
                }
                if(/^>(.+?)\#/i)
                {
                       our $database=$1;
                       $hash_hits{$query}="$database";
                   #    our $flag=1;
                }
#                if($flag == 1)
#                {
#                        if(/query:.+?(\d+).+?(\d+)/i)
#                        {
#                               my $lastsite=$2;
#                               my $firstsite=$1;
#                               our $tmp=$lastsite - $firstsite;
#                        }
#                         if(/sbjct:.+?(\d+).+?(\d+)/i)
#                        {
#                               my $lastsite1=$2;
#                               my $firstsite1=$1;
#                               my $tmp1=$lastsite1-$firstsite1;
#                               if(($tmp > 0 && $tmp1 < 0) || ($tmp1 > 0 && $tmp < 0))
#                               {
#                                        $hash_hits{$query}="$database"."reverse";
#                                       
#                               }
#                               else
#                               {
#                                        $hash_hits{$query}="$database";
#                               }
#                               $flag=0;    
#                        }
#                }
        }
        return %hash_hits;
}
################change the first site with the last site#####################
sub MAX
{
        my @data=@_;
        my $tmp;
        $tmp=$data[1];
        $data[1]=$data[0];
        $data[0]=$tmp;
        return @data;
}
        