#!/usr/bin/perl -w 

###########################################################################################
###########################################################################################
###########################################################################################
#                                                                                         
#     SCANNER: TRANSCRIPTIONAL FACTORS MOTIFS MULTIALIGNMENT.                             
#                                                                                         
#     March 2002 - Guiomar SOLANAS                                                        
#                  Alex VENDRELL                                                          
#                                                                                         
# This program is free software; you can redistribute it and/or modify
# it.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
#
#
###########################################################################################
#
# Creation Date: 10.02.2002
# 
# Last Update: 2002/03/01  
#
# Autors : Guiomar SOLANAS
#          e-mail: guiomar.solanas01@campus.upf.edu
#          Alexandre VENDRELL
#          e-mail: alexandre.vendrell01@campus.upf.edu
# 
#
#
###########################################################################################
###########################################################################################
########################################################################################### 

########################################################################################### 
############# DECLARATION OF THE VARIABLES WE ARE GOING TO USE IN THE PROGRAM ############# 
########################################################################################### 

use strict; 
use Getopt::Std; 

my $fitxer; 
my $matriu; 
my $locussigne; 
my $locus; 
my @totesseq;
my $numero_seq;
my @nomsdeseq;
my @seq;
my $novaseq; 
my $treupunt; 
my @definitivaseq; 
my @vectordeseq;
my $cmesg;
my $compcmesg; 
my $liniamatriu; 
my @nomsdematriu;
my $nommatriu;
my $poslletra;
my $num_matriu_actual;
my $yo;
my @tu;
my $seqactual;
my $numero_matriu;
my @vectormatriu;
my $posseq; 
my $fila; 
my @filera;
my $columna; 
my @puntuacio; 
my $consens;
my @score; 
my $puntuacio_final; 
my $posmatriu; 
my $sum; 
my $tamany_matriu; 
my $lletra;
my $comptapos; 
my @posicions;
my $impressio;
my $llindar;
my $tot; 
my %options;
my $dibuseq;
my @pintaseq;
my $printa;
my $dibumatrix;
my $nouprinta;
my @filamat;
my @color;

##########CONSTANT DECLARATION########################################################
 
my $ZERO=0.0000000000000000000001;
my $A=0;
my $C=1;
my $G=2;
my $T=3;

####################################################################################### 
######################### DECLARATION OF THE GETOPTS AND ROUTINES ##################### 
#######################################################################################

# Getopts hash.

(getopts('-h:vmst:',\%options)) or die ("Problems reading the options, try  again\n"); 

# Rutine for print_mess.

sub print_mess
{
    my @mess=@_;
    print STDERR ("\n\t- SCANNING: @mess \n") if (exists($options{v}));
} 
sub print_error
{
    my @mess=@_;
    print STDERR ("\n\t- SCANNER: @mess \n") and exit(1);
} 
####################################################################################### 
######################## OPENING OF OUTFILES ##########################################
#######################################################################################

# Opening for writing (>).

open (OUTFILE,">outfile.html");

print OUTFILE "<html>   <head>    <title>    OUTFILE SCANNER    </title>
   </head>    <BODY    TEXT color ='483D8B'    BGCOLOR='EEE8AA' link='black' vlink='black' alink='black'><br><br><font color='8B0000'> <font size='10'><center><blockquote><b> SCANNER: TRANSCRIPTIONAL FACTORS MOTIFS MULTIALIGNMENT RESULTS.</font></b></font></center><br><p></blockquote><font color='8B0000'><font size='6'> RESULTS. <hr> </font></font><br><p>\n";


#######################################################################################
################################# OPTION HELP ######################################### 
#######################################################################################

