#!/usr/bin/perl -w
#
#  
#   LookUpstream.pl
#
#
#  Eulàlia Salichs i Carme Gubern, Febrer-Març 2002;
#
##################################################################
##################################################################

use strict;

# Emprem aquest modul per utilitzar les funcions getopts: 

use Getopt::Std;

# Declarem totes les constants:

my $ZERO = 0.0000000000000000000000001;
my $THRESHOLD=-3;

# Declarem totes les variables#############################################

my $fitxer_seq_DNA;
my $fitxer_matriu;
my %matrix; # hash per la matriu de pesos
my %opt; # hash per les opcions

###variables per comptar,avançar###
my $contar;
my $posicio;
my $posicions;
my $posicio_candidat;
my $endavant_candidats;
my $endavant_posicions;
my $recorrer_sequencies;
my $recorrer_matrius;
my $recorrer_candidats;
my $recorrer_posicions;
###################################

my $fila; # lectura del fitxer_matriu
my $linia; # lectura del fitxer_seqDNA
my $long_vector;
my $finestra; # nombre de files de la matriu = nombre de nucleotids de cada candidat
my $pos_final;
my $posicio_seq;
my $score_candidat; # calculada la raó de versemblança
my $nucleotid;
my $suma_total;
my $seq_sencera;
my $nombre_sequencies;
my $nombre_matrius;
my $nombre_candidats;
my $identificador; # de la seqüència
my $nom_matriu;
my $nou_identificador; # nom de la següent seqüència de DNA (no l'actual)
my $nou_nom_matriu;
my $seq_entrada;
my @DNA;
my @mess;
my @sequencies; 
my @identificadors; 
my @matrius;
my @puntuacio; 
my @fila; # introduïm $fila  de la matriu amb que estem treballant
my @candidats; # nombre de candidats per cada sequencia respecte cada matriu
my @posicio_inicial; # de cada candidat respecte cada matriu
my @posicio_final; # de cada candidat respecte cada matriu
my @seq; # conté individualment les sequencies de @sequencies
my $long_seq_sencera;
my $candidat;
my $contingutCG;
my $percentCG;
my $consens;
my $ultim; 
my $llindar; 
my $resultat;
my @colors;

##############################################################
################### PAS 0: INICI #############################
##############################################################

# Declarem i llegim les opcions getopts:

(getopts('vhmst:',\%opt)) or die("Opcio desconeguda\n");


# Assegurem haver cridat tots els fitxers necessaris 
# per executar el programa. Si no és així, ens ho indicarà:

if (scalar(@ARGV) < 2) { 
    print_error("ARGS: treball.pl fitxer_seqDNA fitxer_matriu");
}

# help option
if (exists ($opt{h}))
{
	print("LookUpstream: a program to predict TF binding sites on DNA sequences\n");
	print("\toption -v: verbose\n");
	print("\toption -m: Show information about the matrices\n");
	print("\toption -s: Show information about the sequences\n");
	print("\toption -t x.x: Fix a quality threshold for the results\n");
	exit(0);
}


print_mess ("1. Trobats els fitxers necessaris per l'execucio del programa.");

# Assignem els noms del vector d'arguments (paràmetres) 
# a dues variables amb noms més significatius:

($fitxer_seq_DNA) = shift(@ARGV);
($fitxer_matriu)  = shift(@ARGV);

# Obrim el fitxer on s'emmagatzemaran els resultats per mostrar amb format html:

if (!open(CANDIDAT,">candidats.html")){
    print_error("treball.pl: impossible obrir el fitxer resultat.html");
}

print_mess ("2. Obert correctament el fitxer que contidrà els resultats");

print CANDIDAT"<html>\n";
print CANDIDAT"<head>\n";
print CANDIDAT"<title>\n";
print CANDIDAT"Promotors Potencials\n";
print CANDIDAT"</title>\n";
print CANDIDAT"</head>\n";
print CANDIDAT"<body background=./teal_sand1.gif>\n";
print CANDIDAT"<br>";
print CANDIDAT"<img src=./dnaspin2.gif height=25 width=120 align=right>\n";
print CANDIDAT"<p>";
print CANDIDAT"<img src=./dnaspin2.gif height=25 width=120>\n";
print CANDIDAT"<br><br><br>";
print CANDIDAT"<hr><center><h1><b><font face=Arial,Helvetica size=5><font color=FireBrick><blockquote>PROMOTORS&nbsp;&nbsp;&nbsp;POTENCIALS</h1></center>";
print CANDIDAT"</font></font></b></blockquote><hr>\n";
print CANDIDAT"<br>\n";
print CANDIDAT"<blockquote>\n";

