#!/usr/bin/perl -w
#
# localitzador.pl
#
#   filtering snps overlapping genomic features
#
# bmorancho & jplanas - UPF - 2002/02/28
# 

use strict;

# use Data::Dumper;

#
# Setting VARS

my (@SNPS,%snps_in_gene,%snps_distance,$snp_file,$seq_file,$starting_chr,@last_line,$last_genid,$last_ex_ini,$storer,$last_ex_end,$chrnumber);
my (%global_count,%chr_count,%global_len,%chr_len,%global_dist,%chr_dist,%hot_points,%hot_count,%gen_len,%snps_to_pop,%last_len);
my $chrcounter;
my $PROG = "$0";
@SNPS = ();

#
# MAIN

&parse_argvs();
&read_snps_chr_by_chr();
&global();

exit(0);

#
# FUNCTIONS

sub parse_argvs() {
    if (scalar(@ARGV < 2)) {
        print "USAGE: $PROG <taulasnps> <taulaseqs>\n";
        exit(1);
    }
    ($snp_file,$seq_file) = @ARGV;    
    # If last parameter is missing it si 1 by deffect
    if (scalar(@ARGV == 2)) {
	($snp_file,$seq_file) = @ARGV;
	$starting_chr = 1;
    };
    if (scalar(@ARGV > 2)) {
	($snp_file,$seq_file,$starting_chr,undef) = @ARGV;
    };
    $chrcounter = $starting_chr;
    &starter;
    &sorter ($snp_file, 1);
    &sorter ($seq_file, 0);
} # parse_argvs

sub sorter() {
    my ($file, $column) = @_;

    my (%saver, @f, @k);

    open(MYFILE, "< $file") || do {
	print STDERR "$PROG: no es pot carregar $file\n";
	exit(1);
    };
    
    while (<MYFILE>) {
	next if /^\d+\s*\t*$/o;
	next if /random\t/o;
        next if /^\s*$/o; # skip empty lines
	@f = split /\t/o, $_;
        $f[$column] = substr($f[$column],3);	
	if ($f[$column] eq "X") {$f[$column] = 23;};
	if ($f[$column] eq "Y") {$f[$column] = 24;};
	
	push @{ $saver{$f[$column]} }, $_;
	&counter($.)
   };
    close (MYFILE);
			     
    @k = sort { $a <=> $b }(keys(%saver));
    $chrnumber = scalar(@k);

    open(OUTFILE, "> $file");

    my $counter = 0;

    while ($counter < scalar(@k)) {
	print OUTFILE @{ $saver{$k[$counter]} }, "\n";
	$counter = $counter + 1;
    };
    close (OUTFILE);
    return;
} #sorter

sub starter() {
    my ($start,$endsub);
    
    print STDERR"\n*****************************************************************\n";
    print STDERR"************************ SNPs LOCATOR ***************************\n";
    print STDERR"*****************************************************************\n\n";
    print STDERR"This program will take a wile to be finished. Wait, please\n";
    
    $endsub = 0;
    while (!$endsub) {
	print STDERR"If you want to exit press 'q'\nWhen you would be ready to start, press ENTER\n";
	$start = <STDIN>;
	$start =~ s/(\r)?\n//o;

	if ($start eq "") {
	    $endsub = 1;
	};
	if ($start eq "q") {
	    $endsub = 1;
	    exit (0); 
	};
    }; # while
} # starter

