#!/usr/bin/perl -w

use strict;


if (scalar(@ARGV) < 1) {
    print "assaig.pl <fitxerdesnps>\n";
    exit(1);
}

my $nomfitxersnps = $ARGV[0];
my $N = 1; #representa el cromosoma
my $i;     #representa el gen
my @gens;

while ($N <= 24) { # relacionem el numero de cromosoma amb N. Comencem amb N=1 i anem augmentant fins arribar a 24, que és el màxim de cromocomes.Ho fem amb tots el cromosomes, inclosos X i Y per tal de poder introduirlos en la matriu d'acord amb el seu valor de N
    my $nomfitxercrom="chr$N.txt";
    if ($N == 23) {
	$nomfitxercrom= "chrX.txt";
    }
    if ($N == 24) {
        $nomfitxercrom= "chrY.txt";
    }
    if (!open(CHR,"< $nomfitxercrom")) {
	print "assaig.pl: impossible obrir fitxer $nomfitxercrom\n";
    }

    $i= 0;
    while (<CHR>) {#obrim el fitxer del cromosoma i l'encaixem en aquesta estructura
      $_ =~ m/(\w)\s+(\w+)\s+([\+\-])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\,]+\s+[\d\,]+)/;
      my $nom = $1;                  # nom (identificador) del gen
      my $chr = $2;                  # cromosoma on es troba el gen
      my $cad = $3;                  # ''strand'' (+ o -)
      my $ini = $4;                  # posicio inicial de la transcripcio
      my $fin = $5;                  # posicio final de la transcripcio
      my $codi = $6;                 # posicio inicial de la regio codificant
      my $codf = $7;                 # posicio final de la regio codificant
      my $exoC = $8;                 # nombre d'exons
      my $resumexons = $9;           # coordenades d'inici i final dels diferents exons
      my @regioexons= ($9 =~ m/(\d+)/g);   # vector: la primera meitat son inicis d'exo i la segona meitat son finals.
    
      $gens[$N][$i][0]= $ini;
      $gens[$N][$i][1]= $fin;
      $gens[$N][$i][2]= $codi;
      $gens[$N][$i][3]= $codf;
      $gens[$N][$i][4]= $exoC;
      $gens[$N][$i][5]= [@regioexons];
      $gens[$N][$i][6]= $cad;
	
      $i = $i + 1; #anem augmentant el valor de $i per poder enregistrar tots els gens d'1 cromosoma (files)
    }
    close(CHR);#tanquem el fitxer del cromosoma quan s'ha acabat de llegir tot, ara el tindrem ja enregistrat en la amtriu (ja tenim una matriu de 2D, quan li afegim un altre cromosoma passara a ser de 3D)
    @{$gens[$N]}= sort{${$a}[0] <=> ${$b}[0]}  (@{$gens[$N]}); # ordenem els inicis dels gens!
    $N = $N + 1;
}
printf "hem llegit els gens\n"; #per a tenir cert control del procés.


if (!open (SNPS,"< $nomfitxersnps")) {#si no s'obre el fitxer SNP avisa'ns!!
    print "assaig.pl: impossible obrir $nomfitxersnps\n";
    exit(1);
}

my $SNPtotal= 0;
my $SNPcodif= 0;
my $SNPenexo= 0;
my $trobat;
my @regioexons;
my $x;
my $y;
my $trobat2;
my $donor=0;
my $acceptor=0;

