#!/usr/bin/perl
# Compute gEUMA for each gene from the gUtable file and insertlendis file and print out only genes whose gEUMA is larger than 0.

if(@ARGV<3) { print "usage: $0 gUtable_file insertlendis_file READLENGTH\n"; exit; }

my ($gUtable_file,$insertlendis_file,$READLENGTH) = @ARGV;


my $totalfreq=0;
open LENDIS, $insertlendis_file or die "Can't open insertlendis file $insertlendis_file\n";
while(<LENDIS>){
 chomp;
 my ($insert_size, $freq) = split/\t/;
 $insertlendis{"$insert_size"}=$freq;
 $totalfreq+=$freq;
}
close LENDIS;



print "gene\tgEUMA\n";  #header
open GUTABLE, $gUtable_file or die "Can't open gUtable file $gUtable_file\n";
chomp($header = <GUTABLE>); 
@distances = split/\t/,$header;
shift @distances;   #'gene'
while(<GUTABLE>){
 chomp;
 my @line = split/\t/;
 my $gene = shift @line;
 my $EUMA=0;
 for my $i (0..$#distances){
    $d=$distances[$i];
    $pUMA = $line[$i] * ($insertlendis{$d}/$totalfreq) ;
    $EUMA += $pUMA;
 }
 if($EUMA>=0.005) { printf "$gene\t%.2f\n",$EUMA; }  # use >=0.005  instead of >0, because I'm plotting two decimals after point. (in bp)
}
close GUTABLE;