Mercurial > repos > chawhwa > neuma

comparison NEUMA-1.2.1/calculate_gEUMA.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 gUtable file and insertlendis file and print out only genes whose gEUMA is larger than 0.
3
4 if(@ARGV<3) { print "usage: $0 gUtable_file insertlendis_file READLENGTH\n"; exit; }
5
6 my ($gUtable_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 GUTABLE, $gUtable_file or die "Can't open gUtable file $gUtable_file\n";
23 chomp($header = <GUTABLE>);
24 @distances = split/\t/,$header;
25 shift @distances; #'gene'
26 while(<GUTABLE>){
27 chomp;
28 my @line = split/\t/;
29 my $gene = shift @line;
30 my $EUMA=0;
31 for my $i (0..$#distances){
32 $d=$distances[$i];
33 $pUMA = $line[$i] * ($insertlendis{$d}/$totalfreq) ;
34 $EUMA += $pUMA;
35 }
36 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)
37 }
38 close GUTABLE;
39
40
41