#!/usr/bin/perl -w

use strict;

use Data::Dumper;
local $Data::Dumper::Purity   = 0;
local $Data::Dumper::Deepcopy = 1;

############################################################################################
############################### DECLARACIO DE VARIABLES ####################################
############################################################################################

my $debug = 0; # si $debug == 1 llavors print Dumper de variables hash...

my %refGene = ();

my %snpNih = ();

my %RESULTAT = ();
	    
if (scalar(@ARGV) != 2) {
	print STDERR "USAGE: programa.pl refseq_file snp_file \n";
	exit(1);
};
my ($reffile,$snpfile) = @ARGV; 

print STDERR "## REFSEQ FILE: $reffile\n";
print STDERR "## SNPs FILE: $snpfile\n";


#############################################################################################
##############################     MAIN LOOP        #########################################
#############################################################################################

&obrir_fitxer_refseq($reffile);
&obrir_fitxer_SNPs($snpfile);
&ordenar_array_SNPs();
&tractament_dades();
&imprimir_resultats();
exit(0);

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

#############################################################################################
##############################OBRIM EL PRIMER FITXER#########################################
#############################################################################################
sub obrir_fitxer_refseq() {
    my $file = $_[0];

	open (FITXER1, "< $file") || do {
		print STDERR "### ERROR!!! NO PUC OBRIR FITXER $file \n";
        exit(1);
    };
	
    print STDERR "##\n## LLEGINT FITXER REFSEQS: $file\n##\n";

	while (<FITXER1>) {

		my ($gid,$chr,$strand,$rnao,$rnae,$cdso,$cdse,$exnum,$exono,$exone);

		next if /^\s*$/o;
		next if /^\#/o;

		chomp($_);
		
		($gid,$chr,$strand,$rnao,$rnae,$cdso,$cdse,$exnum,$exono,$exone) = split /\s+/o, $_;
		


        ##################################################################################
        #################CREACIO I DEFINICIO DEL HASH %REFGENE############################
        ##################################################################################
		
		unless (defined($refGene{$chr})) { %{$refGene{$chr}} = () };
		
        #en el cas de que el hash no estigui definit li assignem una paraula clau que sera la segona 
        #posicio 1 del vector definit, es a dir, el cromosoma
		
		unless (defined($refGene{$chr}{$gid})) { %{$refGene{$chr}{$gid}} = () };

        #el segon hash si no te definida una paraula clau li assignem la posicio 0 del vector, 
        #es a dir, el ID del gen
		
		
		$refGene{$chr}{$gid}{CDSstart}   = $cdso + 1;
		
		$refGene{$chr}{$gid}{CDSend}     = $cdse;
		
		$refGene{$chr}{$gid}{TRANSstart} = $rnao + 1;
		
		$refGene{$chr}{$gid}{TRANSend}   = $rnae;
		
		$refGene{$chr}{$gid}{EXONstart}  = $exono; 
		
		$refGene{$chr}{$gid}{EXONend}    = $exone;

		$refGene{$chr}{$gid}{NUMexons}   = $exnum;
		
	}; # while

	($debug) && print STDERR Data::Dumper->Dump( [ \%refGene ], [qw( *refGene )] );
} # obrir_fitxer_refseq

#############################################################################################
##############################OBRIM EL SEGON  FITXER#########################################
#############################################################################################

sub obrir_fitxer_SNPs() {

    my $file = $_[0];

    open (FITXER2, "< $file") || do {
	print STDERR "### ERROR!!! NO PUC OBRIR FITXER $file \n";
        exit(1);
    };
	
     print STDERR "##\n## LLEGINT FITXER SNPs: $file\n##\n";

    while (<FITXER2>) {
	
	my ($chr,$posicioS,$posicioE,$id);

	next if /^\s*$/o;
	next if /^\#/o;

	chomp($_);
	
	($chr,$posicioS,$posicioE,$id) = split /\s+/o, $_;
	



	##################################################################################
        #################CREACIO I DEFINICIO DEL HASH SNPNIH  ############################
        ##################################################################################
		
		unless (defined($snpNih{$chr})) { @{$snpNih{$chr}} = () };

		push @{$snpNih{$chr}}, [$id,$posicioE];



        #en el cas de que el hash no estigui definit li assignem una paraula clau que sera la primera 
        #posicio del vector definit, es a dir, el cromosoma aquest hash apunta a un array (vector)




	    };#while

	($debug) && print STDERR Data::Dumper->Dump( [ \%snpNih ], [qw( *snpNih )] );
} # obrir_fitxer_SNPs


