Mercurial > repos > chawhwa > neuma
comparison NEUMA-1.2.1/calculate_iEUMA.pl @ 0:c44c43d185ef draft default tip
NEUMA-1.2.1 Uploaded
author | chawhwa |
---|---|
date | Thu, 08 Aug 2013 00:46:13 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:c44c43d185ef |
---|---|
1 #!/usr/bin/perl | |
2 # Compute iEUMA for each isoform from the iUtable file and insertlendis file and print out only isoforms whose iEUMA is larger than 0. | |
3 | |
4 if(@ARGV<2) { print "usage: $0 iUtable_file insertlendis_file READLENGTH\n"; exit; } | |
5 | |
6 my ($iUtable_file,$insertlendis_file,$READLENGTH) = @ARGV; | |
7 | |
8 open IUTABLE, $iUtable_file or die "Can't open iUtable file $iUtable_file\n"; | |
9 chomp($header = <IUTABLE>); | |
10 @distances = split/\t/,$header; | |
11 shift @distances; #'gene' | |
12 while(<IUTABLE>){ | |
13 chomp; | |
14 my @line = split/\t/; | |
15 my $gid = shift @line; | |
16 my $transcript = shift @line; | |
17 my $gene = "$gid\t$transcript"; | |
18 for my $i (0..$#distances){ | |
19 $iUtable{$gene}{"$distances[$i]"}=$line[$i]; | |
20 } | |
21 } | |
22 close IUTABLE; | |
23 | |
24 | |
25 my $totalfreq=0; | |
26 open LENDIS, $insertlendis_file or die "Can't open insertlendis file $insertlendis_file\n"; | |
27 while(<LENDIS>){ | |
28 chomp; | |
29 my ($insert_size, $freq) = split/\t/; | |
30 $insertlendis{$insert_size}=$freq; | |
31 $totalfreq+=$freq; | |
32 } | |
33 close LENDIS; | |
34 | |
35 print "gene\ttranscript\tiEUMA\n"; #header | |
36 for my $gene (sort keys %iUtable){ | |
37 my $EUMA=0; | |
38 for my $d (keys %{$iUtable{$gene}}){ | |
39 $pUMA = $iUtable{$gene}{$d} * ($insertlendis{$d}/$totalfreq) ; | |
40 $EUMA += $pUMA; | |
41 } | |
42 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) | |
43 } | |
44 |