diff NEUMA-1.2.1/calculate_iEUMA.pl @ 0:c44c43d185ef draft default tip

NEUMA-1.2.1 Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/NEUMA-1.2.1/calculate_iEUMA.pl	Thu Aug 08 00:46:13 2013 -0400
@@ -0,0 +1,44 @@
+#!/usr/bin/perl
+# Compute iEUMA for each isoform from the iUtable file and insertlendis file and print out only isoforms whose iEUMA is larger than 0.
+
+if(@ARGV<2) { print "usage: $0 iUtable_file insertlendis_file READLENGTH\n"; exit; }
+
+my ($iUtable_file,$insertlendis_file,$READLENGTH) = @ARGV;
+
+open IUTABLE, $iUtable_file or die "Can't open iUtable file $iUtable_file\n";
+chomp($header = <IUTABLE>); 
+@distances = split/\t/,$header;
+shift @distances;   #'gene'
+while(<IUTABLE>){
+ chomp;
+ my @line = split/\t/;
+ my $gid = shift @line;
+ my $transcript = shift @line;
+ my $gene = "$gid\t$transcript";
+ for my $i (0..$#distances){
+   $iUtable{$gene}{"$distances[$i]"}=$line[$i];
+ }
+}
+close IUTABLE;
+
+
+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\ttranscript\tiEUMA\n";  #header
+for my $gene (sort keys %iUtable){
+  my $EUMA=0;
+  for my $d (keys %{$iUtable{$gene}}){
+    $pUMA = $iUtable{$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)
+}
+