Mercurial > repos > bioitcore > splicetrap

comparison bin/ApplyCutoff.jie.pl @ 1:adc0f7765d85 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author bioitcore
date Thu, 07 Sep 2017 15:06:58 -0400
parents
children
comparison
equal deleted inserted replaced
0:d4ca551ca300 1:adc0f7765d85
1 #apply our cuttoff hash table on the IR calculated by Jie
2 #Modified from Martin's code
3 use strict;
4
5 use Cwd;
6 my $PROG = $0;
7 my $CUR_DIR = Cwd::abs_path(Cwd::cwd());
8 my $PROG_ABS_PATH = Cwd::abs_path($PROG);
9 my $SrcFolder=`dirname $PROG_ABS_PATH`;
10 chomp($SrcFolder);
11
12 my %cutoff;
13 my @Exlen;
14
15 my $cutoff_level=$ARGV[1];
16 my $JunctionCut = $ARGV[2];
17 my $noirm = $ARGV[3];
18
19 my $cutoff_level_index=7;
20
21 $cutoff_level_index=8 if $cutoff_level eq "H";
22 $cutoff_level_index=6 if $cutoff_level eq "L";
23
24 open(CUT,"$SrcFolder/../cutoffs/cutoff.pair.0".$cutoff_level_index.".txt") || die "cutoff file not found $!\n";
25
26 while(<CUT>){
27 chomp;
28 my @a=split(/\t/,$_);
29 push @Exlen,$a[0];
30 $cutoff{$a[0]}=$a[1];
31 }
32 close(CUT);
33
34 open(IN,$ARGV[0]);
35
36 while(<IN>){
37 chomp;
38 my @a=split(/\t/,$_);
39 my $Ez='Ez=yes';
40 my $print=$_;
41 if($a[0]=~m/#/g){next}
42 my $eventid=substr($a[0],0,2);
43 my $bir =$a[2];
44 $bir =$a[1] if($noirm eq "noirm");
45 my $j12 = $a[8];
46 my $j23 = $a[9];
47 my $j13 = $a[10];
48 my $cov1=$a[11];
49 my $cov2=$a[12];
50 my $cov3=$a[13];
51 my $siz1=$a[15];
52 my $siz2=$a[16];
53 my $siz3=$a[17];
54
55
56 my $stat1='exon1='.cutoff($siz1,$cov1,\@Exlen,%cutoff);
57 my $stat2='exon2='.cutoff($siz2,$cov2,\@Exlen,%cutoff);
58 my $stat3='exon3='.cutoff($siz3,$cov3,\@Exlen,%cutoff);
59 if($stat1 eq "exon1=yes" and $stat3 eq "exon3=yes")
60 {
61 #$Ez="passed";
62 $Ez=$eventid if $eventid eq "AA";
63 $Ez=$eventid if $eventid eq "AD";
64 $Ez=$eventid if $eventid eq "IR";
65 if ($eventid eq "CS" or $eventid eq "CA" or $eventid eq "ME")
66 {
67 if($bir >0.9)
68 {
69 $Ez = "CS";
70 }
71 else
72 {
73 $Ez = "CA";
74 }
75 }
76
77 }
78 else
79 {
80 #$Ez="declined";
81 $Ez = "na";
82 }
83 if( ($j12<$JunctionCut or $j23<$JunctionCut) and $j13 <$JunctionCut)
84 {
85 $Ez = "na";
86 }
87 print $print,"\t",$stat1,"\t",$stat2,"\t",$stat3,"\t",$Ez,"\n";
88 }
89 close(IN);
90 ####################################################################
91
92 sub cutoff{
93 my($s,$c,$E,%cutoff)=@_;
94 my @Exlen=@$E;
95 if($c eq 'NA'){return('NA')}
96 my $range=$Exlen[$#Exlen];
97 foreach my $l(@Exlen){if($s<$l){$range=$l;last}}
98 if($c<$cutoff{$range}){return('no')}
99 return('yes')
100 }
101
102
103 sub contain{
104 my ($a,@a)=@_;
105 foreach(@a){if($a eq $_){return(1)}}
106 return(0)
107 }
108