|
6
|
1 use strict;
|
|
|
2
|
|
|
3
|
|
|
4 my $inputfile = $ARGV[0];
|
|
|
5 my $usemaxent = 1;
|
|
|
6
|
|
8
|
7 my $modelpath = "/home/wuxbl/rowley/tools/maxentscan/splicemodels/";
|
|
|
8 my @metables = &makemaxentscores($modelpath);
|
|
6
|
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{
|
|
8
|
54 my $dir = shift;
|
|
6
|
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 }
|
|
8
|
74
|
|
6
|
75 sub maxentscore{
|
|
|
76 my $seq = shift;
|
|
|
77 my $table_ref = shift;
|
|
|
78 my @metables = @$table_ref;
|
|
|
79 my @sc;
|
|
|
80 $sc[0] = $metables[0]{&hashseq(substr($seq,0,7))};
|
|
|
81 $sc[1] = $metables[1]{&hashseq(substr($seq,7,7))};
|
|
|
82 $sc[2] = $metables[2]{&hashseq(substr($seq,14,7))};
|
|
|
83 $sc[3] = $metables[3]{&hashseq(substr($seq,4,7))};
|
|
|
84 $sc[4] = $metables[4]{&hashseq(substr($seq,11,7))};
|
|
|
85 $sc[5] = $metables[5]{&hashseq(substr($seq,4,3))};
|
|
|
86 $sc[6] = $metables[6]{&hashseq(substr($seq,7,4))};
|
|
|
87 $sc[7] = $metables[7]{&hashseq(substr($seq,11,3))};
|
|
|
88 $sc[8] = $metables[8]{&hashseq(substr($seq,14,4))};
|
|
|
89 my $finalscore = $sc[0] * $sc[1] * $sc[2] * $sc[3] * $sc[4] / ($sc[5] * $sc[6] * $sc[7] * $sc[8]);
|
|
|
90 return $finalscore;
|
|
|
91 }
|
|
|
92
|
|
|
93
|
|
|
94
|
|
|
95 sub getrest{
|
|
|
96 my $seq = shift;
|
|
|
97 my $seq_noconsensus = substr($seq,0,18).substr($seq,20,3);
|
|
|
98 return $seq_noconsensus;
|
|
|
99 }
|
|
|
100
|
|
|
101 sub scoreconsensus{
|
|
|
102 my $seq = shift;
|
|
|
103 my @seqa = split(//,uc($seq));
|
|
|
104 my %bgd;
|
|
|
105 $bgd{'A'} = 0.27;
|
|
|
106 $bgd{'C'} = 0.23;
|
|
|
107 $bgd{'G'} = 0.23;
|
|
|
108 $bgd{'T'} = 0.27;
|
|
|
109 my %cons1;
|
|
|
110 $cons1{'A'} = 0.9903;
|
|
|
111 $cons1{'C'} = 0.0032;
|
|
|
112 $cons1{'G'} = 0.0034;
|
|
|
113 $cons1{'T'} = 0.0030;
|
|
|
114 my %cons2;
|
|
|
115 $cons2{'A'} = 0.0027;
|
|
|
116 $cons2{'C'} = 0.0037;
|
|
|
117 $cons2{'G'} = 0.9905;
|
|
|
118 $cons2{'T'} = 0.0030;
|
|
|
119 my $addscore = $cons1{$seqa[18]} * $cons2{$seqa[19]}/ ($bgd{$seqa[18]} * $bgd{$seqa[19]});
|
|
|
120 return $addscore;
|
|
|
121 }
|
|
|
122 sub log2{
|
|
|
123 my ($val) = @_;
|
|
|
124 return log($val)/log(2);
|
|
|
125 }
|
