Mercurial > repos > xuebing > splicesite_max_entropy

comparison splicesite/splicesitescore/score3.pl~ @ 4:cf2c4930dc9c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Sun, 01 Apr 2012 01:48:03 -0400
parents
children
comparison
equal deleted inserted replaced
3:2abecae46b7d 4:cf2c4930dc9c
1 use strict;
2
3
4 my $inputfile = $ARGV[0];
5 my $usemaxent = 1;
6
7
8 my @metables = &makemaxentscores;
9
10 open (FILE,"<$inputfile") || die "can't open!\n";
11
12 while(<FILE>) {
13 chomp;
14 if (/^\s*$/) { #discard blank lines;
15 next;
16 }
17 elsif (/^>/) { #discard comment lines;
18 print $_."\t";
19 next;
20 }
21 elsif (/[NQWERYUIOPLKJHFDSZXVBM]/) {
22 next;
23 }
24 else {
25 $_ =~ s/\cM//g; #gets rid of carriage return
26 my $str = $_;
27 print $str."\t";
28 $str = uc($str);
29 if ($usemaxent) {
30 print sprintf("%.2f",&log2(&scoreconsensus($str)*&maxentscore(&getrest($str),\@metables)))."\n";
31 }
32 }
33 }
34
35 sub hashseq{
36 #returns hash of sequence in base 4
37 # &hashseq('CAGAAGT') returns 4619
38 my $seq = shift;
39 $seq = uc($seq);
40 $seq =~ tr/ACGT/0123/;
41 my @seqa = split(//,$seq);
42 my $sum = 0;
43 my $len = length($seq);
44 my @four = (1,4,16,64,256,1024,4096,16384);
45 my $i=0;
46 while ($i<$len) {
47 $sum+= $seqa[$i] * $four[$len - $i -1] ;
48 $i++;
49 }
50 return $sum;
51 }
52
53 sub makemaxentscores{
54 my $dir = "/Users/xuebing/galaxy-dist/tools/mytools/splicesitescore/splicemodels/";
55 my @list = ('me2x3acc1','me2x3acc2','me2x3acc3','me2x3acc4',
56 'me2x3acc5','me2x3acc6','me2x3acc7','me2x3acc8','me2x3acc9');
57 my @metables;
58 my $num = 0 ;
59 foreach my $file (@list) {
60 my $n = 0;
61 open (SCOREF,"<".$dir.$file) || die "Can't open $file!\n";
62 while(<SCOREF>) {
63 chomp;
64 $_=~ s/\s//;
65 $metables[$num]{$n} = $_;
66 $n++;
67 }
68 close(SCOREF);
69 #print STDERR $file."\t".$num."\t".$n."\n";
70 $num++;
71 }
72 return @metables;
73 }
74 sub makewmmscores{
75 my $dir = "/bionet/geneyeo_essentials/MaxEntropy/webserver/splicemodels/";
76 my @list = ('me1s0acc1','me1s0acc2','me1s0acc3','me1s0acc4',
77 'me1s0acc5','me1s0acc6','me1s0acc7','me1s0acc8','me1s0acc9');
78 my @metables;
79 my $num = 0 ;
80 foreach my $file (@list) {
81 my $n = 0;
82 open (SCOREF,"<".$dir.$file) || die "Can't open $file!\n";
83 while(<SCOREF>) {
84 chomp;
85 $_=~ s/\s//;
86 $metables[$num]{$n} = $_;
87 $n++;
88 }
89 close(SCOREF);
90 #print STDERR $file."\t".$num."\t".$n."\n";
91 $num++;
92 }
93 return @metables;
94 }
95 sub makemmscores{
96 my $dir = "/bionet/geneyeo_essentials/MaxEntropy/webserver/splicemodels/";
97 my @list = ('me2s0acc1','me2s0acc2','me2s0acc3','me2s0acc4',
98 'me2s0acc5','me2s0acc6','me2s0acc7','me2s0acc8','me2s0acc9');
99 my @metables;
100 my $num = 0 ;
101 foreach my $file (@list) {
102 my $n = 0;
103 open (SCOREF,"<".$dir.$file) || die "Can't open $file!\n";
104 while(<SCOREF>) {
105 chomp;
106 $_=~ s/\s//;
107 $metables[$num]{$n} = $_;
108 $n++;
109 }
110 close(SCOREF);
111 #print STDERR $file."\t".$num."\t".$n."\n";
112 $num++;
113 }
114 return @metables;
115 }
116 sub maxentscore{
117 my $seq = shift;
118 my $table_ref = shift;
119 my @metables = @$table_ref;
120 my @sc;
121 $sc[0] = $metables[0]{&hashseq(substr($seq,0,7))};
122 $sc[1] = $metables[1]{&hashseq(substr($seq,7,7))};
123 $sc[2] = $metables[2]{&hashseq(substr($seq,14,7))};
124 $sc[3] = $metables[3]{&hashseq(substr($seq,4,7))};
125 $sc[4] = $metables[4]{&hashseq(substr($seq,11,7))};
126 $sc[5] = $metables[5]{&hashseq(substr($seq,4,3))};
127 $sc[6] = $metables[6]{&hashseq(substr($seq,7,4))};
128 $sc[7] = $metables[7]{&hashseq(substr($seq,11,3))};
129 $sc[8] = $metables[8]{&hashseq(substr($seq,14,4))};
130 my $finalscore = $sc[0] * $sc[1] * $sc[2] * $sc[3] * $sc[4] / ($sc[5] * $sc[6] * $sc[7] * $sc[8]);
131 return $finalscore;
132 }
133
134
135
136 sub getrest{
137 my $seq = shift;
138 my $seq_noconsensus = substr($seq,0,18).substr($seq,20,3);
139 return $seq_noconsensus;
140 }
141
142 sub scoreconsensus{
143 my $seq = shift;
144 my @seqa = split(//,uc($seq));
145 my %bgd;
146 $bgd{'A'} = 0.27;
147 $bgd{'C'} = 0.23;
148 $bgd{'G'} = 0.23;
149 $bgd{'T'} = 0.27;
150 my %cons1;
151 $cons1{'A'} = 0.9903;
152 $cons1{'C'} = 0.0032;
153 $cons1{'G'} = 0.0034;
154 $cons1{'T'} = 0.0030;
155 my %cons2;
156 $cons2{'A'} = 0.0027;
157 $cons2{'C'} = 0.0037;
158 $cons2{'G'} = 0.9905;
159 $cons2{'T'} = 0.0030;
160 my $addscore = $cons1{$seqa[18]} * $cons2{$seqa[19]}/ ($bgd{$seqa[18]} * $bgd{$seqa[19]});
161 return $addscore;
162 }
163 sub log2{
164 my ($val) = @_;
165 return log($val)/log(2);
166 }