#!/usr/bin/perl -w

use strict;

system ("clear");
print STDERR "----PATTERNmatching RESULTS----","\n";
print STDERR "By Belloc E, Duch A, Lois S, Moldon A.","\t","UPF -Universitat Pompeu Fabra (Barcelona)-","\n","\n";
print STDERR "WAIT A MOMENT PLEASE...","\n";
########################################################################################
################################ VARIABLE DECLARATION ##################################
########################################################################################
 
my $seqfile;
my $patfile;
my $exonfile;
my $n;
my $i;
my $j;   
my $f;
my $m;
my $first_line;
my $seq;
my $count_in;
my $count_out;
my $long;
my $percent_in;
my $percent_out;
my $add;

my @pat;
my @seq_matrix;
my @pat_matrix;
my @exon;
my @exon_matrix;
my @v;
my @results;
my @match_matrix;


########################################################################################
################################# OPENING FILES ########################################
########################################################################################

# If the number of arguments introduced to run the program is smaller than 3 more files are needed to run the program


if (scalar(@ARGV) < 3) {
  print STDERR "USAGE: patternmatching.pl refseq.aa.fa prosite.tbl refseq.aa.exons\n";
  exit(1);
};

$seqfile = $ARGV[0];
$patfile = $ARGV[1];
$exonfile = $ARGV[2];

# The specified file in the first argument is opened from the shell and passed into the program 

open(SEQ,"< $seqfile") || do {
  print STDERR "patternmatching.pl: Not found sequence file $seqfile\n";
  exit(1);
};

# The specified file in the second argument is opened from the shell and passed into the program

open(PAT,"< $patfile") || do {
  print STDERR "patternmatching.pl: Not found pattern file $patfile\n";
  exit(1);
};

# The specified file in the third argument is opened from the shell and passed into the program

open(EXON,"< $exonfile") || do {
  print STDERR "patternmatching.pl: Not found exon file $exonfile\n";
  exit(1);
};

# Open a file to write the output (matches found)
open(OUTFILE,"> outfile.txt");

########################################################################################
############################ PREPARING DATA ############################################
########################################################################################

# Create a matrix of sequences from a file in fasta format (in $seqfile) in order to have:
#   -First column: sequence identificator
#   -Second column: string (protein sequence)

$/ = '>';                                                  
$first_line = <SEQ>;                                                 
$i = 0;
while (<SEQ>) {
     chomp;
     @v = split(/\n/, $_);
     $seq_matrix[$i][0] = $v[0];
     $j = 1;
     $seq_matrix[$i][1] = '';
     while ($j < scalar(@v)) {
         $seq_matrix[$i][1] = $seq_matrix[$i][1].$v[$j];
   
	 $j = $j + 1;
     }
     $i = $i + 1;
}
close (SEQ);

# Create a matrix of patterns from a file downloaded from prosite database (in $patfile): 
#   -First column: pattern identificator
#   -Second column: regular expression defining an specific pattern

$/ = "\n";
$i = 0;
while (<PAT>) {
    $_=~s/^  //g;
    chomp ($_);
    @pat = split(/ /,$_);
    $pat_matrix[$i][0] = $pat[0];
    $pat_matrix[$i][1] = $pat[1];
    $i = $i + 1;
}
close (PAT);

# Create a reference exon matrix (in $exon file):
#   -First column: sequence identificator
#   -Second column: number of exon
#   -Third column: initial position of the exon
#   -Fourth column: last position of the exon 

$i = 0;
while (<EXON>) {
    chomp;
    @exon = split (/ /, $_);
    $exon_matrix[$i][0] = $exon[0];
    $exon_matrix[$i][1] = $exon[1];
    $exon_matrix[$i][2] = $exon[2];
    $exon_matrix[$i][3] = $exon[3];
    $exon_matrix[$i][1] =~s/NM_[0-9]*\.//g;
    $i = $i + 1;
}
close (EXON);

########################################################################################
################################# MAIN ALGORITHM #######################################
########################################################################################

# Searches patterns in all the sequences and creates a matching matrix with:
#   -First column: pattern identificator
#   -Second column: sequence identificator
#   -Third column: first position of the found match in the sequence    
#   -Fourth column: last position of the found match in the sequence
# Create file containing a table with: pattern identificator, sequence identificator, total number of matches per each pattern, initial position of the match, length of the match, sequence of the found match
 

$j=0;

$n = 0;

