5 #include "scatter_ps.h"
7 /*** To print a point do ". NEWCOIL_SCORE PAIRCOIL_SCORE\n" **/
8 /*** Then end list of all points with "showpage"**/
9 /*** The parameters are:
10 PairCoil likelihood linear region has equation like(x)=ax+b
11 where the linear region runs from likelihood lower_like to upper_like.
13 For computing NewCoil likelihoods we use gaussian approximations for their
14 scores on the negatives in genbnk with parameter mean-g and sd-g
15 and for the positives with parameter mean-cc and sd-cc.
16 Then their likelihood is gauss-cc(x)/[gauss-cc(x) + 30*gauss-g(x)]
17 (assuming that negatives are 30 times more prevalent than coils.
20 void scatter_trailer(FILE *scatter)
22 fprintf(scatter,"showpage\n");
28 /** Read in scores output by "scboth -v" which are in format:
29 "seq# stock_score paircoil_score seq_code seq_name".
30 Then output them as points in scatter plot as "stock_score paircoil_score ."
32 void scatter_points(FILE *scatter, FILE *scores)
35 double stock_score, paircoil_score;
37 while (fgets(line,300,scores)) {
38 sscanf(line,"%*d %lf %lf", &stock_score, &paircoil_score);
39 fprintf(scatter,"%lf %lf .\n", stock_score, paircoil_score);
45 void scatter_header(FILE *scatter, double a, double b, double lower_like, double upper_like)
47 fprintf(scatter,"%%!PS-Adobe-\n");
48 fprintf(scatter,"%%%%Center and scale the figure\n");
49 fprintf(scatter,"8.5 36 mul 11 36 mul translate\n");
50 fprintf(scatter,"0.8 0.8 scale\n");
51 fprintf(scatter,"-8.5 36 mul -11 36 mul translate\n");
53 fprintf(scatter,"\n%%%%The Figure Title.\n");
54 fprintf(scatter,"/name () def\n");
55 fprintf(scatter,"\n%%For printing on vertical and horizontal axis.\n");
56 fprintf(scatter,"/rshow {dup stringwidth pop neg 0 rmoveto show} def\n");
57 fprintf(scatter,"/cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");
58 fprintf(scatter,"/label {/Times-Roman findfont 20 scalefont setfont cshow} def\n");
59 fprintf(scatter,"/ttlabel {/Courier findfont 20 scalefont setfont cshow} def\n");
60 fprintf(scatter,"%%%%Center and print the Title.\n");
61 fprintf(scatter,"4.5 72 mul 9.5 72 mul moveto name label\n");
63 fprintf(scatter,"\n/adjust {ymin sub 10 mul exch xmin sub 10 mul exch} def\n");
66 fprintf(scatter,"\n%%%%The following is used for printing out the scatter plot points.\n");
67 /** Unnecessary Stuff??
68 fprintf(scatter,"/choose { pop true } def %% choose will be redefined for selective plotting\n");
69 fprintf(scatter,"/. {choose {adjust .25 0 360 arc fill} {pop pop true} ifelse} def\n");
71 fprintf(scatter,"/. {1 index 1 index ymin gt exch xmin gt and {adjust .9 0 360 arc fill} {pop pop}ifelse} def\n\n");
74 fprintf(scatter,"/xmin -2 def\n");
75 fprintf(scatter,"/ymin -26 def\n");
76 fprintf(scatter,"/xmax 30 def\n");
77 fprintf(scatter,"/ymax 28 def\n");
78 fprintf(scatter,"/vline 1.35 ln 28 mul def\n");
80 fprintf(scatter,"/e 2.7182818284590452353 def\n\n");
82 fprintf(scatter,"%%%%For computing down*gaussian(x) with input x mean sd\n");
83 fprintf(scatter,"/tgauss { \n");
84 fprintf(scatter," /stddev exch def sub stddev div dup mul -2 div %% -[(x-mean)/sd]^2/2\n");
85 fprintf(scatter," e exch exp %% e^<above>\n");
86 fprintf(scatter," 2.5066282746310002 div stddev div down mul\n");
87 fprintf(scatter,"} def\n\n");
90 stock_likelihood_setup(scatter);
93 define_paircoil_likelihoodline(scatter, a, b, lower_like, upper_like);
97 fprintf(scatter,"\n/axes { %% Do the axes\n");
98 fprintf(scatter," 90 100 translate\n");
99 fprintf(scatter,"%%%% Horizontal\n");
100 fprintf(scatter," 0.5 setlinewidth\n");
101 fprintf(scatter," xmin ymin adjust moveto xmax ymin adjust lineto stroke\n");
102 fprintf(scatter," xmax ymin adjust moveto -8 4 rlineto stroke\n");
103 fprintf(scatter," xmax ymin adjust moveto -8 -4 rlineto stroke\n");
104 fprintf(scatter," vline ymin adjust moveto 0 -66 rmoveto\n");
105 fprintf(scatter," (NEWCOILS) ttlabel\n");
106 fprintf(scatter,"%%%% Vertical\n");
107 fprintf(scatter," xmin ymin adjust moveto xmin ymax adjust lineto stroke\n");
108 fprintf(scatter," xmin ymax adjust moveto 4 -8 rlineto stroke\n");
109 fprintf(scatter," xmin ymax adjust moveto -4 -8 rlineto stroke\n");
110 fprintf(scatter," xmin hline adjust moveto -80 0 rmoveto\n");
111 fprintf(scatter," gsave 90 rotate (PairCoil) ttlabel grestore\n");
112 fprintf(scatter,"%%%% Vertical Ticks\n");
113 fprintf(scatter," -24 2 26 {xmin exch adjust moveto -3 0 rmoveto 6 0 rlineto stroke} for\n");
114 fprintf(scatter,"%%%% Horizontal Ticks\n");
115 fprintf(scatter," [10 11 12 13 14 15 16 18 20 22 24 26 28]\n");
116 fprintf(scatter," {10 div ln 28 mul ymin adjust moveto 0 -3 rmoveto 0 6 rlineto stroke} forall\n");
117 fprintf(scatter,"%%%% Horizontal Numbers\n");
118 fprintf(scatter," /Times-Roman findfont 15 scalefont setfont\n");
119 fprintf(scatter," [10 11 12 13 14 16 18 20 24 28]\n");
120 fprintf(scatter," {10 div dup ln 28 mul ymin adjust -22 add moveto\n");
121 fprintf(scatter," 5 string cvs cshow} forall\n");
124 stock_percentages(scatter);
127 fprintf(scatter,"%%%% Vertical Numbers (paircoil scores)\n");
128 fprintf(scatter," -24 2 26 {dup xmin exch adjust moveto -10 -4 rmoveto\n");
129 fprintf(scatter," 4 string cvs rshow} for\n");
130 fprintf(scatter,"%%%% Vertical Likelihoods (percentages)\n");
131 fprintf(scatter," -10 2 18 {dup xmin exch adjust moveto -10 -4 rmoveto -60 0 rmoveto sfit show (%%) show} for\n");
132 fprintf(scatter," newpath\n");
133 fprintf(scatter,"} def %%Ends the Definition of the Axes\n\n\n");
137 fprintf(scatter,"0 setgray\n");
138 fprintf(scatter,"axes\n");
139 fprintf(scatter,"newpath\n");
140 fprintf(scatter,"vline dup\n");
141 fprintf(scatter,"ymin adjust moveto ymax adjust lineto\n");
142 fprintf(scatter,"stroke\n");
143 fprintf(scatter,"hline dup\n");
144 fprintf(scatter,"xmin exch adjust moveto xmax exch adjust lineto\n");
145 fprintf(scatter,"stroke\n");
146 fprintf(scatter,"newpath\n");
154 /***This prints out post script stuff for the least square line ***/
155 /***between 20% and 80% determined by the rest of the likelihood program ***/
157 void define_paircoil_likelihoodline(FILE *scatter, double a, double b, double lower_like, double upper_like)
159 /** x sfit should ouput the likelihood score for x in the linear **/
160 /** region of the likelihood graph (between lower_like and upper_like) **/
161 /** where the curve is like=ax+b. **/
163 fprintf(scatter,"\n%%%%x sfit should ouput the likelihood score for x in the linear region of the\n");
164 fprintf(scatter,"%%%%likelihood graph (between lower_like and upper_like) where the curve is like=ax+b.\n");
166 fprintf(scatter,"/a %lf def\n", a);
167 fprintf(scatter,"/b %lf def\n", b);
168 fprintf(scatter,"%%%%hline is used to draw the 50 percent likelihood line for PairCoil\n");
169 fprintf(scatter,"/hline .5 b sub a div def\n");
170 fprintf(scatter,"/lower_like %.0lf def\n",lower_like);
171 fprintf(scatter,"/upper_like %.0lf def\n",upper_like);
172 fprintf(scatter,"/sfit {\n");
173 fprintf(scatter," a mul b add 100 mul\n");
174 fprintf(scatter," .5 add cvi %%Rounding to nearest percentage.\n");
175 fprintf(scatter," dup upper_like gt %%If likelihood > upper_like\n");
176 fprintf(scatter," {(>) show upper_like 3 string cvs}\n");
177 fprintf(scatter," %%then print > and make string upper_like.\n");
178 fprintf(scatter," {dup lower_like lt %%else if < lower_like \n");
179 fprintf(scatter," {(<) show lower_like 3 string cvs}\n");
180 fprintf(scatter," %%then print < and make string lower_like.\n");
181 fprintf(scatter," {( ) show 3 string cvs} ifelse %%else print a space and\n");
182 fprintf(scatter," %%make likelihood a string and convert to text\n");
183 fprintf(scatter," } ifelse %%ends else for when upper_like>likelihood(x)\n");
184 fprintf(scatter,"} def\n\n");
189 void stock_likelihood_setup(FILE *scatter)
191 fprintf(scatter,"%%%%To compute the likelihood of STOCK use:\n");
192 fprintf(scatter,"%%%% gauss_positives(x)/[gauss_pos(x)+30*gauss_genbnk(x)]\n");
194 fprintf(scatter,"%%%%NewCoil parameters for computing likelihoods.\n");
195 fprintf(scatter,"/mean-g 0.77 def\n");
196 fprintf(scatter,"/sd-g 0.20 def\n");
197 fprintf(scatter,"/mean-cc 1.63 def\n");
198 fprintf(scatter,"/sd-cc 0.24 def\n");
199 fprintf(scatter,"/down 1 def\n\n");
201 fprintf(scatter,"/prob { dup mean-g sd-g tgauss 30 mul exch\n");
202 fprintf(scatter," mean-cc sd-cc tgauss dup 3 1 roll add div\n");
203 fprintf(scatter,"} def\n");
204 fprintf(scatter,"/prob-log {28 div 2.7182818284590452353 exch exp prob} def\n\n");
208 void stock_percentages(FILE *scatter)
209 { fprintf(scatter,"%%%% Stock Percentages at 1.2, 1.3, 1.4.\n");
210 fprintf(scatter," 1.1 ln 28 mul ymin adjust -44 add moveto (< 1%) cshow\n");
211 fprintf(scatter," 1.6 ln 28 mul ymin adjust -44 add moveto (>99%) cshow\n");
212 fprintf(scatter," [12 13 14]\n");
213 fprintf(scatter," {10 div ln 28 mul dup ymin adjust -44 add moveto\n");
214 fprintf(scatter," prob-log 100 mul cvi 12 string cvs cshow} forall\n\n");