#!/usr/bin/perl -w

use strict;

if (scalar(@ARGV) < 2) {                                             #avisar si no has escrit el nom del fitxers
   print STDERR "T'has oblidat de posar el nom dels fitxers!!!\n";
   exit(1);
}

my %anotacions = llegir_anotacions($ARGV[0]);                        #funció llegir anotacions refGene



my $nomfitxersnps = $ARGV[1];


if (!open (SNPS,"< $nomfitxersnps")) {                               #avisar si no s'obre el fitxer dels SNPs
    print STDERR "Impossible obrir $nomfitxersnps\n"; 
    exit(1);
}
   

my $intergenica = 0;
my $comptasnps = 0;
my $contaexons = 0;
my $containtrons = 0;
my $cds = 0;
my $utr3 = 0;
my $utr5 = 0;
my $donor = 0;
my $acceptor = 0;



while (<SNPS>) {                                                     #llegim fitxer SNPs
   my $chro;
   my $inicisnp;


  if ($_ =~ m/\A\d+\s+(\w+)\s+(\d+)/) {                              #si encaixa la línia amb la següent expressió regular

    $chro = $1;                   #cromosoma on es troba el SNP
    $inicisnp = $2;               #inici SNP
 
   # print "buscant snp de cromosoma $chro a la posicio $inicisnp\n"; #xivato, per controlar funcionament

    my $g;

    if (exists($anotacions{$chro})){
    $g = buscam_un_gen($anotacions{$chro},$inicisnp);  #funció per a buscar SNPs en gens o regions intergèniques
  
    }else{

        $g = -1;

        }
       
  
    if ($g == -1) {  # no l'hem trobat a un gen

     $intergenica = $intergenica + 1 ;

    } elsif ($g != -1) { # l'hem trobat a un gen

       # print "l'snp cau al gen amb identificador $anotacions{$chro}[$g][0]\n"; 

       my @v = busca_exons_introns($anotacions{$chro}[$g],$inicisnp); #funció per veure si el snp es troba en UTR o CDS

       $contaexons = $contaexons + $v[0];
       $containtrons = $containtrons + $v[1];
       $cds = $cds + $v[2];
       $utr3 = $utr3 + $v[3];
       $utr5 = $utr5 + $v[4];
       $donor = $donor + $v[5];
       $acceptor = $acceptor + $v[6];
     
    }

  }
    $comptasnps = $comptasnps + 1;

      if ($comptasnps % 100 == 0){

      print "$comptasnps SNPs analitzats\n"

      }
}


print "Tenim $intergenica SNPs en regions intergeniques\n";
print "Tenim $contaexons SNPs en exons\n";
print "Tenim $containtrons SNPs en introns\n";
print "Tenim $cds SNPs en regions CDs\n";
print "Tenim $utr3 SNPs en regions UTR 3'\n";
print "Tenim $utr5 SNPs en regions UTR 5'\n";  
print "Trobem $donor SNPs en regions donor\n";
print "Trobem $acceptor SNPs en regions acceptor\n";
print "En total hem analitzat $comptasnps SNPs\n";


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

#####################Funció per trobar SNPS en exons i introns#################################################################

sub busca_exons_introns  {
   
     my @gen = @{$_[0]};
     my $pos = $_[1];
     my $contaexons = 0;
     my $containtrons = 0;
     my $z = 0;
     my $cds = 0;
     my $utr3 = 0;
     my $utr5 = 0;
     my $donorp = 0;
     my $donorn = 0;
     my $donor = 0;
     my $acceptorp = 0;
     my $acceptorn= 0;
     my $acceptor = 0;



    while ($z < $gen[6]) {

     if ($pos >= $gen[7][$z][0] && $pos <= $gen[7][$z][1]) {

         $contaexons = $contaexons + 1; #per saber el nombre total d'exons
           
         if ($pos >= $gen[4] && $pos <= $gen[5]){

	     $cds = $cds + 1;
         }

	 if  ($pos < $gen[4]){
	   
	     $utr5 = $utr5 + 1;

         }

         if ($pos > $gen[5]){

	     $utr3 = $utr3 + 1;
	 }
     }


    if ($z > 0 && $pos > $gen[7][$z-1][1] && $pos < $gen[7][$z][0]) {

      $containtrons = $containtrons + 1;

 
      if ($gen[1] eq '+') {

         if ($pos <= ($gen[7][$z-1][1] + 6) && $gen[7][$z][0]-$gen[7][$z-1][1] > 6) {
	   $donorp = $donorp + 1;
          }

     
          if ($pos >= ($gen[7][$z][0] - 20) && $gen[7][$z][0]-$gen[7][$z-1][1] > 20){
	    $acceptorp = $acceptorp + 1;
          }
          
	
      }
 


      if ($gen[1] eq '-') {

        if ($pos >= ($gen[7][$z][0] - 6 && $gen[7][$z][0]-$gen[7][$z-1][1] > 6)){

	    $donorn = $donorn + 1;
        }

        if ($pos <= ($gen[7][$z-1][1] + 20) && $gen[7][$z][0]-$gen[7][$z-1][1] > 20){
	    $acceptorn = $acceptorn + 1;
        }

       
     }

	$donor = $donorp + $donorn;
        $acceptor = $acceptorp + $acceptorn;
     
 
 }
    

      #aqui sabrem si l'snp es troba en un intró!!!!!!!!!!!

   $z = $z + 1;

 }

  
     return ($contaexons,$containtrons,$cds,$utr3,$utr5,$donor,$acceptor);

 }


