-
Notifications
You must be signed in to change notification settings - Fork 6
/
TPM.pl
47 lines (44 loc) · 1 KB
/
TPM.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
47
#!/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,
'featurecount_summary|B=s' => \$featurecount_summary,
'outfile|G=s' => \$outfile,
);
open A,"< $featurecount" or die"$!";
open B,"< $featurecount_summary" or die"$!";
open OUT,"> $outfile" or die"$!";
print OUT "geneName\tTPM\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]/$a[5];
}
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])*1000000)/$to;
print OUT "\n";
}
close OUT;