#!/usr/bin/perl -w

use strict;


my @genoma;
my $n = 1; #numero de cromosoma
my $i; #lectura de gens

#Enregistrem els fitxers dels cromosomes en una matriu, que despres utilitzarem per esbrinar on cau cada SNP.
while ($n <= 24){
    my $nomfitxer ="chr$n.txt";
    open (CROMOSOMA,"<$nomfitxer"); #Obrim cada fitxer.
    $i = 0; #per anar llegint gens 
    while (<CROMOSOMA>){
	 
	$_=~m/(\w+)\s+(\w+)\s+([\+\-])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\,]+\s+[\d\,]+)/o;
	
#encaixar l'expressio regular de definir cada fila del fitxer de cada cromosoma.
	my $chr = $1; 	#numero de cromosoma
	my $iden = $2;	#identificador de RefSeq
	my $strand = $3;	#strand + o -
	my $itrans = $4;	#posicio d'inici de transcripcio
	my $ftrans = $5;	#posicio de final de transcripcio
	my $itrad = $6;	#posicio d'inici de traduccio
	my $ftrad = $7;	#posicio de final de traduccio
	my $nexons = $8;	#numero d'exons que te el gen
	my $ife2 = $9; 	#posem en un vector inici i final de cada exo. La primera meitat son inicis i la segona meitat son finals.
	my @ife = split(/,/,$ife2);	#Li diem que cada posicio del vector contindra els numeros d'inici i final dels exons, que es troben separats per comes dins el fitxer. Sino posaria tots els inicis o tots els finals en una sola posicio del vector.
	
#definim posicions dins la matriu genoma

	$genoma[$n][$i][0]= $iden;
	$genoma[$n][$i][1]= $strand;
	$genoma[$n][$i][2]= $itrans;
	$genoma[$n][$i][3]= $ftrans;
	$genoma[$n][$i][4]= $itrad;
	$genoma[$n][$i][5]= $ftrad;
	$genoma[$n][$i][6]= $nexons;
	$genoma[$n][$i][7]= \@ife;

	$i = $i + 1; #per anar llegint gens de cada cromosoma
	
    }
    	
    close (CROMOSOMA);	#Tanquem el fitxer.

    @{$genoma[$n]} = sort{ ${$a}[2] <=> ${$b}[2] }(@{$genoma[$n]}); #Ordenem els gens per numero d'identificador RefSeq un cop estan enregistrats dins la matriu.

    $n = $n + 1; #serveix per entrar a cada cromosoma
}

#print "hem llegit el genoma\n";


#anem a llegir fitxer snps 
my $fitxersnps;
my $j = 0; #per anar llegint files del fitxer de SNPs.Passarem a una altre linia quan ho haguem comparat amb el genoma.
my $intro = 0; #per calcular el total d'SNPS que cauen en introns.
my $donor = 0; #per calcular el total d'SNPS que cauen en donors.
my $acceptor = 0; #per calcular el total d'SNPS que cauen en acceptors.
my $cds = 0; #per calcular el total d'SNPS que cauen en CDS. 
my $utr5 = 0; #per calcular el total d'SNPS que cauen en 5'UTR.
my $utr3 = 0; #per calcular el total d'SNPS que cauen en 3'UTR. 
my $intergenic = 0; #SNPs que cauen fora de gens.
my $comptat = 0; #Si hem comptat el SNP sortim.


open (SNP,"<$ARGV[0]"); #Obrirem el fitxer que nosaltres cridem a l'hora d'executar el programa, perque el fitxer de SNPs es tan gran que hi hem d'accedir des de la seva ubicacio inicial.