#############################################################################################
##############################    ORDENAR ARRAYS        #####################################
#############################################################################################



sub ordenar_array_SNPs() {

    my $chr;

    print STDERR "##\n## ORDENAR SNPs PER COORDENADES...\n##\n";

    foreach $chr (keys %snpNih){

	@{$snpNih{$chr}} = map {$_->[1]}
                           sort {$a->[0] <=> $b->[0]}
                           map { [ $_->[1], $_ ] } @{$snpNih{$chr}};

    }; # tancar el foreach

    ($debug) && print STDERR Data::Dumper->Dump( [ \%snpNih ], [qw( *snpNih )] );
} # ordenar_array_SNPs


##############################################################################################
###########################TRACTAMENT DE LES DADES############################################
##############################################################################################

sub tractament_dades() {

    my $chr;
    my $id_gene;
    my $snps;

    print STDERR "##\n## ANALITZANT DADES...\n##\n";

    foreach $chr (keys (%refGene)) {

#en aquest pas assignem a la variable $chr totes les paraules claus del hash  %refGene, que serien els
#diferents cromosomes: chr1, chr2...

	unless (defined($RESULTAT{$chr})) {%{$RESULTAT{$chr}}= ()};
	    

	foreach $id_gene (keys %{$refGene{$chr}}) {
            
	    my ($CDSs, $CDSe, $TRANSs, $TRANSe, @es, @ee, $numex);
	    my ($p, $exons, $introns, $UTR5, $UTR3, $UTRtotal, $donore, $donori, $donor, $acceptore, $acceptori, $acceptor, $Splicing, $CDS, $TOTAL);

	    $CDSs = $refGene{$chr}{$id_gene}{CDSstart};    
            $CDSe = $refGene{$chr}{$id_gene}{CDSend};    
            $TRANSs = $refGene{$chr}{$id_gene}{TRANSstart};    
            $TRANSe = $refGene{$chr}{$id_gene}{TRANSend};    
            @es = split(/,/o,$refGene{$chr}{$id_gene}{EXONstart});
            @ee = split(/,/o,$refGene{$chr}{$id_gene}{EXONend});
            $numex = $refGene{$chr}{$id_gene}{NUMexons};
	    
	    $exons = 0;
	    $introns = 0;
	    $CDS = 0;
	    $UTR5 = 0;
	    $UTR3 = 0;
	    $UTRtotal = 0;
	    $donore = 0;
	    $donori = 0;
	    $donor = 0;
	    $acceptore = 0;
	    $acceptori = 0;
	    $acceptor = 0;
	    $Splicing = 0;
	    $TOTAL = 0;
	    
	    foreach $snps (@{$snpNih{$chr}}) {
		
		my ($snpid,$pos);

		$pos = $snps->[1];

		next if ($pos < $TRANSs);
                last if ($pos > $TRANSe);
		
		
		if (($CDSs <= $pos) && ($pos <= $CDSe)) {
		    
	    
		    $p = 0; # recorre el vector @es

		    while ($p < $numex) {

			if (($es[$p] <= $pos) && ($pos <= $ee[$p])) {
			    			
			    $exons = $exons + 1;
			    
			    if (($ee[$p]-2 == $pos) || ($ee[$p]-1 == $pos) || ($ee[$p] == $pos)) {

				$donore = $donore + 1;
			    };

			    if (($es[$p] == $pos) || ($es[$p]+1 == $pos) || ($es[$p]+2 == $pos)) {
			       
				$acceptore = $acceptore + 1;
				
			    }
	    
			};
			
			
			if (($p + 1 < $numex) && ($ee[$p] < $pos) && ($pos < $es[$p + 1])) {
			    
  			    $introns = $introns + 1;
			    
 			    if (($ee[$p]+1 == $pos) || ($ee[$p]+2 == $pos) || ($ee[$p]+3 == $pos) || ($ee[$p]+4 == $pos) || ($ee[$p]+5 == $pos) || ($ee[$p]+6 == $pos)) {
				
 				$donori = $donori + 1;
			    };
			   
			    
			    
			    if (($es[$p + 1]-1 == $pos) || ($es[$p + 1]-2 == $pos) || ($es[$p + 1]-3 == $pos) || ($es[$p + 1]-4 == $pos) || ($es[$p + 1]-5 == $pos) || ($es[$p + 1]-6 == $pos) || ($es[$p + 1]-7 == $pos) || ($es[$p + 1]-8 == $pos) || ($es[$p + 1]-9 == $pos) || ($es[$p + 1]-10 == $pos) || ($es[$p + 1]-11 == $pos) || ($es[$p + 1]-12 == $pos) || ($es[$p + 1]-13 == $pos) || ($es[$p + 1]-14 == $pos) || ($es[$p + 1]-15 == $pos) || ($es[$p + 1]-16 == $pos)) {
				
				$acceptori = $acceptori + 1;
			    };
			};
			
			
			$donor = $donore + $donori;

			$acceptor = $acceptore + $acceptori;
			
			
			$p = $p + 1;
			
		    } # while ($p < $numex)
		    
		    $Splicing = $donor + $acceptor;
		    
		    
		    $CDS = $exons + $introns;
		    

		} else { # ! ($CDSs <= $pos) && ($pos <= $CDSe)

		    $UTRtotal = $UTRtotal + 1;
		    
		    if ($pos < $CDSs) {
			
			$UTR5 = $UTR5 + 1;
			
		    } else {
			
			$UTR3 = $UTR3 + 1;
			
		    };
		    		    
		};
		
		$TOTAL = $UTRtotal + $CDS;
		
	    } #foreach SNPs
	    
	    
	    unless (defined($RESULTAT{$chr}{$id_gene})) { @{$RESULTAT{$chr}{$id_gene}} = () };
	    
	    push @{$RESULTAT{$chr}{$id_gene}}, [$UTR5,$UTR3,$UTRtotal,$exons,$introns,$CDS,$donor,$acceptor,$Splicing,$TOTAL];

	
	    
	}; #foreach IDgene

    }; #foreach chr

    ($debug) && print STDERR Data::Dumper->Dump( [ \%RESULTAT ], [qw( *RESULTAT )] );
} # tractament_dades