###############################################################
################### PAS 1 :  SEQUENCIES ####################### 
###############################################################

# Obrim el fitxer que conté la seqüència de DNA:

if (!open(SEQUENCIA,"< $fitxer_seq_DNA")) {
    print_error( "treball.pl: impossible obrir $fitxer_seq_DNA");
}

print_mess ("3. Obert correctament el fitxer que conte la sequencia de DNA");

# Llegim la primera línia del fitxer_seq_DNA:

$linia = <SEQUENCIA>;
chomp($linia);
$identificador = $linia;

# Eliminem el símbol ">" característic del format FASTA:

$identificador =~s/>//;    

# Inicialitzem el vectors per emmagatzemar la informació necessaria per la representació gràfica:

@identificadors=();
@sequencies=();
@candidats=();
@posicio_inicial=();
@posicio_final=();

while(){

    print_mess("\n");print_mess ("4.Processant:*************** $identificador *******************");
    
    if (!exists ($opt{s}) && !exists ($opt{m})){ 
		
		print "\n***************** $identificador **********************\n";
		print CANDIDAT"<font face=Courier size=4><font color=FireBrick><b>SEQUENCIA: </b></font></font>\n";
		print CANDIDAT"<font face=Courier size=3><font color=DarkSlateGray>$identificador</font></font><br><p>\n";
    }
    
    $seq_sencera="";
    
    push @identificadors, $identificador; #construïm vector on cada posició contidrà l'identificador de les sequencies
	
    # Llegim la seqüència actual:
    
    while ($linia=<SEQUENCIA>){
		
		last if ($linia =~ />/);
		
		chomp($linia);
		$seq_sencera = $seq_sencera.$linia;
    }
    
    push @sequencies, $seq_sencera; #construïm vector on cada posició contindrà les sequències a analitzar
    
    if (defined($linia)){
		
		# Si queda alguna seqüència  més, processem l'identificador:
		
		chomp($linia);
		$nou_identificador = $linia;
		$nou_identificador =~ s/>//;
		
    }
    
    @DNA = split (//,$seq_sencera);
    $long_seq_sencera = scalar(@DNA);
    
    # Calculem el contingut C+G de la seqüència de DNA:
    
    $posicio=0;
    $contingutCG=0;
    while ($posicio < $long_seq_sencera){
		if (($DNA[$posicio] eq 'C') or ($DNA[$posicio] eq 'G')){
			$contingutCG++;
		}
		$percentCG = ($contingutCG/$long_seq_sencera)*100;
		$posicio++;
    }
    
    # Informació sobre la seqüència:
    
    if (exists($opt{s}) ){
		
		print "\nSequencia: ******$identificador******\n";
		print "Nombre nucleotids: ",$long_seq_sencera,"\n";
		print "Contingut C+G: "; printf("%.2f",$percentCG,); print " % \n";
		
		print_mess("4.1 OPCIO s: seleccionada!");
		print_mess("4.2 OK! Executat correctament");
		
    }
    
    ####################################################
    ################# PAS 1.1 : MATRIU #################
    ####################################################
	
    
    # Obrim la matriu. Si no es possible ens avisara:
    
    if (!open(MATRIU,"< $fitxer_matriu")) {
		print_error("treball.pl: impossible obrir $fitxer_matriu");
    }
    
    print_mess ("5. Obert correctament el fitxer que conte la matriu"); 
    
    
    # Llegim la primera línia de la matriu: 
    
    $fila = <MATRIU>;
    chomp ($fila);
    $nom_matriu = $fila;
    
    # Eliminem el símbol "#" que encapçala cada matriu:
    
    $nom_matriu =~ s/\#//;
 
    @matrius=();
    
    while(){
		
		push @matrius , $nom_matriu; 
#construim vector on  cada posicio contindra el nom de les matrius analitzades
		
		print_mess ("6.Processant: **************  $nom_matriu *************");
		
		if (!exists ($opt{m}) && !exists($opt{s})) { 
			
			print "\n************* $nom_matriu ***************\n";
			print CANDIDAT"<blockquote>\n";
			print CANDIDAT"<font face=Courier size=4><font color=FireBrick><b>MATRIU: </b></font></font>\n";
			print CANDIDAT"<font face=Courier size=3><font color=DarkSlateGray>$nom_matriu</font></font><br><p>\n";
		}
		
        # Llegim la resta de files:
		
		$posicio=0;
		$consens ="";
		
		while($fila=<MATRIU>){
			
			last if ($fila =~ /\#/);
			chomp($fila);
			
			# Calculem el nombre de símbols d'aquesta fila:
			
			@fila = split (/ +/,$fila);
			$ultim = scalar (@fila);
			$consens = $consens.$fila[$ultim-1]; # conté la seqüència consens de la matriu.
			shift(@fila); # eliminem la primera columna del vector @linia
			pop(@fila); # eliminem l'última columna del vector @linia
			$long_vector = scalar (@fila);
			
			# Construïm el hash %matrix:
			
			$matrix{'A'}[$posicio] = $fila[0];
			$matrix{'C'}[$posicio] = $fila[1];
			$matrix{'G'}[$posicio] = $fila[2];
			$matrix{'T'}[$posicio] = $fila[3];
			
			# Calculem el nombre total de seqüències utilitzades:
			
			$suma_total = $matrix{'A'}[$posicio] + $matrix{'C'}[$posicio] + $matrix{'G'}[$posicio] + $matrix{'T'}[$posicio];
			
			# Normalitzarem el hash per trobar el pes de cada nucleòtid (càlcul de la 
			# freqüència relativa de cada nucleòtid respecte el total):
			
			for  $nucleotid ('A','C','G','T'){
				
				if ($matrix{$nucleotid}[$posicio]==0) {
					
					$matrix{$nucleotid}[$posicio]=$ZERO;
				}
				
				$matrix{$nucleotid}[$posicio] = (log($matrix{$nucleotid}[$posicio]) - (log($suma_total))) - (log (0.25));
				
			}
			
			$posicio++;
		}
		
		if (defined($fila)){
			
			# Si queda alguna matriu  més, processem el nom de la matriu:
			
			chomp($fila);
			$nou_nom_matriu = $fila;
			$nou_nom_matriu =~ s/\#//;
			
		}
		
		$finestra=$posicio;
		
		print_mess ("7. Construïda i normalitzada la matriu: $nom_matriu");
		
        # Informació sobre la matriu
		
		if (exists($opt{m}) && !defined($linia) ) { 
			
			print"\n\t*********** $nom_matriu **********\n";
			$posicio=0;
			print "\t","  A\t","  C\t","  G\t","  T\t","\n"; 
			
			while ($posicio<$finestra) {
				print $posicio+1,"\t";
				for $nucleotid ('A','C','G','T'){ 
					
					printf("%.2f\t",$matrix{$nucleotid}[$posicio]);	    
					
				}
				
				print"\n";
				$posicio++;
			}
			
			print "\nS'han utilitzat ", $suma_total, " sequencies de ", $finestra," nucleotids per construir$nom_matriu\n";
			print "Sequencia consens: ", $consens,"\n";
			print_mess("7.1 OPCIO m: seleccionada!");
			print_mess("7.2 OK! Executat correctament");
			
		}
		
		####################################################
        ########### PAS 1.2:AVALUACIO CANDIDATS ############
        ####################################################
		
		print_mess ("8.Processant: ********** $identificador amb $nom_matriu ***********");
		
		
		if (exists($opt{t})){
			
			$llindar=$opt{t}; # threshold introduït per l'ususari
			print_mess("8.1 OPCIO t: el threshold seleccionat es $opt{t}");
			
		}
		
		else {
			
			$llindar=$THRESHOLD; # threshold per defecte
			print_mess("8.1 OPCIO t: el threshold seleccionat per defecte es $THRESHOLD");
		}
		
		if (!exists($opt{s}) and !exists($opt{m})) {  
			
			print CANDIDAT"<table border=0 cellpadding=5 cellspacing=5><tr>\n";
			print CANDIDAT"<th bgcolor=CadetBlue><font face=Courier size=4 color=Beige>Posició Inicial</font></th>\n";
			print CANDIDAT"<th bgcolor=CadetBlue><font face=Courier size=4 color=Beige>Posició Final</font></th>\n";
			print CANDIDAT"<th bgcolor=CadetBlue><font face=Courier size=4 color=Beige>Puntuació</font></th>\n";
			print CANDIDAT"<th bgcolor=CadetBlue><font face=Courier size=4 color=Beige>Candidat</font></th>\n";
			print CANDIDAT"<tr>";
			
			$posicio_seq=0;
			$contar=0; # comptabilitzar el nombre de candidats, que superen el threshold, per cada sequencia analitzada respecte cada matriu
	    
			while ($posicio_seq <= $long_seq_sencera - $finestra) { 
				
				# Avaluem el candidat actual:
				
				$posicio_candidat = 0;
				$score_candidat =0;
				$candidat ="";
				
				# Determinem l'score de cada candidat i la seva seqüència:
				
				while ($posicio_candidat < $finestra) {
					
					$score_candidat = $score_candidat + $matrix{$DNA[$posicio_seq + $posicio_candidat]}[$posicio_candidat];
					$candidat = $candidat.$DNA[$posicio_seq + $posicio_candidat];
					$posicio_candidat++;
					
				}
				
				# Filtrat d'aquells que assoleixen el threshold escollit o l'assignat per defecte:
				
				if ($score_candidat >= $llindar){
					
					print $posicio_seq, "\t", $pos_final=$posicio_seq + $finestra,"\t"; 
					printf("%.2f\t", $score_candidat, "\n");
					print $candidat,"\n";
					
					push @posicio_inicial,$posicio_seq; #emmagatzemem la posició inicial de cada candidat
					push @posicio_final,$pos_final; #emmagatzemem la posició final de cada candidat
					
					$contar++;
					
                    #Perquè en el fitxer de sortida (candidats.html) els resultats es mostrin amb dos decimals: 
					
					@puntuacio = split(//,$score_candidat);
					
					$posicio=0;
					$resultat="";
					while ($puntuacio[$posicio]){
						last if ($puntuacio[$posicio]=~ /\./);
						$resultat=$resultat.$puntuacio[$posicio];
						$posicio++;
					}
					
					$resultat=$resultat.$puntuacio[$posicio];
					$posicio++;
					$resultat=$resultat.$puntuacio[$posicio];
					$posicio++;
					
					if ($puntuacio[$posicio+1]>=5){
						
						$puntuacio[$posicio]=$puntuacio[$posicio]+1;
						$resultat=$resultat.$puntuacio[$posicio];
					}
					
					else {
						
						$resultat=$resultat.$puntuacio[$posicio];
					}
					
					
					
					###############################################################################################
					
					print CANDIDAT"<td align=center><font face=Courier size=3 color=DarkSlateGray>$posicio_seq</font></td>\n";
					print CANDIDAT"<td align=center><font face=Courier size=3 color=DarkSlateGray>$pos_final</font></td>\n";
					print CANDIDAT"<td align=center><font face=Courier size=3 color=DarkSlateGray>$resultat</font></td>\n";
					print CANDIDAT"<td align=center><font face=Courier size=3 color=DarkSlateGray>$candidat</font></td>\n";
					print CANDIDAT"<tr>";
				}
				
				$posicio_seq++;
			}
			
			push @candidats, $contar; #vector que conté el nombre de candidats per cada sequencia, respecte les diferents matrius
			print "\n"; 
			print CANDIDAT "</table>\n"; 
			print CANDIDAT"</blockquote>\n";
			
		}
		
		
		$nom_matriu= $nou_nom_matriu;
		
		
		# Si ja no hi ha més matrius sortirem del bucle:
		
		if (!defined($fila)) {
			last;
		}
		
    }
    
    $identificador = $nou_identificador;
    close(MATRIU);
    print_mess ("9. Tancat el fitxer que conté la matriu");
	
	
    # Si ja no hi ha més seqüències sortirem del bucle:
    
    if (!defined($linia)) {
		last;
    }
    
}

print CANDIDAT"<br><p><br><br>";
print CANDIDAT"</blockquote>\n";


#############################################################################################################################
################################## Representació gràfica dels resultats #####################################################

$colors[0] = "Cyan";
$colors[1] = "Red";
$colors[2] = "Orange";
$colors[3] = "Yellow";
$colors[4] = "Green";


print CANDIDAT "<br><p><br><p>";
print CANDIDAT "<blockquote>\n";
print CANDIDAT"<hr><center><font color=FireBrick><font face= Arial,Helvetica size=5><b> REPRESENTACIO GRAFICA </font></font></center></b><hr><br><p>\n";

$endavant_posicions=0; # recorrer tot el vector de posicions,final i inicial
$endavant_candidats=0; # recorrer tot el vector de candidats
$recorrer_sequencies=0; # recorrer el vector de seqüències

print CANDIDAT"\n<table border =0>\n";

if (!exists($opt{s}) && !exists($opt{m})){
    
    # Per totes les sequencies:
    
    while ($recorrer_sequencies<scalar(@sequencies)){
	
	@seq=split(//,$sequencies[$recorrer_sequencies]);
	
	print CANDIDAT"<th align=right><font color=FireBrick face=Courier size=4>Sequencia";
	print CANDIDAT"<font color=DarkSlateGray face=Courier size=3>:$identificadors[$recorrer_sequencies]&nbsp;</font></font></th>\n";
        
	$posicio=0; # recorrer tota la sequencia actual
	while($posicio<scalar(@seq)){
	    
	    print CANDIDAT "<th><font color=DarkSlateGray face=Courier size=3>$seq[$posicio]</font></th>\n";
	    $posicio++;
	}
	
	print CANDIDAT "<tr><tr>";    
	
        # Per totes les matrius:
	
	$recorrer_matrius=0; # recorrer el vector de matrius
	while ($recorrer_matrius<scalar(@matrius)){
	    
	    print CANDIDAT "<td align=right><font color=DarkSlateGray face=Courier size=4>$matrius[$recorrer_matrius]&nbsp</font></td>\n"; 
	    
	    $recorrer_candidats=0; # execució del bucle tantes vegades com candidats té la sequencia respecte la matriu actual
	    while($recorrer_candidats<$candidats[$endavant_candidats]){
		
		$recorrer_posicions=0; #per contar fins a la posicio inicial i final
		
                # Fins arribar a la posicio inicial 
		while($recorrer_posicions<$posicio_inicial[$endavant_posicions]){
		    
		    print CANDIDAT "<td>&nbsp;</td>\n";
		    $recorrer_posicions++;
		}	    
		
                # Fins arribar a la posicio final
		while($recorrer_posicions>=$posicio_inicial[$endavant_posicions] && $recorrer_posicions<$posicio_final[$endavant_posicions]){
		    
		    print CANDIDAT "<td bgcolor= $colors[$recorrer_matrius] >&nbsp;</td>\n";
		    $recorrer_posicions++;
		}
		
		print CANDIDAT "<tr>";
		print CANDIDAT "<td align=right>&nbsp&nbsp</td>\n";
		$recorrer_candidats++; 
		$endavant_posicions++; # pas a la següent posicio inicial i final del proper candidat
		
	    }
	    
	    $recorrer_matrius++;
	    $endavant_candidats++; 
	    print CANDIDAT "<tr><tr><tr>";
	}
	print CANDIDAT "<tr><tr><tr>";
	$recorrer_sequencies++;
	
    }
    
}

print CANDIDAT "</table>\n";
print CANDIDAT "</blockquote>\n";
print CANDIDAT"<br><br>";

print CANDIDAT"<center>";
print CANDIDAT"<table border=0>";
print CANDIDAT"<blockquote><td align=left><a href=./programa.html#candidats><img src=./larrow.gif  border=0></a><font face=Courier><font color=FireBrick><font size=1><b>Enrera</b></font></font></td>\n";
print CANDIDAT"<td width=600>&nbsp;</td>\n";
print CANDIDAT"<td align= right><font face=Courier><font color=FireBrick><font size=1><b>Inici</b>\n";
print CANDIDAT"<a name=PROMOTORS POTENCIALS></a><a href=# PROMOTORS POTENCIALS><img src=./uarrow.gif  border=0></a></font></font></font></td></blockquote>\n";
print CANDIDAT"</table>";
print CANDIDAT"</center>";
print CANDIDAT"<p><br>";
print CANDIDAT "</body>\n";
print CANDIDAT "</html>\n";

close (CANDIDAT);
close(SEQUENCIA);
print_mess("\n");print_mess("10. Tancat el fitxer que conté la sequencia de DNA");
print_mess ("11. Avaluats tots els candidats");

###########################################################
################## PAS 3 : FINAL ##########################
###########################################################

print_mess("12. OK! Executat correctament");
exit(0);


###########################################################
################## SUBRUTINES #############################
###########################################################

 
# Subrutina per la impressió de missatges d'informació sobre el procés
# a la pantalla mitjançant l'opció -v:

sub print_mess
{
    @mess = @_;
    
    print STDERR ("\t: @mess \n") if (exists($opt{v}));
}

# ERROR: imprimim el missatge i sortim del programa:

sub print_error
{
    my @mess = @_;
    
    print STDERR ("\t: @mess \n") and exit();
}