while (<SNP>){
  
    my @v = split(/\s+/,$_); #Posem en cada posicio del vector cadascuna de les columnes del fitxer, separades per un o mes espais.

    my $chrsnp = $v[1]; #Columna on hi ha el cromosoma on cau el SNP.
    my $pos = $v[3];	#Posicio del SNP dins el genoma.
    my $str = $v[6];	#Strand on cau, + o -.
      
    $chrsnp =~ s/^chr//oi; #Agafem nomes el numero del cromosoma on cau el SNP i l'enregistrem en una variable.



    if ($chrsnp eq "X"){             #ho fem perque el fitxer SNPs no te annotacions 23 i 24. Necessitem tenim numerats igual els dos fitxers. Cada numero de cromosoma quedara enregistrat en la variable $n, com en @genoma.
	$n = 23;
    }else {
	if ($chrsnp eq "Y"){
	    $n = 24;
	}else {
	    $n = $chrsnp;
        }
    }
    
#Donarem valor 1 a $n i quan li sumem 1 anirem canviant de cromosoma. Com que $n significa el mateix en @genoma que en el fitxer de SNPs, anirem a buscar on cau el SNP directament al cromosoma on esta.
    
    if ($n =~ /\A\d+\Z/) {	#Amb aquesta expressio regular descartem aquells SNPs que no sabem en quin cromosoma cauen, i que per tant no hi ha cap numero despres de chr.
    
      my @numcromosoma = @{$genoma[$n]};    #Ara treballem amb matrius de 2D, per cromosoma.
      my $indexgens = &buscagen(\@numcromosoma,$pos,$str); #Cridem la funcio buscagens.Ens donara la posicio del gen dins la matriu @numcromosoma.
  
    
      $i = $indexgens; #Fila del gen on cau el SNP. Ho retorna la funcio buscagens.
    
      if ($i == -1) { #Si no hem trobat cap gen significa que el SNP es inntergenic.
	$intergenic = $intergenic +1;
	$comptat = 1
	
      } elsif ($i != -2) {

         my @infiex = @{$genoma[$n][$i][7]}; #Desreferenciem el vector de inici i final d'exons per poder-lo utilitzar
    
         #Aquestes variables ens serviran per recorrer el vector @infiex.
         my $in = 0; #La primera posicio del vector sera  l'inici del primer exo.
         my $fin = $genoma[$n][$i][6]; #El final del primer exo coincideix amb la posicio corresponent al num d'exons, perque tenim igual num d'inici d'exons que num d'exons.

         #Mentre no sabem on cau el SNP, entrem dins d'aquest while i anem definint cadascuna de les regions del gen i anem mirant si el SNP cau en alguna d'aquestes regions.
#En el moment que trobem on cau, $comptat=1 i sortim del bucle.
         
         $comptat = 0;    #Per saber si hem trobat on cau el SNP.
         while($in < $genoma[$n][$i][6] && $comptat == 0){

	    
		  if ( $pos > ($infiex[$fin] + 36)  && $pos < ($infiex[$in+1] - 3)) {
		      $intro = $intro + 1;
		      $comptat = 1;
		  }
	    
		  if ($pos >= ($infiex[$fin] - 3) && $pos <= ($infiex[$fin] + 36)){
		
		      $donor = $donor + 1;
		      $comptat = 1;
		  }
	    
		  if ($pos >= ($infiex[$in] - 20) && $pos <= ($infiex[$in] + 3)){
		      $acceptor = $acceptor + 1;
		      $comptat = 1;
		  }
	    
		  if ($pos > ($infiex[$in] + 3) && $pos < ($infiex[$fin] - 3) && $comptat == 0) {
		      
		      if ($pos >= $numcromosoma[$i][4] && $pos <= $numcromosoma[$i][5]){
			  $cds = $cds + 1;
			  $comptat = 1;
		      }    
		
		      if ($pos >= $numcromosoma[$i][2] && $pos < $numcromosoma[$i][4]){
			  $utr5 = $utr5 + 1;
			  $comptat = 1;
		      }
		
		      if ($pos > $numcromosoma[$i][5] && $pos <= $numcromosoma[$i][3]){
			  $utr3 = $utr3 + 1;
			  $comptat =  1;
		      }
		
		
		  }
		#Per mourens dins de les posicions d'inici i final d'exo.
		  $fin = $fin + 1; 
		  $in = $in + 1;
	      }
	
        }
   
        $j = $j + 1;	#Canviem de SNP.

        if ($j % 100000 == 0) {	#Per tenir un control del programa. Cada 100000 SNPs que llegeix ho printa per pantalla i aixi podem saber que el programa va llegint.
	
        }
    }
}

close (SNP);



# anem a calcular el  % d'snps que cauen a cada regio dels gens.
my $snpstotals = 9156175; # wc del fitxer d'snps   
my $intergenicstotals = ($intergenic / $snpstotals) * 100;
my $intronstotals = ($intro / $snpstotals) * 100;
my $donorstotals = ($donor / $snpstotals) * 100;
my $acceptorstotals = ($acceptor / $snpstotals) * 100;
my $cdstotals = ($cds / $snpstotals) * 100;
my $utr5totals = ($utr5 / $snpstotals) * 100;
my $utr3totals = ($utr3 / $snpstotals) * 100;

print " el % de snps que cau en INTRONS es $intronstotals\n";
print " el % de snps que cau en DONORS es $donorstotals\n";
print " el % de snps que cau en ACCEPTORS es $acceptorstotals\n";
print " el % de snps que cau en CDS es $cdstotals\n";
print " el % de snps que cau en 5'UTR es $utr5totals\n";
print " el % de snps que cau en 3'UTR es $utr3totals\n";
print " el % de snps que cau en REGIO INTERGENICA es $intergenicstotals\n";
 
#########################################################################################################################################
####################################################   FUNCIONS    ###########################################################################################################################################################################################################################################################################################################################################################
#Funcio per buscar el gen que conte el SNP dins la matriu @numcromosoma. I ens retornara la fila del el gen a $indexgens. 

sub buscagen {

    my $intergenic = 0;	#Si el SNP no cau en cap gen.
    my $indexgens = 0;
    my @numcromosoma = @{$_[0]};    #Ara treballem amb matrius de 2D, per cromosoma.
    my $max = scalar(@numcromosoma) - 1;
    my $min = 0;
    my $mig;
    my $gentrobat = 0;
    my $snp = $_[1];
    my $str = $_[2];
   
    while ($min <= $max && !$gentrobat) {
	$mig = int ( ($min + $max)/2);
        
	if ($snp <  $numcromosoma[$mig][2] && $snp <  $numcromosoma[$mig][3]) {
	    
	    $max = $mig - 1;
	}
	else {
	    
	    if ($snp > $numcromosoma[$mig][2] && $snp > $numcromosoma[$mig][3]) {
		
		$min = $mig + 1;

	    }  else {

		if ($str eq $numcromosoma[$mig][1]) { #Nomes ens interessa si el SNP i el gen que hem trobat tenem la mateixa strand. Sino el descartem.

		    $gentrobat = 1;
		    $indexgens = $mig;
		} else {
		    $gentrobat = 1;
		    $indexgens = -2;
		}
	    }
        }  
    }
    
    if ($gentrobat == 0){
        $indexgens = -1;
    
    }
    
    return $indexgens;
}