|
6
|
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 }
|
