|
0
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 if(@ARGV<2) { print "usage: $0 gNIR_file EUMA_file > output.gFVKM\n"; exit; }
|
|
|
4 my ($gNIR_file,$EUMA_file) = @ARGV;
|
|
|
5
|
|
|
6 open IN, $gNIR_file or die "Can't open gNIR_file $gNIR_file\n";
|
|
|
7 <IN>; # skip header
|
|
|
8 while(<IN>){
|
|
|
9 chomp;
|
|
|
10 split/\t/;
|
|
|
11 $RNAseq{$_[0]}=$_[1];
|
|
|
12 }
|
|
|
13 close IN;
|
|
|
14
|
|
|
15 open IN, $EUMA_file or die "Can't open EUMA file $EUMA_file\n";
|
|
|
16 while(<IN>){
|
|
|
17 chomp;
|
|
|
18 split/\t/;
|
|
|
19 $EUMA{$_[0]}=$_[1];
|
|
|
20 }
|
|
|
21 close IN;
|
|
|
22
|
|
|
23
|
|
|
24 print "gene\tgFVK\tread.count\tgEUMA(bp)\n";
|
|
|
25 for my $gene (sort keys %EUMA){
|
|
|
26 next if $EUMA{$gene} ==0;
|
|
|
27 if(!exists $RNAseq{$gene}){ $RNAseq{$gene} = 0; }
|
|
|
28 $gFVKM = $RNAseq{$gene} / $EUMA{$gene} * 1000;
|
|
|
29 print "$gene\t$gFVKM\t$RNAseq{$gene}\t$EUMA{$gene}\n";
|
|
|
30 }
|
|
|
31
|
|
|
32
|