##############################################################################################
###########################  IMPRIMINT RESULTATS  ############################################
##############################################################################################


sub imprimir_resultats() {
    
    my $chr;
    my $id_gene;
    my ($CDSs, $CDSe, $TRANSs, $TRANSe, @es, @ee, $numex);
    my $valors;
    my ($UTR5,$UTR3,$UTRtotal,$EXONs,$INTRONs,$CDS,$DONORs,$ACCEPTORs,$Splicing,$TOTAL);
  
    
    print STDERR "##\n## IMPRIMINT DADES...\n##\n";
    
    foreach $chr ((keys %RESULTAT ) && (keys %refGene )) {
	
	foreach $id_gene ((keys %{$RESULTAT{$chr}}) && (keys % {$refGene{$chr}}))  {
	    
	    $CDSs = $refGene{$chr}{$id_gene}{CDSstart};    
            $CDSe = $refGene{$chr}{$id_gene}{CDSend};    
            $TRANSs = $refGene{$chr}{$id_gene}{TRANSstart};    
            $TRANSe = $refGene{$chr}{$id_gene}{TRANSend};    
            @es = split(/,/o,$refGene{$chr}{$id_gene}{EXONstart});
            @ee = split(/,/o,$refGene{$chr}{$id_gene}{EXONend});
            $numex = $refGene{$chr}{$id_gene}{NUMexons};
	    
	    
	    
	    foreach $valors (@{$RESULTAT{$chr}{$id_gene}}) {
		
		$UTR5 = $valors -> [0];

		$UTR3 = $valors -> [1];

		$UTRtotal = $valors -> [2];

		$EXONs = $valors -> [3];
		
		$INTRONs = $valors -> [4];

		$CDS = $valors -> [5];

		$DONORs = $valors -> [6];
		
		$ACCEPTORs = $valors -> [7];

		$Splicing = $valors -> [8];

		$TOTAL = $valors -> [9];
		
		
		print "> $id_gene ($chr: $TRANSs\-$TRANSe)";
		print "\t $UTR5";
		print "\t $UTR3";
		print "\t $UTRtotal";
		print "\t $EXONs";
		print "\t $INTRONs";
		print "\t $CDS";
		print "\t $DONORs";
		print "\t $ACCEPTORs";
		print "\t $Splicing";
		print "\t $TOTAL\n";
		
		
	    } # foreach $valors
	    
	} # foreach $id_gene
	
    } # foreach $chr 
    
    ($debug) && print STDERR Data::Dumper->Dump( [ \%RESULTAT ], [qw( *RESULTAT )] );
    
} #tancar funcio imprimir dades