print OUTFILE "----PATTERNmatching RESULTS----","\n";
print OUTFILE "By Belloc E, Duch A, Lois S, Moldon A.","\t","UPF -Universitat Pompeu Fabra (Barcelona)-","\n","\n"; 
while ($j < scalar(@pat_matrix)) {
    $i = 0;
    while ($i < scalar(@seq_matrix)) {
   	@results = &match($seq_matrix[$i][1], $pat_matrix[$j][1]);
        if (scalar(@results) > 0) {
	       print OUTFILE "Id. PATTERN","\t",$pat_matrix[$j][0],"\n", "Id. SEQUENCE","\t",$seq_matrix[$i][0],"\n";
             print OUTFILE "Total MATCHES","\t",scalar(@results)/3,"\n";
	    $f = 0;
	    while ($f < scalar(@results)) {
                 print OUTFILE  "MATCH",$f/3+1,"\t","\t", $results[$f]."(".$results[$f+1].")","\t","...".$results[$f+2]."...","\n";
            $match_matrix[$n][0]=$pat_matrix[$j][0];
		$match_matrix[$n][1]=$seq_matrix[$i][0];
		$match_matrix[$n][2]=$results[$f];
		$match_matrix[$n][3]=$results[$f]+$results[$f+1]-1;

		$n = $n + 1;
		$f = $f + 3;		
            }
	    
        }

        $i = $i + 1;
    }

    $j = $j + 1;
}

# Compares the matching matrix with the exon matrix in order to determine whether the match is produced intraexonically ($count_in) or interexonically ($count_out)

print STDERR "MATCHES FOUND: ",scalar(@match_matrix),"\n";

$count_in = 0;
$count_out = 0;
$long = 0;
$percent_in = 0;
$percent_out = 0;
$add = 0;
$i = 0;
$n = 0;

# In first place we read along the matching matrix 

my $OUTFORMAT = "%-8s %5d    %6d (%5.1f)  %6d (%5.1f)  %10d\n";
            #          1         2         3         4         5         6         7         8         9        10
            # 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
            # ssssssss ddddd    dddddd (fff.f)  dddddd (fff.f)  dddddddddd
print STDOUT "Pattern   Length  IntraExonic(%)  InterExonic(%)  TotalMatches\n" ;
print STDOUT "**************************************************************\n";
while ($n < scalar(@match_matrix)) {

    $j = 0;

    # Second, we read along the matching matrix until the sequence Id ($match_matrix[$n][1]) matches with the sequence Id of the exon matrix ($exon_matrix[$j][0]) 

    while ($j < scalar(@exon_matrix) && $match_matrix[$n][1] ne $exon_matrix[$j][0]) {
	$j = $j + 1;
    }

    if ($j < scalar(@exon_matrix)) {

	    # When the Ids match, the algorithm will search in which exon positions the pattern matches are located.
            # $count_in: counts how many times the pattern match fits exactly in an exon 
            # $count_out: counts how many times the pattern match does not fit exactly in an exon    

	    while ($j < scalar(@exon_matrix) && $match_matrix[$n][1] eq $exon_matrix[$j][0]) {

		if (($match_matrix[$n][2] > $exon_matrix[$j][3])){
		}
		else{
		    if (($match_matrix[$n][2] >= $exon_matrix[$j][2])  && ($match_matrix[$n][3] <= $exon_matrix[$j][3])){
			$count_in = $count_in + 1;
		    }
		    else{
			if (($match_matrix[$n][2] >= $exon_matrix[$j][2])  && ($match_matrix[$n][3] > $exon_matrix[$j][3])) {
			    $count_out = $count_out + 1;
			}
		    }
		}
	        $j = $j + 1;
	    }
    }
    else {
	print STDERR " $match_matrix[$n][1] PROTEIN NOT FOUND IN EXON MATRIX!!!\n";
        $m = $m + 1;
    }
    

    if ($n+1 >= scalar(@match_matrix)) {
	$long=$match_matrix[$n][3]-$match_matrix[$n][2]+1;
	$long= $match_matrix[$n][3]-$match_matrix[$n][2]+1;
	$percent_in= ($count_in/($count_in+$count_out)*100);
	$percent_out=($count_out/($count_in+$count_out)*100);
	$add= $count_in + $count_out;
	printf STDOUT $OUTFORMAT, $match_matrix[$n][0], $long, $count_in, $percent_in, $count_out, $percent_out, $add;
    }
    else {
	if ($match_matrix[$n][0] ne $match_matrix[$n+1][0]) {
	    $long= $match_matrix[$n][3]-$match_matrix[$n][2]+1;
	    $percent_in= ($count_in/($count_in+$count_out)*100);
	    $percent_out= ($count_out/($count_in+$count_out)*100);
	    $add= $count_in + $count_out;
	    printf STDOUT $OUTFORMAT, $match_matrix[$n][0], $long, $count_in, $percent_in, $count_out, $percent_out, $add;
	    $count_in = 0;
	    $count_out = 0;
        }
    }
    $n = $n + 1;   
}


exit(0);

########################################################################################
################################## FUNCTION ############################################
########################################################################################

# function which finds matches between patterns and sequences
# returns: first position of the match ($a), length of the match ($b), and sequence match

sub match() {
    my ($seq,$pat) = @_;
    my @out = ();
    while ($seq =~ m/$pat/g) {
        my ($a,$b,$c);
        $a = length($`) + 1;
        $b = length($&);
        $c = $&;
        push @out, ($a, $b, $c);
    };
    return @out;
}