if (exists($options{h})) { 
    
    if ($options{h} == 1) { 
	print "\nOPCIONES:\n\t -h: si el usuario quiere recibir ayuda acerca del funcionamineto del programa.\n\t -v: si el usuario quiere recibir  informacion mientras el programa esta corriendo (en ingles).\n\t -m: si el usuario necesita mas informacion sobre la matriz de pesos (PWM) y el consenso de esta.\n\t -s: si el usuario necesita mas informacion sobre la secuencia (nombre, contenido en G+C, longitud)\n\t -txx: si el usuario quiere fijar un valor umbral en los resultados (si no, por defecto sera 80).\n\n";
    } 
    
    if ($options{h} == 2) { 
	print "\nOPCIONS:\n\t -h: si l'usuari vol rebre ajuda sobre el funcionament del programa.\n\t -v: si l'usuari vol rebre informacio mentre el programa corre (en anglès)\n\t -m: si l'usuari vol més informacio sobre la matriu de pesos (PWM) i el consens d'aquesta.\n\t -s: si l'usuari vol més informacio sobre la sequencia (nom, contingutC+G, longitud)\n\t -txx: si l'usuari vo fixar un valor llindar en els resultats (si no, per defecte sera 80).\n\n"; 
    } 
    
    if ($options{h} == 3) { 
	print"\nOPTIONS:\n\t -h: if the user wants to receive information about the program.\n\t -v: if the user wants to receive information while the program is running (in English).\n\t -m: if the user needs further information about the PWM (and consensus).\n\t -s: if the user needs further information about the sequence (name, C+G content, length)\n\t -txx: if the user wants to fix a threshold on the results (if not, threshold=80 will be used).\n\n"; 
    }
    exit (1);
} 

######################################################################################
################# CHECKING NECESSARY PARAMETRES INTRODUCTION #########################
######################################################################################

if (scalar(@ARGV)!=2) { 
    
    print_error("Correct USAGE: scaner.pl fileseq filematrix"); 
} 

# And assignation to two variables of the two parametres.

$fitxer=$ARGV[0]; 
$matriu=$ARGV[1]; 

###################################################################################### 
################################ SEQUENCES READING ################################### 
###################################################################################### 

# Opening of the sequence file and puting in a filehandle (and error message if it's impossible).

if (!open(SEQUENCIA,"< $fitxer")) { 
    
    print_error("Eerror while opening  $fitxer"); 
} 

print_mess ("Correct sequence/s opening!");

# Splited  by "#" (it doesn't exist in the sequence file, so we put all the file in $tot). 
# Close of the sequence filehandle.

$/="\#";
$tot = <SEQUENCIA>;
close (SEQUENCIA);

# In the split, divided  by > (every element of @totesseq will be one sequence).

@totesseq= split (/\>/,$tot);

# Delete of the first element (because if the empty word).

shift @totesseq;

# Begining of reading sequence loop. $numero_seq is the counter of sequences. 
# Every time, it takes one element of @toteseq = one sequence.

$numero_seq=0;

