Mercurial > repos > chawhwa > neuma

comparison NEUMA-1.2.1/calculate_iEUMA.2.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 gEUMA for each gene from the iUtable file and insertlendis file and print out only genes whose gEUMA is larger than 0.
3
4 if(@ARGV<3) { print "usage: $0 iUtable_file insertlendis_file READLENGTH\n"; exit; }
5
6 my ($iUtable_file,$insertlendis_file,$READLENGTH) = @ARGV;
7
8
9 my $totalfreq=0;
10 open LENDIS, $insertlendis_file or die "Can't open insertlendis file $insertlendis_file\n";
11 while(<LENDIS>){
12 chomp;
13 my ($insert_size, $freq) = split/\t/;
14 $insertlendis{"$insert_size"}=$freq;
15 $totalfreq+=$freq;
16 }
17 close LENDIS;
18
19
20
21 print "gene\tgEUMA\n"; #header
22 open IUTABLE, $iUtable_file or die "Can't open iUtable file $iUtable_file\n";
23 chomp($header = <IUTABLE>);
24 @distances = split/\t/,$header;
25 shift @distances; #'gene'
26 shift @distances; #'transcript
27 while(<IUTABLE>){
28 chomp;
29 my @line = split/\t/;
30 my $gid = shift @line;
31 my $transcript = shift @line;
32 my $gene = "$gid\t$transcript";
33
34 my $EUMA=0;
35 for my $i (0..$#distances){
36 $d=$distances[$i];
37 $pUMA = $line[$i] * ($insertlendis{$d}/$totalfreq) ;
38 $EUMA += $pUMA;
39 }
40 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)
41 }
42 close IUTABLE;