while (<SNPS>) {#obrim el fitxer de SNPs que hem indicat en la comanda per cridar el programa, i l'encaixem en una estructura
  $_ =~ m/(\d+)\s+(\w+)\s+(\d+)\s+(\d+)\s+(\w+)/;
  my $desc = $1;                  # desconegut
  my $crom = $2;                  # cromosoma on es troba el SNP
  my $pos1 = $3;                  # poscio inicial SNP
  my $pos2 = $4;                  # posicio final SNP
  my $ident = $5;                 # identificdor SNP

  $crom =~ /chr(\w+)/; #perque el programa obri el fitxer del cromosoma requerit cal tornar a dir-li que 23 i 24 són X i Y (en matriu estan com a numeros pero en fitxer de SNPs com a lletres)
  if ($1 eq "X") {
      $N = 23;
  } else {
      if ($1 eq "Y") {
	  $N = 24;
      } else {
	  $N = $1;
      }
  }
  
  if ($N =~ /\A\d+\Z/) {#només quan la $crom sigui un cromosoma, eliminem aquells SNPs que no sabem en quin cromosoma es troben.

    eval{#punt de control (si la funció va bé)
        $i = busca_snp($gens[$N], $pos1);
   };
    if ($i== -1){
	print"SNP fora de gen\n";
    }#per si esta fora
 
    
    @regioexons = @{$gens[$N][$i][5]};
    $x= 0;                     #posicio dels inicis d'exons
    $y= $gens[$N][$i][4];      #posicio dels finals d'exons (la diferencia entre inici i final d'exons ve determinada pel nombre dels mateixos)
    $trobat2 =0; #quan esta en exo
   
    while ($trobat2==0 && $x<$gens[$N][$i][4]) {#vale, aqui posarem per a detectar regions d'spl: incloem regio d'splicing 30 nucleotids per davant i per darrere de l'inici o fi d'exo: entre x-30 i x sera una cosa, i y+30 i y sera el contrari. Creem una variable que sumi cadascun i no cal que surti.mb!
    
	if ($gens[$N][$i][6] eq '+') { #si la cadena es positiva es mira aixo
	         if($pos1 >= $regioexons[$x]-30 && $pos1<= $regioexons[$x]) {
		     $acceptor= $acceptor +1;
		 }
		 if($pos1>= $regioexons[$y] && $pos1 <= $regioexons[$y]+30) {
		     $donor= $donor +1;
	         }
	}
	if ($gens[$N][$i][6] eq '-') { #si la cadena es negativa
		     if($pos1 >= $regioexons[$x]-30 && $pos1<= $regioexons[$x]) {
		     $donor= $donor +1;
		 }
		 if($pos1>= $regioexons[$y] && $pos1 <= $regioexons[$y]+30) {
		     $acceptor= $acceptor +1;
		 }
	}
     
	if($pos1 >= $regioexons[$x] && $pos1 <= $regioexons[$y] ) { #mirem si SNP esta en algun exo, si ho estÃ  q passi al seguent SNP
			$SNPenexo= $SNPenexo + 1;
			$trobat2= 1;
		        		  
		if ($pos1>=$gens[$N][$i][2] && $pos1<=$gens[$N][$i][3]){  #si forma part de la regio codificant (veure etapa en la que s'ha encaixat el fitxer en la matriu)
			    $SNPcodif= $SNPcodif + 1;
 		}
	} 
	else {  #si SNP no esta en l'exo que mirem d'aquest gen, canvia d'exo, augmentant un valor x i y. Aixi anem rastrejant tots els exons possibles
	      $x= $x + 1;
	      $y= $y + 1; 
	}
		   
    }
   
    $SNPtotal= $SNPtotal + 1; #cada vegada que canviem de SNP passem per aqui i sumem 1 a SNPtotal

     if ($SNPtotal % 100000 == 0) { #per anar controlant imprimirem cada vegada que augmenti 100.000
       print "processats $SNPtotal SNPs\n";
     }

  } # fi del if de si el cromosoma es realment un numero 

}

close (SNPS);#tanquem fitxer
my $noexo= $SNPtotal-$SNPenexo;
my $utr= $SNPenexo-$SNPcodif;
my $perc1= ($SNPenexo/$SNPtotal)* 100;
my $perc2= ($noexo/$SNPtotal)* 100;
my $perc3= ($SNPcodif/$SNPtotal)* 100;
my $perc4= ($utr/$SNPtotal)* 100;
my $splicing= $acceptor +$donor;
my $perc5;
my $perc6;
my $perc7;

print "el nombre de SNPs totals analitzats es $SNPtotal\n";
printf "el nombre de SNPs en exons es $SNPenexo i el percentatge es %.2f%%\n",$perc1;
printf "el nombre de SNPs que es troben en intro es $noexo i el percentatge es %.2f%%\n",$perc2;
printf "el nombre de SNPs pertanyents a exons i de regio codificant es $SNPcodif, i per tant el nombre de UTRs es $utr. El percentatge de SNPs d'exons codificants es %.2f%%",$perc3;
printf "i el percentatge de UTRs es %.2f%%\n",$perc4;
print "el nombre de SNPs en splicing es $splicing\n";

#if ($splicing != 0) { #nomes en cas que hi hagi SNPs en slpicing, perque sino no pot dividir per 0 i es torna boig
   $perc5= ($donor/$splicing)* 100;
   $perc6= ($acceptor/$splicing)* 100;
   $perc7= ($splicing/$SNPtotal)* 100;
   printf "que es un %.2f%%\n",$perc7; printf "del total de SNPs (%.2f%%\n",$perc5;printf"en regio donadora i %.2f%%\n",$perc6; printf"en acceptora)\n";

#}

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

# FUNCIO: busca_snp
# PROPOSIT: es busca un snp dins d'un conjunt de gens d'un cromosoma

sub busca_snp {
    my @matgens = @{$_[0]};   #primer element entrat (@gens)
    my $snp = $_[1]; #segon element que crides ($pos1)
    my $min = 0; 
    my $max = scalar(@matgens)-1; #total de gens del cromosoma que mirem
    my $trobada = 0; #variable per quan ha trobat el gen
    my $gen = -1; #variable a retornar: si no canvia el valor voldra dir que no hem trobat gen
    while ($min <= $max && !$trobada) { 
	my $mig = int( ($min+$max)/2 ); #el valor mig de la llista de gens inicialment sera la meitat sense decimals		
	        if ($snp > $matgens[$mig][1]) { #si SNP es superior a la posicio final del gen mig, llavors fes que el minim ara sigui mig +1
		    $min = $mig + 1;
		}
	        else { 
		    if ($snp < $matgens[$mig][0]) { #si SNP es inferior a la posicio d'inici del gen mig
			$max = $mig - 1;
		    }
		    else { #si SNP esta contingut en el gen mig
			$trobada = 1; 
			$gen = $mig;
		    }
		}

        }  
    return $gen;#torna l aposicio del gen, que sera la coordenada $i de la matriu
}