sub read_snps_chr_by_chr() {
    my $check = 1;   
    open(SNP_FILE, "< $snp_file") || do {
        print STDERR "$PROG: impossible to load $snp_file\n";
        exit(1);
    };

    my $previous_chr = $starting_chr - 1;

    while (<SNP_FILE>) {
        my (@f,$snpstart,$snpend,$snpid,$current_chr);
        next if /^\s*$/o; # skip empty lines
        # next if /^#/o; # only if commented lines in source
        s/(\r)?\n//o; # eliminate "\n"
	s/(^M)?\n//o;
        @f = split /\t/o, $_; # suposing taht each line has only 5 fields...
        $f[1] = substr($f[1],3);  # Eliminates the three leading characters of the column "chrXX"

        $current_chr = $f[1];
        $snpstart = $f[2];
        $snpend = $f[3];
        $snpid = $f[4];

	next if $f[1] !~m/^(\d)?\d$/o;

	#CHECKING
	# If it is the first time chacks that $starting_chr (from @ARGV) is the same than
	# the first crhomosome in SNPs table
	if ($current_chr > $starting_chr + 1 && $check) {
	    die "imput files don't match with starting chromosome\n",
                "starting chromosome arg lower than starting chromosome in files\n\n";
	};
	$check = 0; # CHECKING is not needed any more.
	#If the present chromosome is the same than in the last loop cycle, ACUMULATE DATA IN @SNPS
        if  ($current_chr == $previous_chr) {
            push @SNPS, ($snpstart,$snpend,$snpid);
        }

	# If is not, chromosome has changed:
	  # if that is the first time: ONLY WRITE CHROMOSOME NAME
	  # if other times: 
	    # enter in  read_seqs amb with the last cycle loop crhomosome as a parameter 
	    # write results for that chromosome
	    # escriu el nom del cromosoma seguentwrite next chromosome name
	    # initialize variables for new chromosome
            # acumulate data from this cycle in @SNPS
	    # changes anterior chromosome to current chromosome
        else {
	    if ($previous_chr == $starting_chr - 1) {
		print "\n";
		&separador();
		print "CHROMOSOME $current_chr:\n";
		$previous_chr = $current_chr;
		next;};
            &read_seqs ($previous_chr);
	    &print_results ();
	    print "\n";
	    &separador();
            print "CHROMOSOME $current_chr:\n";
            @SNPS = ();
	    %snps_in_gene = ();
	    %gen_len = ();
	    %hot_points = ();
	    %hot_count = ();
            push @SNPS, ($snpstart,$snpend,$snpid);
            $previous_chr = $current_chr;
        }; 
    }; #while SNP_FILE

    &read_seqs($previous_chr); # to enter the last time
    close(SNP_FILE);
    &read_seqs (0); # to close SEQS
    &print_results (); # of the last chromosome
} # read_snps