#####################Funció per trobar SNPs en gens o regions intergeniques#######################################################

sub buscam_un_gen {

  my @dic = @{$_[0]};
  my $pos = $_[1];
  my $min = 0;
  my $max = scalar(@dic)-1;
  my $trobada = 0;
  my $mig;

  while ($min <= $max && !$trobada) {

    $mig = int( ($min+$max)/2 );

    if ($dic[$mig][2] < $pos && $dic[$mig][3] < $pos) {
      $min = $mig + 1;
    }
    else {

      if ($dic[$mig][2] > $pos && $dic[$mig][3] > $pos) {
        $max = $mig - 1;
      }
      else {
        $trobada = 1;
      }

    }

  }

  if ($trobada == 0) {
      $mig = -1;
  }

  return $mig;

}


######################################Funció per llegir anotacions#############################################################

sub llegir_anotacions {

  my $nomfitxeranot = $_[0];
  my %matanotacions;
  my %numgens;

  if (!open(GENOMA,"< $nomfitxeranot")) {                      #avisar si no s'obre el fitxer refGene 
      print STDERR "Impossible obrir $nomfitxeranot\n";
      exit(1);
  }


  while (<GENOMA>) {                  # mentre el fitxer tingui linies per llegir

                                      # si encaixa la linia amb la següent expressió regular
    if ($_ =~ m/(\w+)\s+(\w+)\s+([\+\-])\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+([^\n]+)/) {

      my $nom  = $1;                  # nom (identificador) del gen
      my $chro = $2;                  # cromosoma on es troba el gen
      my $stra = $3;                  # ''strand'' (+ o -)
      my $init = $4;                  # posicio inicial de la transcripcio
      my $fint = $5;                  # posicio final de la transcripcio
      my $cdsS = $6;                  # posicio inicial de la regio codificant
      my $cdsE = $7;                  # posicio final de la regio codificant
      my $numexo = $8;                # nombre d'exons
      my @exon = ($9 =~ m/(\d+)/g);   # coordenades d'inici i final dels exons

      if (!exists($numgens{$chro})) {
	  $numgens{$chro} = 0;
      }

      my $g = $numgens{$chro};
      my @inifiexo;
      my $j = 0;

      while ($j <= $numexo){

       $inifiexo[$j][0] = $exon[$j];
       $inifiexo[$j][1] = $exon[$j+$numexo];

       $j = $j + 1;

      }
 
      $matanotacions{$chro}[$g][0] = $nom;         #enregistrem les coses
      $matanotacions{$chro}[$g][1] = $stra;
      $matanotacions{$chro}[$g][2] = $init;
      $matanotacions{$chro}[$g][3] = $fint;
      $matanotacions{$chro}[$g][4] = $cdsS;
      $matanotacions{$chro}[$g][5] = $cdsE;
      $matanotacions{$chro}[$g][6] = $numexo;
      $matanotacions{$chro}[$g][7] = [@inifiexo];

      $numgens{$chro} = $numgens{$chro} + 1;


    } else {

      print STDERR "pasapas2.pl: la seguent linia no encaixa\n$_\n";

    }

  }

  my @claus = keys(%matanotacions);

  my  $m =0;
 
  while ($m < scalar(@claus)) {                     #ordenem gens dins cromosoma segons inici transcripció
      
    @{$matanotacions{$claus[$m]}} = sort{ ${$a}[2] <=> ${$b}[2]}(@{$matanotacions{$claus[$m]}});
    $m = $m + 1;
  
  }

  return %matanotacions;
}


