Mercurial > repos > chawhwa > neuma
comparison NEUMA-1.2.1/merge_LVKM_readcount.pl @ 0:c44c43d185ef draft default tip
NEUMA-1.2.1 Uploaded
author | chawhwa |
---|---|
date | Thu, 08 Aug 2013 00:46:13 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:c44c43d185ef |
---|---|
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 |