sub read_seqs() {
    my ($snps_chr) = @_;

    if ( ! $snps_chr) {
        close(SEQS);
        return;
    };

    if ($snps_chr eq $starting_chr) {
       open(SEQS, "< $seq_file") || do {
       print STDERR "$PROG: impossible to load $seq_file\n";
       exit (1);
        };
   }; # if ! $snps_chr

    my (@f,$genid,$ex_ini,$ex_end,$current_chr,%lastintron,$previous_genid,$previous_ex_ini,$previous_ex_end,$ant_genid,$ant_ex_ini,$ant_ex_end,$flag);

    foreach my $sec (qw/ intron exon pre post /) {
	if ($sec eq 'pre' || $sec eq 'intron') {
	    { $chr_count{$sec} = 1 };
	}
	else {
	    { $chr_count{$sec} = 0 };
	}; # if/else
    }; # foreach
    
    # si exists a single SEQS line from the previous chromosome because the chromosome has changed, process it
    if ($last_line[1] ) {

	@f = split /\t/o, $last_line[0]; # suposing that each line has only  9 fields...
	$f[0] = substr ($f[0] , 3);
        $current_chr = $f[0];
        $ex_ini = $f[3];
        $ex_end = $f[4];
        $genid = $f[8];

	next if $f[0] !~m/^(\d)?\d$/o;
	
	print STDERR "$current_chr\t$. : $genid $ex_ini $ex_end : ";
	# Starting new chromosome:
	    # first exon of this crhomosome  without acceptor box
	    print STDERR ":";
	    &checkcoords($genid,$ex_ini,($ex_end - 4),'exon',0);
	    # donor box of first exon of this crhomosome	    
	    print STDERR "::";
	    &checkcoords($genid,($ex_end - 3),($ex_end + 3),'post',0);
	    print STDERR "\n";

	    $previous_genid = $genid;
	    $previous_ex_ini = $ex_ini;
	    $previous_ex_end = $ex_end;

	    $lastintron{$genid} = $ex_end + 4;
	
	$last_line[1] = 0;
    }; #if LAST

    while (<SEQS>) {
	next if /^\s*$/o; # skip empty lines
        # next if /^#/o; # only f commented lines in source
        s/(\r)?\n//o; # treiem "\n"
	s/(^M)?\n//o;
        @f = split /\t/o, $_; # suposing that each line has only 9 fields...
	$f[0] = substr ($f[0] , 3);
        $current_chr = $f[0];
        $ex_ini = $f[3];
        $ex_end = $f[4];
        $genid = $f[8];
        print STDERR "$current_chr\t$. : $genid $ex_ini $ex_end : ";	

	# Recovering the first line of current chromosome from last line of previous &read_seqs
        if ($current_chr > $snps_chr ) {
	    @last_line = ($_,1);
	    $last_genid = $genid;
	    $last_ex_ini = $ex_ini;
	    $last_ex_end = $ex_end;
	    last;
	};
	
        if ($current_chr < $snps_chr && $current_chr ne $starting_chr) {
            print STDERR "fatal error in snps file\n";
            exit (1);
        };

        if (! defined($snps_in_gene{$genid})) { 
            %{ $snps_in_gene{$genid} } = ();
        };

	# deletes bad data type of in the last exon of the previous gene & sets first exon and donor of the new gene:
	if (defined($previous_genid) && $previous_genid ne $genid) {
	    my $down_counter;
	    foreach my $sec (qw/ intron exon pre post /) {
		$hot_count{$sec} = ();
		if (! defined(@{ $snps_distance{$previous_genid}{$sec} })) { @{ $snps_distance{$previous_genid}{$sec} }  = (); }; 
		&snps_distance_mean($previous_genid, $sec);
	    }; 

	    # special case: singletons:
	    if ($flag) {
		$gen_len{$previous_genid}{'post'} = ($gen_len{$previous_genid}{'post'} - $last_len{'post'});
		$chr_len{'post'} = ($chr_len{'post'} - $last_len{'post'});

		if ($snps_to_pop{$previous_genid}{'post'}) {
		    my $popping;
		    $down_counter = ($snps_to_pop{$previous_genid}{'post'} - 1);
		    while ($snps_to_pop{$previous_genid}{'post'}) {
			pop(@{ $snps_in_gene{$previous_genid}{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			$snps_to_pop{$previous_genid}{'post'} = $snps_to_pop{$previous_genid}{'post'} - 1;
		    };
		    while ($down_counter) {
			$popping = 1;
			pop(@{ $snps_distance{$previous_genid}{'post'} });
			$down_counter = $down_counter - 1;
		    };
		    if ($popping) { pop(@{ $snps_distance{$previous_genid}{'post'} }); };
		}; # if
		$gen_len{$previous_genid}{'exon'} = ($gen_len{$previous_genid}{'exon'} - $last_len {'exon'});
		$chr_len{'exon'} = ($chr_len{'exon'} - $last_len {'exon'});
		if ($snps_to_pop{$previous_genid}{'exon'}) {
		    my $popping;
		    $down_counter = ($snps_to_pop{$previous_genid}{'exon'} - 1);
		    while ($snps_to_pop{$previous_genid}{'exon'}) {
			pop(@{ $snps_in_gene{$previous_genid}{'exon'} });
			$snps_to_pop{$previous_genid}{'exon'} = $snps_to_pop{$previous_genid}{'exon'} - 1;
		    };
		    while ($down_counter) {
			$popping = 1;
			pop(@{ $snps_distance{$previous_genid}{'exon'} });
			$down_counter = $down_counter - 1;
		    };
		    if ($popping) { pop(@{ $snps_distance{$previous_genid}{'exon'} }); };
		}; # if
		&checkcoords($previous_genid,$previous_ex_ini,$previous_ex_end,'exon',1);
		&snps_distance_mean($previous_genid, 'exon');
	    } #if flag
	    else {
	    
		# deletes record of donor box in last exon of previous gene
		$gen_len{$previous_genid}{'post'} = ($gen_len{$previous_genid}{'post'} - $last_len{'post'});
		$chr_len{'post'} = ($chr_len{'post'} - $last_len{'post'});

		if ($snps_to_pop{$previous_genid}{'post'}) {
		    my $popping;
		    $down_counter = ($snps_to_pop{$previous_genid}{'post'} - 1);
		    while ($snps_to_pop{$previous_genid}{'post'}) {
			pop(@{ $snps_in_gene{$previous_genid}{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			pop(@{ $hot_points{'post'} });
			$snps_to_pop{$previous_genid}{'post'} = $snps_to_pop{$previous_genid}{'post'} - 1;
		    };
		    while ($down_counter) {
			$popping = 1;
			pop(@{ $snps_distance{$previous_genid}{'post'} });
			$down_counter = $down_counter - 1;
		    };
		    if ($popping) { pop(@{ $snps_distance{$previous_genid}{'post'} }); };
		}; # if

		# deletes record of last exon of previous gene (it will be replaced 3 lines below)
		$gen_len{$previous_genid}{'exon'} = ($gen_len{$previous_genid}{'exon'} - $last_len {'exon'});
		$chr_len{'exon'} = ($chr_len{'exon'} - $last_len {'exon'});
		if ($snps_to_pop{$previous_genid}{'exon'}) {
		    my $popping;
		    $down_counter = ($snps_to_pop{$previous_genid}{'exon'} - 1);
		    while ($snps_to_pop{$previous_genid}{'exon'}) {
			pop(@{ $snps_in_gene{$previous_genid}{'exon'} });
			$snps_to_pop{$previous_genid}{'exon'} = $snps_to_pop{$previous_genid}{'exon'} - 1;
		    };
		    while ($down_counter) {
			$popping = 1;
			pop(@{ $snps_distance{$previous_genid}{'exon'} });
			$down_counter = $down_counter - 1;
		    };
		    if ($popping) { pop(@{ $snps_distance{$previous_genid}{'exon'} }); };
		} # if
		# replacing the last exon of previous gene without acceptor box
		&checkcoords($previous_genid,($previous_ex_ini + 4),$previous_ex_end,'exon',1);
		# first exon of gene without acceptor box
		print STDERR ":";
		&checkcoords($genid,$ex_ini,($ex_end - 4),'exon',0);
		# donor box of first exon of gene
		print STDERR "::";
		&checkcoords($genid,($ex_end - 3),($ex_end + 3),'post',0);
		print STDERR "\n";
		
		$ant_genid = $previous_genid;
		$ant_ex_ini = $previous_ex_ini;
		$ant_ex_end = $ant_ex_end;
		
		$previous_genid = $genid;
		$previous_ex_ini = $ex_ini;
		$previous_ex_end = $ex_end;
		
		$lastintron{$genid} = $ex_end + 4;
		
		$flag = 1;
		next;
	    }; #if/else
	}; #if (gene change)
	
	$flag = 0;
	
	if(! defined($previous_genid)) {
	    print STDERR ":";
	    # first exon of gene without donor box
	    &checkcoords($genid,$ex_ini,($ex_end - 4),'exon',0);
	    print STDERR "::";
	    # donor box of first exon of gene	    
	    &checkcoords($genid,($ex_end - 3),($ex_end + 3),'post',0);
	    print STDERR "\n";
	    
	    $previous_genid = $genid;
	    $previous_ex_ini = $ex_ini;
	    $previous_ex_end = $ex_end;
	    
	    $lastintron{$genid} = $ex_end + 4;
	    
	    next;
	};
	
	# General cases
	
	if (defined($lastintron{$genid}) && $lastintron{$genid} < $ex_ini-4) {
	    &checkcoords($genid,$lastintron{$genid},($ex_ini - 4),'intron',0);
	}
	
	$lastintron{$genid} = $ex_end + 4;
	
	print STDERR ":";
	&checkcoords($genid,$ex_ini+4,$ex_end-4,'exon',0);
	print STDERR ":";
	&checkcoords($genid,$ex_ini-3,$ex_ini+3,'pre',0);
	print STDERR ":";
	&checkcoords($genid,$ex_end-3,$ex_end+3,'post',0);
	print STDERR "\n";	    
	
	if ($.>1) {
	    $ant_genid = $previous_genid;
	    $ant_ex_ini = $previous_ex_ini;
	    $ant_ex_end = $previous_ex_end;
	}
	
	$previous_genid = $genid;
	$previous_ex_ini = $ex_ini;
	$previous_ex_end = $ex_end;
	
	$storer = $_;
	next;
    
    }; # while SEQS

    @f = split /\t/o, $storer; # suposing that each line has only 9 fields...
    $f[0] = substr ($f[0] , 3);
    $current_chr = $f[0];
    $ex_ini = $f[3];
    $ex_end = $f[4];
    $genid = $f[8];
    
    # deletes bad data type of in the last exon of a crhomosome:
    # deletes record of donor box in last exon of a crhomosome
    $gen_len{$genid}{'post'} = ($gen_len{$genid}{'post'} - $last_len {'post'});
    $chr_len{'post'} = ($chr_len{'post'} - $last_len {'post'});
    my $down_counter;
    if ($snps_to_pop{$genid}{'post'}) {
	my $popping;
	$down_counter = ($snps_to_pop{$genid}{'post'} - 1);
	while ($snps_to_pop{$genid}{'post'}) {
	    pop(@{ $snps_in_gene{$genid}{'post'} });
	    pop(@{ $hot_points{'post'} });
	    pop(@{ $hot_points{'post'} });
	    pop(@{ $hot_points{'post'} });
	    pop(@{ $hot_points{'post'} });
	    $snps_to_pop{$genid}{'post'} = $snps_to_pop{$genid}{'post'} - 1;
	};
	while ($down_counter) {
	    $popping = 1;
	    pop(@{ $snps_distance{$genid}{'post'} });
	    $down_counter = $down_counter - 1;
	};
	if ($popping) { pop(@{ $snps_distance{$previous_genid}{'post'} }); };
    }; # if
      # deletes record of last exon of this crhomosome (it will be replaced some lines below)
	$gen_len{$genid}{'exon'} = ($gen_len{$genid}{'exon'} - $last_len {'exon'});
	$chr_len{'exon'} = ($chr_len{'exon'} - $last_len {'exon'});
   
    if ($snps_to_pop{$genid}{'exon'}) {
	my $popping;
	$down_counter = ($snps_to_pop{$genid}{'exon'} - 1);
	while ($snps_to_pop{$genid}{'exon'}) {
	    pop(@{ $snps_in_gene{$genid}{'exon'} });
	    $snps_to_pop{$genid}{'exon'} = $snps_to_pop{$genid}{'exon'} - 1;
	};
	while ($down_counter) {
	    $popping = 1;
	    pop(@{ $snps_distance{$genid}{'exon'} });
	    $down_counter = $down_counter - 1;
	};
	if ($popping) { pop(@{ $snps_distance{$previous_genid}{'exon'} }); };
    }; # if
      &checkcoords($genid,($ex_ini +4),$ex_end,'exon',1,1);
      &snps_distance_mean($genid, 'exon');

    return;
} # read_seqs

sub checkcoords() {
    my ($genid, $ini, $end, $sec, $backwards, $last) = @_;
    my $i = 0;
    my $current_snp = 0; 
    my $previous_snp = 0;
    
    $snps_to_pop{$genid}{$sec} = ();
    
    if (! defined($chr_count{$sec})) {
    }; # if
    $chr_count{$sec} ++;
    if (! defined($hot_count{$sec})) { 
	if ($sec eq 'pre' || $sec eq 'intron') {
	    { $hot_count{$sec} = 1 };
	}
	else {
	    { $hot_count{$sec} = 0 };
	}; # if/else
    };
    $hot_count{$sec} ++;
    if ($ini >= $end) { return };
    if (!$backwards) {$backwards = 0};
    if (!$last) {$last = 0};
    
    if (! defined($chr_len{$sec})) { $chr_len{$sec} = 0 };
    $chr_len{$sec} = ($chr_len{$sec} + ($end - $ini) ); 
    $last_len{$sec} = ($end - $ini);

    if (! defined($gen_len{$genid}{$sec})) { $gen_len{$genid}{$sec} = 0 };
    $gen_len{$genid}{$sec} = ($gen_len{$genid}{$sec} + ($end - $ini) ); 

    while ($i < scalar(@SNPS)) {
	if ($backwards && $ini > $SNPS[$i]) {
	    $i = $i + 3;
	    next; }; 
	if (!$backwards && $ini > $SNPS[$i+1]) {
	    $i = $i + 3;
	    next; }; 
	if (! $backwards && $end < $SNPS[$i]) { last; };
	if ($backwards && $end < $SNPS[$i]) { last; };
	if (! defined($snps_in_gene{$genid}{$sec})) { 
	    @{ $snps_in_gene{$genid}{$sec} } = ();
	}; 
	if (!$backwards) { 
	    print STDERR "X";
	};

	if (!defined $snps_to_pop{$genid}{$sec}) { $snps_to_pop{$genid}{$sec} = () };
	$snps_to_pop{$genid}{$sec} ++;

	if ($backwards)  { print STDERR "\n"; };
	push @{ $snps_in_gene{$genid}{$sec} }, $SNPS[$i+2];
	
	$current_snp = $SNPS[$i];
	if ($previous_snp) {
	    push @{ $snps_distance{$genid}{$sec} },($current_snp - $previous_snp);
	};
	
	if ($sec eq 'post') { 
	    push @{ $hot_points{$sec} },($genid,$hot_count{$sec},$sec,$SNPS[$i+2]);
	};
	
	if ($sec eq 'pre') {
	    push @{ $hot_points{$sec} }, ($genid,($hot_count{$sec} ),$sec,$SNPS[$i+2]);
	};
	
	
	$previous_snp = $current_snp;
        $i = $i + 3;
    }; # while

    return;
} # checkcoords

sub snps_distance() {
    my ($genid, $sec) = @_;
    my $counter = 0;
    my $acumulator = 0;

    my $ref = \@{ $snps_distance{$genid}{$sec} };
    
    while ($counter < scalar(@{ $ref })) {
	$acumulator = $acumulator + $ref->[$counter];
	$counter = $counter + 1;;
    }; 

    return;
} # snps_distance

sub snps_distance_mean() {
    my ($genid, $sec) = @_;
    my $counter = 0;
    my $acumulator = 0;

    my $ref = \@{ $snps_distance{$genid}{$sec} };

    while ($counter < scalar(@{ $ref })) {
	$acumulator = $acumulator + $ref->[$counter];
	$counter = $counter + 1;;
    };        

    if (scalar(@{ $ref }) > 1) {
        push @{ $ref }, $acumulator/(scalar(@{ $ref } - 1));
    };    

    if (! defined(@{ $chr_dist{$sec} })) { @{ $chr_dist{$sec} } = (0,0) };
    my $tref = \@{ $snps_distance{$genid}{$sec} };
    my $sref = \@{ $chr_dist{$sec} };
    if (! defined ($tref->[$#{ $tref }])) { push @{$tref }, 0};
    
    $sref->[0] = ($sref->[0] + $tref->[$#{ $tref }]);
    $sref->[1] ++;

    return;
} # snps_distance_mean

sub print_results() {
    # print Data::Dumper->Dump([ \%snps_in_gene ], qw( *snps_in_gene ));   
    my @ary = keys %snps_in_gene;
    my $n = 0;	
    my $output_file_name = "Chr".$chrcounter."_raw.out";	
    open(OUTPUT, "> $output_file_name");
    while ($n<scalar(@ary)) {
        my $gene = $ary[$n];
        if (! defined($snps_in_gene{$gene})) { next; };
	if (! defined($snps_distance{$gene})) { next; };
        &separador();
	foreach my $key (qw/ intron exon pre post /) {
	    if (! defined($snps_in_gene{$gene}{$key})) { 
		@{ $snps_in_gene{$gene}{$key} } = ();
	    };
	    my $ref = \@{ $snps_distance{$gene}{$key} };
	    my $sref = \@{ $snps_in_gene{$gene}{$key} };
	    my $tref = \@{ $chr_dist{$key} };
	    my $uref = $gen_len{$gene}{$key};
	    if (! defined ($ref->[$#{ $ref }])) { push @{$ref }, 0};
	    
	    printf STDOUT "%-15s %-6s", $gene, $key;
	    print OUTPUT "$gene\t$key\t";
	    if ($uref) {
		 printf STDOUT "%8d %5d (%.2f)\t", 
		 $uref, scalar(@{ $sref }),$ref->[$#{ $ref }];
		 print OUTPUT $uref,"\t", scalar(@{ $sref }),"\t",$ref->[$#{ $ref }],"\t";
	     }
	    else {
		 printf STDOUT "%s %5d (%.2f)\t", 
		"NO $key\t", scalar(@{ $sref }),$ref->[$#{ $ref }];
		 print OUTPUT "NO $key\\t", scalar(@{ $sref }),"\t",$ref->[$#{ $ref }],"\t";
	     }; #if/else
	    if ($uref) {
		printf STDOUT "%.2f\t%s\n",
		(scalar(@{ $sref })/($uref/1000)) ,"@{ $sref }\n"; 
		print OUTPUT (scalar(@{ $sref })/($uref/1000)) ,"\t",(defined(@{ $sref }) ? "@{ $sref }" : ""),"\n";
	    }
	    else {
		print STDOUT "%s\t%s\n",
		"NO $key\t" ,"@{ $sref }\n"; # join(" ", @{ $snps_in_gene{$gene}{$key} });
		print OUTPUT "NO $key\t",(defined(@{ $sref }) ? "@{ $sref }" : ""),"\n";
	    };# if/else
	    	    
	    if (! defined($global_count{$key})) { $global_count{$key} = 0 };
	    $global_count{$key} = ($global_count{$key} + $chr_count{$key});
	    
	    if (! defined($global_len{$key})) { $global_len{$key} = 0 };
	    $global_len{$key} = ($global_len{$key} + $chr_len{$key});
	    
	    if (! defined(@{ $global_dist{$key} })) { @{ $global_dist{$key} } = (0,0) };
	    @{ $global_dist{$key} }[0] = (@{ $global_dist{$key} }[0] + @{ $chr_dist{$key} }[0]); 
	    @{ $global_dist{$key} }[1] = (@{ $global_dist{$key} }[1] + @{ $chr_dist{$key} }[1]); 
	};# foreach
	    $n = $n + 1;
    };# while

    print OUTPUT"\n";
    foreach my $key (qw/ intron exon pre post /) {
	my $ref = \@{ $chr_dist{$key} };
	my $tref = $chr_count{$key};
	print OUTPUT $key,"\t",$tref,"\t",$chr_len{$key},"\t";
	if (!defined($tref)) {
	    print OUTPUT "NONE\t";
	}
	else { print OUTPUT $tref/($chr_len{$key}/1000),"\t" }; #if/else

	if ($ref->[0] == 0 || $ref->[1] == 0) {
	    print OUTPUT "NONE\n";
	}
	else {print OUTPUT ($ref->[0]/$ref->[1]),"\n"; }; #if/else
    }; #foreach
    close(OUTPUT);
                   
    my $hot_file_name = "Chr".$chrcounter."_hot.out";
    open(HOTFILE, "> $hot_file_name");
    print HOTFILE "CROMOSOMA $chrcounter\n";
    print HOTFILE "GEN\tEXO\tLLOC\tID\n";
    foreach my $key (qw/ pre post /) {
	my $ref = \@{ $hot_points{$key} };
	if (! defined(@{ $ref} )) { @{ $ref } = () };
	#print " TOTS ELS HOTPOINTS: @{ $ref }\n";
	my $hot_counter = 0;
	while ($hot_counter < scalar(@{ $ref })) {
	    print HOTFILE @{ $ref }[$hot_counter],"\t",@{ $ref }[$hot_counter + 1],"\t",
	    @{ $ref }[$hot_counter + 2],"\t",@{ $ref }[$hot_counter + 3],"\n";
	    $hot_counter = $hot_counter + 4;
	};
    }; # foreach
    close(HOTFILE);

    $chrcounter ++;
    &separador();
} # print_results

sub separador() {
    print STDOUT ("-" x 75)."\n";
}

sub counter($) {
    my $c = shift(@_);
    if ($c % 100 == 0) {
	print STDERR ".";
	if ($c % 5000 == 0) {print STDERR"[$c]\n";};
    };
    return;
} # counter

sub global() { 
    open(GLOBAL, "> global.out");
    foreach my $sec (qw/ intron exon pre post /) {
	
	print GLOBAL $sec,"\t",$global_count{$sec},"\t",$global_len{$sec},"\t";
	if ($global_count{$sec} == 0 || $global_len{$sec} == 0) { print GLOBAL "NONE\t"; }
	else { print GLOBAL $global_count{$sec}/($global_len{$sec}/1000),"\t"; }; # if/else
	if (@{ $global_dist{$sec} }[0] == 0 || @{ $global_dist{$sec} }[1] == 0) {  print GLOBAL "NONE\n"; }
	else { print GLOBAL (@{ $global_dist{$sec} }[0]/@{ $global_dist{$sec} }[1]),"\n"; } #if/else
    };
    close(GLOBAL);
} #global

#
# This program was checked by Bea Morancho & Joan Planas on 2002 March 15th
# IT DO WORK.
#
#END
#