|
0
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 if(@ARGV<4) { print "usage: $0 genewise/isoformwise[g/i] EUMAcut LVKM_out_dir NR(1/0) > output.gNIR(iNIR)\n"; exit; }
|
|
|
4 my $option = shift @ARGV; #either g or i
|
|
|
5 my $EUMAcut = shift @ARGV;
|
|
|
6 my $LVKM_dir = shift @ARGV;
|
|
|
7 my $NR_option = shift @ARGV;
|
|
|
8
|
|
|
9 my $suffix = $option."LVKM";
|
|
|
10
|
|
|
11 my @sample_list = split/\n/,`ls -1 $LVKM_dir`;
|
|
|
12 for my $i (0..$#sample_list){
|
|
|
13 if($NR_option eq '1'){
|
|
|
14 if($sample_list[$i] =~ /\.\d+.-NR.[gi]LVKM$/) { $samples{$`}=1; }
|
|
|
15 }
|
|
|
16 else {
|
|
|
17 if($sample_list[$i] =~ /\.\d+.[gi]LVKM$/) { $samples{$`}=1; }
|
|
|
18 }
|
|
|
19 }
|
|
|
20
|
|
|
21 @sample_list = keys %samples;
|
|
|
22
|
|
|
23
|
|
|
24 for my $sample (@sample_list){
|
|
|
25 next if($sample eq '');
|
|
|
26 my $file;
|
|
|
27 if($NR_option eq '1'){
|
|
|
28 $file = "$LVKM_dir/$sample.$EUMAcut.-NR.$suffix";
|
|
|
29 }
|
|
|
30 else {
|
|
|
31 $file = "$LVKM_dir/$sample.$EUMAcut.$suffix";
|
|
|
32 }
|
|
|
33
|
|
|
34 print STDERR "$file\n"; #DEBUGGING
|
|
|
35 open IN,$file or die "Can't open LVKM file $file\n";
|
|
|
36 <IN>;
|
|
|
37 while(<IN>){
|
|
|
38 chomp;
|
|
|
39
|
|
|
40 my $gene;
|
|
|
41 if($option eq 'g'){
|
|
|
42 my ($gid,$symbol,$origGFVK,$gEUMA) = (split/\t/)[0,1,5,7];
|
|
|
43 $gene = "$gid\t$symbol";
|
|
|
44 $temp = $origGFVK * $gEUMA;
|
|
|
45 if($temp == 0){
|
|
|
46 $EXPR{$gene}{$sample} = 0;
|
|
|
47 }else{
|
|
|
48 $EXPR{$gene}{$sample} = int(($temp + 0.5) /1000);
|
|
|
49 }
|
|
|
50 }
|
|
|
51 elsif($option eq 'i'){
|
|
|
52 my ($gid,$symbol,$NM,$finalIFVK,$iEUMA) = (split/\t/)[0,1,2,5,8];
|
|
|
53 $gene = "$gid\t$symbol\t$NM";
|
|
|
54 $temp = $finalIFVK * $iEUMA;
|
|
|
55 if($temp == 0){
|
|
|
56 $EXPR{$gene}{$sample} = 0;
|
|
|
57 }else{
|
|
|
58 $EXPR{$gene}{$sample} = int(($temp + 0.5) /1000);
|
|
|
59 }
|
|
|
60 }
|
|
|
61 }
|
|
|
62 close IN;
|
|
|
63 }
|
|
|
64
|
|
|
65
|
|
|
66
|
|
|
67 $"="\t";
|
|
|
68 if($option eq 'g') { print "gene.id\tgene.symbol\t@sample_list\n"; }
|
|
|
69 elsif($option eq 'i') { print "gene.id\tgene.symbol\tisoform\t@sample_list\n"; }
|
|
|
70 for my $gene (keys %EXPR){
|
|
|
71 print "$gene";
|
|
|
72 for my $sample (@sample_list){
|
|
|
73 if(exists ${$EXPR{$gene}}{$sample}){
|
|
|
74 print "\t$EXPR{$gene}{$sample}";
|
|
|
75 }
|
|
|
76 else { print "\tNA"; }
|
|
|
77 }
|
|
|
78 print "\n";
|
|
|
79 }
|
|
|
80
|
|
|
81
|
|
|
82
|
