|
0
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 if(@ARGV<6) { print "usage: $0 final.iFVKM_file mapping_stat_file gene2symbol_file factor_col sample_name iLVKM_file\n"; exit; }
|
|
|
4 my ($iFVKM_file,$mapping_stat_file,$gene2symbol_file,$factor_col,$sample_name,$iLVKM_file) = @ARGV;
|
|
|
5
|
|
|
6 open STAT, $mapping_stat_file or die "Can't open mapping stat file $mapping_stat_file\n";
|
|
|
7 <STAT>; #discard header
|
|
|
8 while(<STAT>){
|
|
|
9 chomp;
|
|
|
10 next if !/\S/;
|
|
|
11 my ($sample,$factor) = (split/\t/)[0,$factor_col];
|
|
|
12 $norm_factor1{$sample} = $factor/1000000;
|
|
|
13 if($factor==0) { die "Error: Total reads is zero. Please check datatype option. Use option 'E' if your reference doesn't contain 'NM' prefix for mRNA.\n"; }
|
|
|
14
|
|
|
15 }
|
|
|
16 close STAT;
|
|
|
17
|
|
|
18 open GENE2SYM, $gene2symbol_file or die "can't open gene2symbol file $gene2symbol_file\n";
|
|
|
19 while(<GENE2SYM>){
|
|
|
20 chomp;
|
|
|
21 split/\t/;
|
|
|
22 $gene2symbol{$_[0]}=$_[1];
|
|
|
23 }
|
|
|
24 close GENE2SYM;
|
|
|
25
|
|
|
26
|
|
|
27 $"="\t";
|
|
|
28 open OUT1, ">$iLVKM_file" or die "Can't write to norm1outfile $iLVKM_file\n";
|
|
|
29 print OUT1 "gene.id\tgene.symbol\tisoform.id\tiLVKM\tiFVKM\tfinal.iFVK\torig.iFVK\tiFVKM.factor\tiEUMA\t#isoforms\t#measured_isoforms\n";
|
|
|
30
|
|
|
31 open IFVKM, $iFVKM_file or die "Can't open iFVKM file $iFVKM_file\n";
|
|
|
32 <IFVKM>; #discard header
|
|
|
33 while(<IFVKM>){
|
|
|
34 chomp;
|
|
|
35 my @line = split/\t/;
|
|
|
36 my $gene_id = shift @line;
|
|
|
37 my $isoform_id = shift @line;
|
|
|
38 my $gene_symbol = $gene2symbol{$gene_id};
|
|
|
39 $norm1=$line[0]/$norm_factor1{$sample_name};
|
|
|
40 $norm2=log($norm1+1)/log(2);
|
|
|
41 printf OUT1 "$gene_id\t$gene_symbol\t$isoform_id\t%.3f\t%.3f\t@line\n",$norm2,$norm1,;
|
|
|
42 }
|
|
|
43 close IFVKM;
|
|
|
44
|
|
|
45
|
|
|
46 close OUT1;
|
|
|
47
|
|
|
48
|
|
|
49
|
