--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/NEUMA-1.2.1/calculate_gEUMA.pl	Thu Aug 08 00:46:13 2013 -0400
@@ -0,0 +1,42 @@
+#!/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;
+
+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;
+ for my $i (0..$#distances){
+   $gUtable{$gene}{"$distances[$i]"}=$line[$i];
+ }
+}
+close GUTABLE;
+
+
+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
+for my $gene (sort keys %gUtable){
+  my $EUMA=0;
+  for my $d (keys %{$gUtable{$gene}}){
+    $pUMA = $gUtable{$gene}{$d} * ($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)
+}
+