|
0
|
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;
|
