Mercurial > repos > chawhwa > neuma

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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