forked from HuangLab-Fudan/ASJA
-
Notifications
You must be signed in to change notification settings - Fork 0
/
FPKM-featurecounts.pl
46 lines (43 loc) · 997 Bytes
/
FPKM-featurecounts.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#!/usr/bin/perl
# usege: perl TPM.pl -A featurecount -B featurecount.summary -O gene_TPM.txt
use strict;
use warnings;
use Cwd;
use File::Basename;
use Getopt::Long;
use List::Util qw/sum/;
my $featurecount;
my $featurecount_summary;
my $outfile;
&GetOptions (
'featurecount|A=s' => \$featurecount,
'outfile|G=s' => \$outfile,
);
open A,"< $featurecount" or die"$!";
#open B,"< $featurecount_summary" or die"$!";
open OUT,"> $outfile" or die"$!";
print OUT "geneName\tFPKM-featurecounts\n";
# my $totalReads;
# while(<B>){
# s/\s+$//;
# if ($_ =~ /Assigned\t(\d+)/){$totalReads = $1;last;}
# }
# close B;
my @total;
while(<A>){
s/\s+$//;
next if $_ =~ /[Program][Geneid]/;
my @a = split/\t/,$_;
push @total,$a[6];
}
my $to= sum @total;
close A;
open ARP,"< $featurecount" or die"$!";
while(<ARP>){
s/\s+$//;
next if $_ =~ /[Program][Geneid]/;
my @a = split/\t/,$_;
printf OUT $a[0]."\t"."%.3f",(($a[6]/$a[5])*1000000000)/$to;
print OUT "\n";
}
close OUT;