#!/usr/bin/perl -w

use strict;

if (scalar(@ARGV) < 3){
    print "Consens.pl out_del_buscontp proteina_a_comparar nucleotid\n";
    exit(1); 
};

if (!open(LLISTAT,"< $ARGV[0]")) {
    print "Consens.pl: impossible obrir ATPs.txt\n";
    exit(1);
};


my $prota = $ARGV[1];
my $nucleotid = $ARGV[2];

my @posicions;
my $i=0;

while (<LLISTAT>) {
    
    if(m/(\S+)\s+(\S+)\s+(\d+)\s+(\d+)$/o) {
	$posicions[$i][0]= $1;
	$posicions[$i][1]= $2;
	$posicions[$i][2]= $3;
	$posicions[$i][3]= $4;
	$i = $i + 1;
    };
};
close(LLISTAT);

my @entraconsens;
my $r = 0;    
$i = 0;

while ($i < scalar(@posicions)){
    unless ("$posicions[$i][0].ent" eq "$ARGV[1]"){
	
	system ("cp -r /disc9/DB/pdb/$posicions[$i][0].ent.Z .");
	system ("gunzip $posicions[$i][0].ent.Z");
	
# Crido la funcio radiator q em calcula la distancia entre els atoms de tot el pdb i els del nucleotid
# i ens torna una matriu amb unicament els atoms que estan a una distancia de 5 amstrongs.

	
	my $proteina = $posicions[$i][0]; 
	my @residuspropers = &radiator($proteina);
	
    
	my $z = 0;
	while ($z < scalar(@residuspropers)){
	    
	    $entraconsens[$r][0] = $residuspropers[$z][4];
	    $entraconsens[$r][1] = $residuspropers[$z][1];
	    
	    $r = $r + 1;  
	    $z = $z + 1;  
	};
	
	system("rm $posicions[$i][0].ent");
    };
    $i = $i + 1;
    
};




my %proporcions = &proporcions(\@entraconsens);


my @matrix =  &radiatorprot ($prota,$nucleotid);



&comparator(\%proporcions,\@matrix);


##############################################################
##############################################################
################         FUNCIONS          ###################
##############################################################
##############################################################

# NOM: radiator
# PROPOSIT: calcula la distancia entre els atoms de tot el pdb i els del nucleotid
#           
# PARAMETRES: primer -  matriu amb coordenades dels atoms del pdb
#             segon - matriu amb coordenades del nucleotid
#
# RETORNA: i ens torna una matriu amb unicament els atoms que estan a una distancia de 5 amstrongs.


sub radiator() {
   

    my $proteina = $_[0];
    if (!open(PDB,"< $proteina.ent")) {
	print "buscontp.pl: impossible obrir $proteina.ent\n";
	exit(1);
    };
    my @matriuntp;
    my @matriugral;
    
    my $j = 0; # ambas variables enregistran en los respectivos
    my $linea2 = 0; # hashes los vectores correspondientes
    while(<PDB>){
	
	if(m/\S+\s+(\d+)\s+(\S+)\s+([\w+]{3,3})\s+(\S*\s*\d+)\s+[+-]?\d*\.?(\d+)\s+([+-]?\d*\.?\d+)\s+([+-]?\d*\.?\d+)\s+/o){
	  
	    unless ($3 eq "HOH"){
		if ($3 eq $posicions[$i][1]){
		    
		    $matriuntp[0][$j]=$1;
		    $matriuntp[1][$j]=$5;
		    $matriuntp[2][$j]=$6;
		    $matriuntp[3][$j]=$7;
		    $matriuntp[4][$j]=$2;
		    $j = $j + 1;
		    
		}else{
		    
		    $matriugral[0][$linea2]=$posicions[$i][0];
		    $matriugral[1][$linea2]=$1;
		    $matriugral[2][$linea2]=$3;
		    $matriugral[3][$linea2]=$4;
		    $matriugral[4][$linea2]=$5;
		    $matriugral[5][$linea2]=$6;
		    $matriugral[6][$linea2]=$7;
		    $linea2 = $linea2 + 1;
		};
		
	    };
	};
    };
    close(PDB);
    
    
    my $l = 0; 
    my $m = 0;
    my $n = 0; 
    my @matriures;

    my $resant = -999;
    my $jate = scalar(@{$matriugral[0]});
    my $tu = scalar(@{$matriuntp[0]});
    
    while ($l < $jate) {
	$m = 0;
	

	while ($m < $tu) { 
	    my $r; 
	    $r = sqrt(($matriugral[4][$l] - $matriuntp[1][$m])**2 + ($matriugral[5][$l] - $matriuntp[2][$m])**2 + ($matriugral[6][$l] - $matriuntp[3][$m])**2);
	    if($r < 5 && $matriugral[3][$l] ne "$resant"){ 
		$matriures[$n][0] = $matriugral[0][$l];
		$matriures[$n][1] = $matriugral[2][$l]; 
		$matriures[$n][2] = $matriugral[3][$l];
		$matriures[$n][3] = $matriugral[1][$l];
		$matriures[$n][4] =  $matriuntp[4][$m];
		$matriures[$n][5] = $r;
		
		$n = $n + 1;
		$resant = $matriugral[3][$l];
	    };
	    $m = $m + 1;
	}; 
	
	$l = $l + 1; 
    };
    return (@matriures);
};#radiator