while ($numero_seq<scalar(@totesseq)) {
    
    $_=$totesseq[$numero_seq];
    
# Split of the sequence: every line will be one vector element. 
    
    @seq=split(/\n/,$_);
    
    $locus=$seq[0];  # In $locus we save the sequence name, and after it goes to @nomstheseq (to print it easier).  
    $nomsdeseq[$numero_seq]=$locus;
    shift (@seq); # we delete the first element, the name.
    
# Deleting of the \n between lines (format: fasta).
    
    $novaseq=""; 
    $treupunt=0; 
    while ($treupunt <scalar(@seq)) { 
	$novaseq = $novaseq.$seq[$treupunt]; 
	$treupunt++;  
    }  
    
# In every element of @vectordeseq we save the sequencs in fasta format.
    $vectordeseq[$numero_seq]=$novaseq;
    
# We split again in a definitive vector.
    
    @definitivaseq= split(//,$novaseq); 
    
#GC  COUNTER: 
    
    $cmesg=0;
    
    $compcmesg=0;
    
    while ($compcmesg<scalar(@definitivaseq)) {
	if ($definitivaseq[$compcmesg] eq "C" || $definitivaseq[$compcmesg] eq "G") {
	    $cmesg++;
	}
	$compcmesg++;
    }
    
# Printing of information about the sequence (if you have chosen it with -s).
    
    if (exists($options{s})) { 
	print "\n";
	print "\t\t- NAME of the sequence:", $locus,"\n\n"; 
	print "\t\t- LENGTH of the sequence:", scalar(@definitivaseq),"\n\n"; 
	print "\t\t- SEQUENCE: ",@definitivaseq,"\n\n"; 
        print "\t\t- C+G:$cmesg\n\n";
    }
    $numero_seq++;
}

# If you have chose -s, the program finishes here.

if (exists($options{s})) { 
    exit (0);
}

print_mess("Correct sequence/s reading!");

###################################################################################### 
############################# MATRIXES PROCESSING #################################### 
###################################################################################### 

# Opening of the matrixes file and out in a filehandle (and message of error if it's impossible).

if (!open (MATRIU,"<$matriu")) { 
    print_error("Error while opening $matriu");
}

print_mess("Correct matrix/es opening!");

# Initializing the matrix reading: 

# Splited  by "*" (it doesn't exist in the sequence file, so we put all the file in $yo). 

$/="\*";
$yo = <MATRIU>;

# Splited by "#" (every element of the vetor @tu will be one matrix).

@tu =split(/\#/,$yo);

# Delete of the first element of @tu (because it is the empty word).

shift (@tu);

# Begining of reading matrixes loop. $numero_matriu is the counter of matrixes. 
# Every time, it takes one element of @tu = one matrix.

$numero_matriu=0;

while ($numero_matriu < scalar (@tu)) {
    
#we'll save every matrix in $linimatriu.
    
    $liniamatriu=$tu[$numero_matriu];

#Every element of @vectormatriu is one line of the matrix.
    
    @vectormatriu=split(/\n/,$liniamatriu);
    $filamat[$numero_matriu]=scalar(@vectormatriu);
    $nommatriu=$vectormatriu[0]; # We save here the name of the matrix ( and afterwars, it goes to @nomsdematriu.
    $nomsdematriu[$numero_matriu]=$nommatriu;
    
    shift (@vectormatriu); # Delete of the first element, the name.
    
#### Making the PMW #######################################################################
    
    $fila=0;
    
# Printing information about the matrixes (if you have chosen -m).
    
    if (exists($options{m})) { 
	print "\n- NAME OF THE MATRIX:",$nommatriu,"\n";
	print "\nPO\t\tA\t\tC\t\tG\t\tT\t\n"; # Print the first row.
    }
    
    while ($fila < scalar (@vectormatriu)) { 
	$_=$vectormatriu[$fila];
 
	chomp; # Delete \n
	
	$columna=0; 
	@puntuacio=split(/ +/,$_); 
	
	$consens=" ";	
	shift @puntuacio; # Delete the first element of every line: position. 
	$consens=$puntuacio[scalar(@puntuacio)-1];
	pop @puntuacio; 
	
# calculus of each PWM position: 
	# we add the total value of the row to check the number of nucleotides on each row.
        
	
	$sum=0; 
	
	while ($columna < scalar(@puntuacio)) { 
	    $sum=$sum+$puntuacio[$columna]; 
	    $columna++; 
	} 
	
	# normalise each postion of the PWM :calculating the quocient of the initial
	# value divided by the total of the row, and by 1/4 (a priori probability) 
	# and then making the logarithm.
	# we assign the normalised values to a  new matrix, skipping row 0 by naming row 0 the old row 1.  
	
	$columna=0; 
	
	if (exists($options{m})) {
	    print $fila+1,"\t";
	}
	
	while ($columna < scalar(@puntuacio) ) { 
	    $puntuacio[$columna]=($puntuacio[$columna]+$ZERO)/$sum;
	    $puntuacio[$columna]=$puntuacio[$columna]*4;
	    $puntuacio[$columna]=log($puntuacio[$columna]);
	    $score[$numero_matriu][$fila][$columna]=$puntuacio[$columna]; 
	    
	    if (exists($options{m})) { 
		printf "\t%.2f\t", $score[$numero_matriu][$fila][$columna]; 
	    } 
	    $columna++; 
	} 
	
	# printing information about the matrix...
	
	if (exists($options{m})) { 
	    print "\t",$consens;
	    print "\n"; 
	} 
	
	$fila++; 
    } 
    
    if (exists($options{m})) { 
	print "\nNumber of rows(filas /files/reihen):", $fila,"\n\n"; 
    } 

	$filera[$numero_matriu]=$fila;

    $numero_matriu++;
}

if (exists($options{m})) {
    exit (0);
}



#closing matrix filehandle

close (MATRIU); 

print_mess("Correct matrix/es reading!");


######################################################################################### 
########################## SEQUENCE CANDIDATES SCORING ################################## 
######################################################################################### 
#
# we go through each sequence position up to the maximum matrix value and compare it with the matrix
# we add each position value to the former one.
#

print_mess("Scoring candidates of the sequences..."); 
print_mess("Initializing candidates filtering..."); 

# loop 1: movement through the sequences file, one sequence per loop.

$numero_seq=0;
$comptapos=0;
while ($numero_seq < scalar (@vectordeseq)) {
    print "\n\t\tSEQ.: $nomsdeseq[$numero_seq]\n";
    print OUTFILE "<br><font color='8B0000'> <font size='6'><align='center'><blockquote><blockquote><b>SEQ.: $nomsdeseq[$numero_seq]</font> </font></blockquote></blockquote>\n";
    $seqactual=$vectordeseq[$numero_seq];
    @definitivaseq=split(//,$seqactual);
    $num_matriu_actual=0;
    
    # loop 2: movement through the matrixes  file, one matrix  per loop.
    
    while ($num_matriu_actual < $numero_matriu) {
	print "\n\t\t\tMATRIX:",$nomsdematriu[$num_matriu_actual], "\n";
	print OUTFILE "<br><blockquote><blockquote><blockquote><font size='5'>MATRIX: $nomsdematriu[$num_matriu_actual]</blockquote></blockquote></blockquote><br><center><table border=0 cellpadding=10 cellspacing=10><tr><th bgcolor='8B0000'>BEGINNING</th><th bgcolor='8B0000'>ENDING</th><th bgcolor='8B0000'>SEQUENCE</th><th bgcolor='8B0000'>SCORE</th><tr>\n";
        $posseq=0; 
       	$tamany_matriu=0; 
	
        # loop 3 : movement inside each seq	
	
	while ($posseq<(scalar(@definitivaseq)-$filera[$num_matriu_actual])) { 
	    $posmatriu=0; 
	    $puntuacio_final=0; 
	    
	    # loop 4: movement inside each matrix: score of each position depending on the existing base 
	    
	    $tamany_matriu= $posseq; 
	    
	    while ($posmatriu < $filera[$num_matriu_actual]) { 
		if ($definitivaseq[$tamany_matriu] eq "A"){ 
		    $puntuacio_final=$puntuacio_final+$score[$num_matriu_actual][$posmatriu][$A]; 
		} 
		if ($definitivaseq[$tamany_matriu] eq "C"){ 
		    $puntuacio_final=$puntuacio_final+$score[$num_matriu_actual][$posmatriu][$C]; 
		} 
		if ($definitivaseq[$tamany_matriu] eq "G"){ 
		    $puntuacio_final=$puntuacio_final+$score[$num_matriu_actual][$posmatriu][$G]; 
		} 
		if ($definitivaseq[$tamany_matriu] eq "T"){ 
		    $puntuacio_final=$puntuacio_final+$score[$num_matriu_actual][$posmatriu][$T]; 
		} 
		$posmatriu++; 
		$tamany_matriu++; 
	    } 
	    
	    $llindar=7; # threshold by default
	    $lletra=0; 
	    
	    # assignment of the chosen threshold if option t is chosen.
	  
	    if (exists($options{t})) { 
		$llindar=$options{t}; 
	    } 

            # filtering of the candidates according to the threshold.

	    
            ################################################################################## 
            ################ CANDIDATES FILTERING  ########################################### 
            ################################################################################## 
        
	    if ($puntuacio_final>$llindar) { 
		$lletra=$posseq; 
		$posicions[$comptapos]=$lletra+1;
		$poslletra=$posseq; 
		print "\t\t\tPOSITION:$lletra "; 
		print OUTFILE "<td>$lletra</td>";
		$tot=$lletra+$filera[$num_matriu_actual];

		print OUTFILE "<td>$tot</td>\n";
		print OUTFILE "<td>\n";
	    
		# loop to print the sequences candidates' filtering.
                # If the score is greater than the threshold the seq is printed
                
		while ($poslletra< $lletra+$filera[$num_matriu_actual]) { 
		    print $definitivaseq[$poslletra]; 
		    print OUTFILE "$definitivaseq[$poslletra]"; 
		    $poslletra++; 
		} 
		print OUTFILE "</td>";
		printf  ("\tSCORE:%.5f\n",$puntuacio_final); 
		print OUTFILE "<td>$puntuacio_final</td>\n";
	    } 
	    else 		{
		$posicions[$comptapos]=0;
	    }
	    
	    if ($puntuacio_final>$llindar) {
		print OUTFILE ("<tr>\n");
	    }
	    $posseq++;
	    $comptapos++;
	}
	print OUTFILE "</table></center><hr>\n";

	$num_matriu_actual++;
    } print OUTFILE "<p>\n";
    $numero_seq++;

}
print "\n";


######################################################################################### 
########################## DRAWING A GRAPHICAL PICTURE ################################## 
######################################################################################### 
#
# 
#

$color[0] = "red";
$color[1] = "green";
$color[2] = "orange";
$color[3] = "blue";
$color[4] = "black";
$color[5] = "white";


$dibuseq=0;
$comptapos=0;

print OUTFILE "<font color='8B0000'><font size='6'>VISUAL RESULTS:</font></font>";
while ($dibuseq < scalar(@vectordeseq)){
    $_=$vectordeseq[$dibuseq];
    @pintaseq=split (//,$_);
    $printa=0;
    print OUTFILE "<table border=0><tr><td><font color='8B0000'><font size='6'>$nomsdeseq[$dibuseq] </font></font></td>";
	
    while ($printa< scalar (@pintaseq)){
		print OUTFILE "<td><b>$pintaseq[$printa]</b></td>";
		$printa++;
    } 
    $dibumatrix=0;
    
    while ($dibumatrix<$numero_matriu){
		# $printa counts goes from 0 to SEQLENGTH - MATRIXDIM
		$printa=0;
		
		print OUTFILE "<tr><td><font size='5'>$nomsdematriu[$dibumatrix]</td>";
		
		# traversing a piece of @posicions corresponding to the current matrix
		while ($printa < (scalar(@pintaseq) - $filera[$dibumatrix])) {
			$impressio=0; 

			# a. drawing a candidate
			if ($printa+1 == ($posicions[$comptapos])){
				while ($impressio<$filera[$dibumatrix]) {
					print OUTFILE "<td bgcolor='$color[$dibumatrix]'>&nbsp</td>";
					$impressio++;
				}

				# correction if the counter overlaps the maximum value
				if (($printa + $filera[$dibumatrix]) < (scalar(@pintaseq) - $filera[$dibumatrix]))
				{
					$printa = $printa + $filera[$dibumatrix];
					$comptapos = $comptapos + $filera[$dibumatrix];
				}
				else
				{
					$printa++;
					$comptapos++;
				}
			}   
			# b. drawing an empty cell
			else { 
				print OUTFILE "<td>&nbsp</td>";
				$printa++;
				$comptapos++;
			}
		}

		$dibumatrix++;
	}
    print OUTFILE "</table>";
    $dibuseq++;
}



print OUTFILE "<br><font color='8B0000'><br><hr>
<br><font color='8B0000'><font size='4'>By:   </font>
<font color='483D8'><font size='4'> Guiomar Solanas & Alex Vendrell. .</font></font> <br><font size='4'>For further information...  </font></font><font size='4'><A HREF='./mipagweb.html'>Visit the SCANNER homepage!</a>\n<br><font color='483D8'><br><center>Universitat Pompeu Fabra, March 2002</center></font>";

close (OUTFILE);

print_mess("Finished scoring!"); 


print_mess("Candidates filtering finished!"); 


print_mess("Program successfully finished!\n\n"); 


exit(0); 

 
