#!/usr/bin/perl -w

use strict;

############################################################################
##    With this script we can classify the output of the BLAST analysis   ## 
##        (predicted U12 introns) into five different categories          ##
############################################################################

# define the variables

my $VERBOSE = 0;

my @blast = "";                      # array to store every line using tabs (this will read second line)
                                     # we store every position of the array into a variable
     my $query = 0;                  # (esay to work with, easy to correct)

     my $subject = 1;

     my $identity = 2;

     my $aliglenght = 3;

     my $mismt = 4;

     my $gaps = 5;

     my $qstart = 6;

     my $qend = 7;

     my $sstart = 8;

     my $send = 9;

     my $evalue = 10;

     my $score = 11;

my @previousblast = "";                    # array to read first line of the input

my @out = "";                              # array to store the output

my $SUBJECT = "";                          # used to compare if we are working with the same subject

my $entrada = $ARGV[0];

my $llegint = 0;                           # for go to another line ($reading, in english)

my $evalue_thresh = 1.0;                   # e-value threshold

my $n = 0;                                 # when it's ==1, it means that we have classified one query

$~ = 'OUTPUT_INTRON';                      # to store the output into a file

############################################################################################

my @outname = split (/\./,$entrada);            # output file will be: name of the category plus the 
                                                # name of the input file without the extension.
                                                # we split the name of the input into an array
unlink <*.outname[0]>;

open (FILE,$entrada) || die "No trobo l'arxiu $entrada\n";    # open or die!!


while (<FILE>) {

	if (( m/^\w+/ )&&($n == 1)){                  # skip those hits whose 
	     next;                                    # query has been already classified ($n == 1)
 	}

	elsif ( m/^#/){                               # skip comment lines and restart $n and $llegint
	
     		$n = 0;
     		$llegint = 0;
     		next;
  	}


	chomp $_;

	@blast = split (/\t/,$_);           # split the line of the input into an arrat usin tabs


# some parameters should be acomplished to consider one hit:

	if (($llegint > 0) && ($previousblast[$evalue] < $evalue_thresh ) && ($previousblast[$gaps] < 4)) {
    
		$SUBJECT = $previousblast[$subject];   # compare if we are in the same subject

   		if (($blast[$subject] eq $SUBJECT) 
	   		&& 
       		($previousblast[$identity] >= 97.0)
	   		&& 
       		($blast[$identity] >= 97.0)
	   		&&
      		((($previousblast[$qend] >= 41) && ($previousblast[$qend] <= 60) && ($blast[$qstart] >= 41) && ($blast[$qstart] <= 60))
           		||
      		(($previousblast[$qstart] >= 41) && ($previousblast[$qstart] <= 60) && ($blast[$qend] >= 41) && ($blast[$qend] <= 60)))) 
	   
	     
		{
			open (OUTPUT_INTRON,">>2Exons.$outname[0]");          # put this line into 
			select OUTPUT_INTRON;                                 # corresponding category
	   		$out[0] = $blast[$query];                             # Store the ouput into
	   		$out[1] = $blast[$subject];                           # an array
	   		$out[2] = "Two exons";
	   		$out[3] = (($blast[$identity]) + ($previousblast[$identity]))/200;
	 
		}

    		elsif (($previousblast[$identity] >= 97.0) 
	       		&&
	  	($previousblast[$qstart] > 0) && ($previousblast[$qstart] <= 36) 
	       		&& 
	  	($previousblast[$qend] >= 65) && ($previousblast[$qend] <= 100)) 
		{

                        open (OUTPUT_INTRON,">>aligjunc.$outname[0]");
                        select OUTPUT_INTRON;
 	    		if ($VERBOSE){print STDERR "llegint $blast[1]\n";}
	    		if ($VERBOSE){print STDERR "llegint $previousblast[1]\n";}
	    		$out[0] = $previousblast[$query];
	    		$out[1] = $previousblast[$subject];
	    		$out[2] = "Alignment across junction";
	    		$out[3] = ($previousblast[$identity])/100;
	 
		}


     		elsif (($previousblast[$identity] >= 97.0) 
		 	&&
	    	($previousblast[$qstart] > 0) && ($previousblast[$qstart] <= 36) 
		 	&& 
	    	($previousblast[$qend] >= 41) && ($previousblast[$qend] <= 60)) 
		
		{

                        open (OUTPUT_INTRON,">>1exonleft.$outname[0]");
                        select OUTPUT_INTRON;
	    		if ($VERBOSE){print STDERR "llegint $blast[1]\n";}
	    		if ($VERBOSE){print STDERR "llegint $previousblast[1]\n";}
	    		$out[0] = $previousblast[$query];
	    		$out[1] = $previousblast[$subject];
	    		$out[2] = "One exon flank-left";
	    		$out[3] = ($previousblast[$identity])/100;
	 
		}


    		elsif (($previousblast[$identity] >= 97.0) 
			&&
	  	($previousblast[$qstart] >= 41) && ($previousblast[$qstart] <= 60) 
			&& 
	  	($previousblast[$qend] >= 65) && ($previousblast[$qend] <= 100))

		{

                        open (OUTPUT_INTRON,">>1exonrigth.$outname[0]");
                        select OUTPUT_INTRON;
	   		if ($VERBOSE){print STDERR "llegint $blast[1]\n";}
	   		if ($VERBOSE){print STDERR "llegint $previousblast[1]\n";}
	    		$out[0] = $previousblast[$query];
	    		$out[1] = $previousblast[$subject];
	    		$out[2] = "One exon flank-right";
	    		$out[3] = ($previousblast[$identity])/100;
	    
	   
		} else {
                        open (OUTPUT_INTRON,">>nohit.$outname[0]");
                        select OUTPUT_INTRON;
			$out[0] = $blast[$query];
			$out[1] = " ";
			$out[2] = "No hits found!";
			$out[3] = " ";
		
		}
 		
                    # for detailed info about all the categories see 'scripts' section of the webpage 
	
			write;
			undef (@out);                   # undefine @out when we change the query
		       
    		$n = 1;                                 # to define that one query has been organized
	       
 	}
 
	@previousblast = @blast;                 # another array to read two lines at the same time
 
        $llegint = $llegint + 1;                 # go to next line
 }

close(FILE);

# we will use 'format' command to organize the output in columns

format OUTPUT_INTRON = 
@<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<<
$out[0],                                               $out[1],						$out[2],			$out[3]
.