# NOM: Proporcions
# PROPOSIT: Busca la proporcio de aa per cada un dels atoms del Nucleotid
#           
# PARAMETRES: primer -  matriu amb els atoms del nucleotid i els residus propers
#
# RETORNA: L'output del programa.

sub proporcions() {


    my @entraconsens = @{$_[0]};
    
    if (scalar(@ARGV) < 1){
	print "consens.pl output_del_residuspropers.pl\n";
	exit(1); 
    };
    
    if (!open(LLISTAT,"< $ARGV[0]")) {
	print "consens.pl: impossible obrir output_del_residuspropers.pl\n";
	exit(1);
    };
    
    my $y = 0;
    my %aa;
    my @names;
    my $flag = 0;
    my $reqtvariable;
    my $i;
    $names[0] = "JATE";
    my $l = 0;
    while ($y < scalar(@entraconsens)){
   	$flag = 0;
	$i = 0;
	while ($i < scalar(@names)){
	    if ("$entraconsens[$y][0]" eq "$names[$i]"){
		$flag = 1;
	    };
	    $i = $i + 1;
	};
	
	
	if ($flag == 0){
	    $aa{"$entraconsens[$y][0]"}[0] = $entraconsens[$y][1];
	    $names[$l]="$entraconsens[$y][0]";
	    
	    $l=$l+1;
	}else{
	    $reqtvariable = scalar(@{$aa{"$entraconsens[$y][0]"}});
	    $aa{"$entraconsens[$y][0]"}[$reqtvariable]=$entraconsens[$y][1];
	    
	};
	$y = $y + 1;
    };
 




    my %proporcions;
    my $j=0; #numero de residus trobats per aquest atom
    my %seqcons;
    my $k;
    $i = 0;
    my @names1;
    

    print "------------------------------------------------------------------------------------------------------------------------------\n";
    print "                                       proporcions de residus per cada atom del nucleotid\n";
    print "------------------------------------------------------------------------------------------------------------------------------\n";
    print "       | ALA | ARG | ASN | ASP | CYS | GLN | GLU | GLY | HIS | ILE | LEU | LYS | MET | PHE | PRO | SER | THR | TRP | TYR | VAL |\n";
    print "------------------------------------------------------------------------------------------------------------------------------\n";

    while ($i<scalar(@names)){
	$j = 0;
	$l = 0;
	@names1 = (0);
	
	while ($j<scalar(@{$aa{"$names[$i]"}})){
	    
	    $flag = 0;
	    $k = 0;
	    
	    while($k<$j){
		if ($aa{"$names[$i]"}[$k] eq $aa{"$names[$i]"}[$j]){
		    $flag = 1;
		};
		$k = $k + 1;
	    };

	    if ($flag == 0){
		$seqcons{"$aa{$names[$i]}[$j]"} = 1;
		$names1[$l]=$aa{$names[$i]}[$j];
		$l = $l + 1;
	   
	    }else{
		$seqcons{"$aa{$names[$i]}[$j]"} = $seqcons{"$aa{$names[$i]}[$j]"} + 1;
		
	    };	   
	    $j = $j + 1;	    
	};
	my %taula;
	my $m = 0;
	my $proporcio = "0.00";
	
	$taula{ALA} = "0.00";
	$taula{ARG} = "0.00";
	$taula{ASN} = "0.00";
	$taula{ASP} = "0.00";
	$taula{CYS} = "0.00";
	$taula{GLN} = "0.00";
	$taula{GLU} = "0.00";
	$taula{GLY} = "0.00";
	$taula{HIS} = "0.00";
	$taula{ILE} = "0.00";
	$taula{LEU} = "0.00";
	$taula{LYS} = "0.00";
	$taula{MET} = "0.00";
	$taula{PHE} = "0.00";
	$taula{PRO} = "0.00";
	$taula{SER} = "0.00";
	$taula{THR} = "0.00";
	$taula{TRP} = "0.00";
	$taula{TYR} = "0.00";
	$taula{VAL} = "0.00";
	
	while ($m < $l){
	    my @socorrido = 0;
	    my $arrodonint = 0;

	    $proporcio = $seqcons{"$names1[$m]"}/$j;

	    @socorrido = split(//og,$proporcio);

	    if (scalar(@socorrido)>4){
		$arrodonint = $socorrido[4];
	    };

	    $proporcio = substr($proporcio, 0, 4);

	    if ($arrodonint > 4){
		$proporcio = $proporcio + 0.01;
	    };

	    
	    if ("$names1[$m]" eq "ALA"){
		$taula{ALA} = $proporcio;
	    };
	    if ("$names1[$m]" eq "ARG"){
		$taula{ARG} = $proporcio;
	    };
	    if ("$names1[$m]" eq "ASN"){
		$taula{ASN} = $proporcio;
	    };
	    if ("$names1[$m]" eq "ASP"){
		$taula{ASP} = $proporcio;
	    };
	    if ("$names1[$m]" eq "CYS"){
		$taula{CYS} = $proporcio;
	    };
	    if ("$names1[$m]" eq "GLN"){
		$taula{GLN} = $proporcio;
	    };
	    if ("$names1[$m]" eq "GLU"){
		$taula{GLU} = $proporcio;
	    };
	    if ("$names1[$m]" eq "GLY"){
		$taula{ARG} = $proporcio;
	    };
	    if ("$names1[$m]" eq "HIS"){
		$taula{HIS} = $proporcio;
	    };
	    if ("$names1[$m]" eq "ILE"){
		$taula{ILE} = $proporcio;
	    };
	    if ("$names1[$m]" eq "LEU"){
		$taula{LEU} = $proporcio;
	    };
	    if ("$names1[$m]" eq "LYS"){
		$taula{LYS} = $proporcio;
	    };
	    if ("$names1[$m]" eq "MET"){
		$taula{MET} = $proporcio;
	    };
	    if ("$names1[$m]" eq "PHE"){
		$taula{PHE} = $proporcio;
	    };
	    if ("$names1[$m]" eq "PRO"){
		$taula{PRO} = $proporcio;
	    };
	    if ("$names1[$m]" eq "SER"){
		$taula{SER} = $proporcio;
	    };
	    if ("$names1[$m]" eq "THR"){
		$taula{THR} = $proporcio;
	    };
	    if ("$names1[$m]" eq "TRP"){
		$taula{TRP} = $proporcio;
	    };
	    if ("$names1[$m]" eq "TYR"){
		$taula{TYR} = $proporcio;
	    };
	    if ("$names1[$m]" eq "VAL"){
		$taula{VAL} = $proporcio;
	    };
	    	    
	    $m = $m + 1;
	};
	printf " %-5s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s  %4s \n",
	$names[$i],$taula{ALA},$taula{ARG},$taula{ASN},$taula{ASP},$taula{CYS},$taula{GLN},$taula{GLU},$taula{GLY},$taula{HIS},$taula{ILE},$taula{LEU},$taula{LYS},$taula{MET},$taula{PHE},$taula{PRO},$taula{SER},$taula{THR},$taula{TRP},$taula{TYR},$taula{VAL};
	
	
	$proporcions{$names[$i]}[0] = $taula{ALA};
	$proporcions{$names[$i]}[1] = $taula{ARG};
	$proporcions{$names[$i]}[2] = $taula{ASN};
	$proporcions{$names[$i]}[3] = $taula{ASP};
	$proporcions{$names[$i]}[4] = $taula{CYS};
	$proporcions{$names[$i]}[5] = $taula{GLN};
	$proporcions{$names[$i]}[6] = $taula{GLU};
	$proporcions{$names[$i]}[7] = $taula{GLY};
	$proporcions{$names[$i]}[8] = $taula{HIS};
	$proporcions{$names[$i]}[9] = $taula{ILE};
	$proporcions{$names[$i]}[10] = $taula{LEU};
	$proporcions{$names[$i]}[11] = $taula{LYS};
	$proporcions{$names[$i]}[12] = $taula{MET};
	$proporcions{$names[$i]}[13] = $taula{PHE};
	$proporcions{$names[$i]}[14] = $taula{PRO};
	$proporcions{$names[$i]}[15] = $taula{SER};
	$proporcions{$names[$i]}[16] = $taula{THR};
	$proporcions{$names[$i]}[17] = $taula{TRP};
	$proporcions{$names[$i]}[18] = $taula{TYR};
	$proporcions{$names[$i]}[19] = $taula{VAL};

	
	$i = $i + 1;
    };
    return (%proporcions);

};#proporcions







# NOM: radiatorprot
# PROPOSIT: calcula el mateix que la funcio radiator pero per la proteina problema
#           
# PARAMETRES: primer -  vector en el que la primera posicio es la inicial del NT i la segona la final
#
# RETORNA: dona un output en el que veiem els residus que estan aprop dels atoms del nucleotid.


sub radiatorprot (){

       
    my $nucleotid = $_[1];
    my $proteina = $_[0];

    if (!open(PDB,"< $proteina")) {
	print "buscontp.pl: impossible obrir $proteina\n";
	exit(1);
    };

    my @matriuntp;
    my @matriugral;
    
    my $j = 0; # ambas variables enregistran en los respectivos
    my $linea2 = 0; # hashes los vectores correspondientes

    while(<PDB>){
	
	if(m/\S+\s+(\d+)\s+(\S+)\s+([\w+]{3,3})\s+(\S*\s*\d+)\s+[+-]?\d*\.?(\d+)\s+([+-]?\d*\.?\d+)\s+([+-]?\d*\.?\d+)\s+/o){
	  
	    unless ($3 eq "HOH"){
		if ($3 eq "$nucleotid"){
 
		    $matriuntp[0][$j]=$1;
		    $matriuntp[1][$j]=$5;
		    $matriuntp[2][$j]=$6;
		    $matriuntp[3][$j]=$7;
		    $matriuntp[4][$j]=$2;
		    $j = $j + 1;
		    
		}else{
	
		    $matriugral[0][$linea2]=$proteina;
		    $matriugral[1][$linea2]=$1;
		    $matriugral[2][$linea2]=$3;
		    $matriugral[3][$linea2]=$4;
		    $matriugral[4][$linea2]=$5;
		    $matriugral[5][$linea2]=$6;
		    $matriugral[6][$linea2]=$7;
		    $linea2 = $linea2 + 1;
		};
	    };
	};
    };
    close(PDB);
    
    
    my $l = 0; 
    my $m = 0;
    my $n = 0; 
    my @matriures;

    my $resant = -999;
    my $jate = scalar(@{$matriugral[0]});
    my $tu = scalar(@{$matriuntp[0]});
    
print "La proteina $proteina te aquests residus a una distancia maxima de 5 A dels seus atoms:\n";
		print"------------------------------------------\n";
		print " atom    residu    distancia(A)\n";
		print"------------------------------------------\n";

    while ($l < $jate) {
	$m = 0;
	

	while ($m < $tu) { 
	    my $r; 
	    $r = sqrt(($matriugral[4][$l] - $matriuntp[1][$m])**2 + ($matriugral[5][$l] - $matriuntp[2][$m])**2 + ($matriugral[6][$l] - $matriuntp[3][$m])**2);
	   
	    if($r < 5 && $matriugral[3][$l] ne "$resant"){ 
		$matriures[$n][0] = $matriugral[0][$l];
		$matriures[$n][1] = $matriugral[2][$l]; 
		$matriures[$n][2] = $matriugral[3][$l];
		$matriures[$n][3] = $matriugral[1][$l];
		$matriures[$n][4] =  $matriuntp[4][$m];
		$matriures[$n][5] = $r;
		
		print " $matriures[$n][4]\t";
		print " $matriures[$n][1]\t"; 
		print "  $matriures[$n][5]\n";


		$n = $n + 1;
		$resant = $matriugral[3][$l];
	    };
	    $m = $m + 1;
	}; 
	$l = $l + 1; 
    };
    return (@matriures);
}; #radiatorprot










# NOM: comparator
# PROPOSIT: calcula un numero del 0 a l'1 que ens diu lo frequents que son els residus a una distancia de 5 A
#          dels atoms del nucleotid de la proteina problema respecte a les proporcions trobades anteriorment.
#           
# PARAMETRES: primer -  
#
# RETORNA: dona un output en el que veiem un numero del 0 a l'1 que ens diu lo frequents que son els residus a una distancia de 5 A
#          dels atoms del nucleotid de la proteina problema respecte a les proporcions trobades anteriorment.


sub comparator (){
    
    my %proporcions = %{$_[0]};
    my @matriures = @{$_[1]};
    my $posicio;
    my $n = 0;
    my $suma = 0;
    

    while ($n < scalar(@matriures)){
	
	if ("$matriures[$n][1]" eq "ALA"){
	    $posicio = 0;
	};
	if ("$matriures[$n][1]" eq "ARG"){
	    $posicio = 1;
	};
	if ("$matriures[$n][1]" eq "ASN"){
	    $posicio = 2;
	};
	if ("$matriures[$n][1]" eq "ASP"){
	    $posicio = 3;
	};
	if ("$matriures[$n][1]" eq "CYS"){
	    $posicio = 4;
	};
	if ("$matriures[$n][1]" eq "GLN"){
	    $posicio = 5;
	};
	if ("$matriures[$n][1]" eq "GLU"){
	    $posicio = 6;
	};
	if ("$matriures[$n][1]" eq "GLY"){
	    $posicio = 7;
	};
	if ("$matriures[$n][1]" eq "HIS"){
	    $posicio = 8;
	};
	if ("$matriures[$n][1]" eq "ILE"){
	    $posicio = 9;
	};
	if ("$matriures[$n][1]" eq "LEU"){
	    $posicio = 10;
	};
	if ("$matriures[$n][1]" eq "LYS"){
	    $posicio = 11;
	};
	if ("$matriures[$n][1]" eq "MET"){
	    $posicio = 12;
	};
	if ("$matriures[$n][1]" eq "PHE"){
	    $posicio = 13;
	};
	if ("$matriures[$n][1]" eq "PRO"){
	    $posicio = 14;
	};
	if ("$matriures[$n][1]" eq "SER"){
	    $posicio = 15;
	};
	if ("$matriures[$n][1]" eq "THR"){
	    $posicio = 16;
	};
	if ("$matriures[$n][1]" eq "TRP"){
	    $posicio = 17;
	};
	if ("$matriures[$n][1]" eq "TYR"){
	    $posicio = 18;
	};
	if ("$matriures[$n][1]" eq "VAL"){
	    $posicio = 19;
	};
	unless ($proporcions{$matriures[$n][4]}[$posicio] == 0){
	    $suma = $suma - log($proporcions{$matriures[$n][4]}[$posicio]);
	
	};
	$n = $n +1;
    };
    print "la probabilitat de trobar els residus concrets d'aquesta proteina aprop dels atoms del nucleotid es: 10 exp -$suma !!!!!!!!!";
};