Mercurial > repos > xuebing > sharplab_seq_motif

comparison splicesitescore/score5.pl @ 15:0e221dbd17b2 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Sat, 31 Mar 2012 08:53:06 -0400
parents
children
comparison
equal deleted inserted replaced
14:d1f0f85ee5bc 15:0e221dbd17b2
1 use strict;
2
3
4 my $inputfile = $ARGV[0];
5 my $usemaxent = 1;
6
7 my $modelpath = "/Users/xuebing/galaxy-dist/tools/mytools/splicesitescore/";
8 my %me2x5 = &makescorematrix($modelpath.'me2x5');
9 my %seq = &makesequencematrix($modelpath.'splicemodels/splice5sequences');
10
11 my %bgd;
12 $bgd{'A'} = 0.27;
13 $bgd{'C'} = 0.23;
14 $bgd{'G'} = 0.23;
15 $bgd{'T'} = 0.27;
16
17
18
19 open (FILE,"<$inputfile") || die "can't open!\n";
20
21 while(<FILE>) {
22 chomp;
23 if (/^\s*$/) { #discard blank lines;
24 next;
25 }
26 elsif (/^>/) { #discard comment lines;
27 print $_."\t";
28 next;
29 }
30 elsif (/[NQWERYUIOPLKJHFDSZXVBM]/) {
31 next;
32 }
33 else {
34 $_ =~ s/\cM//g; #gets rid of carriage return
35 my $str = $_;
36 print $str."\t";
37 $str = uc($str);
38 if ($usemaxent) {
39 print sprintf("%.2f",&log2(&scoreconsensus($str)*$me2x5{$seq{&getrest($str)}}))."\n";
40 }
41 }
42 }
43
44
45 sub makesequencematrix{
46 my $file = shift;
47 my %matrix;my $n=0;
48 open(SCOREF, $file) || die "Can't open $file!\n";
49 while(<SCOREF>) {
50 chomp;
51 $_=~ s/\s//;
52 $matrix{$_} = $n;
53 $n++;
54 }
55 close(SCOREF);
56 return %matrix;
57 }
58 sub makescorematrix{
59 my $file = shift;
60 my %matrix;my $n=0;
61 open(SCOREF, $file) || die "Can't open $file!\n";
62 while(<SCOREF>) {
63 chomp;
64 $_=~ s/\s//;
65 $matrix{$n} = $_;
66 $n++;
67 }
68 close(SCOREF);
69 return %matrix;
70 }
71
72 sub getrest{
73 my $seq = shift;
74 my @seqa = split(//,uc($seq));
75 return $seqa[0].$seqa[1].$seqa[2].$seqa[5].$seqa[6].$seqa[7].$seqa[8];
76 }
77 sub scoreconsensus{
78 my $seq = shift;
79 my @seqa = split(//,uc($seq));
80 my %bgd;
81 $bgd{'A'} = 0.27;
82 $bgd{'C'} = 0.23;
83 $bgd{'G'} = 0.23;
84 $bgd{'T'} = 0.27;
85 my %cons1;
86 $cons1{'A'} = 0.004;
87 $cons1{'C'} = 0.0032;
88 $cons1{'G'} = 0.9896;
89 $cons1{'T'} = 0.0032;
90 my %cons2;
91 $cons2{'A'} = 0.0034;
92 $cons2{'C'} = 0.0039;
93 $cons2{'G'} = 0.0042;
94 $cons2{'T'} = 0.9884;
95 my $addscore = $cons1{$seqa[3]}*$cons2{$seqa[4]}/($bgd{$seqa[3]}*$bgd{$seqa[4]});
96 return $addscore;
97 }
98
99 sub log2{
100 my ($val) = @_;
101 return log($val)/log(2);
102 }