/* Ethan Wolf 1995 */ #include #include #include "scatter_ps.h" /*** To print a point do ". NEWCOIL_SCORE PAIRCOIL_SCORE\n" **/ /*** Then end list of all points with "showpage"**/ /*** The parameters are: PairCoil likelihood linear region has equation like(x)=ax+b where the linear region runs from likelihood lower_like to upper_like. For computing NewCoil likelihoods we use gaussian approximations for their scores on the negatives in genbnk with parameter mean-g and sd-g and for the positives with parameter mean-cc and sd-cc. Then their likelihood is gauss-cc(x)/[gauss-cc(x) + 30*gauss-g(x)] (assuming that negatives are 30 times more prevalent than coils. ******/ void scatter_trailer(FILE *scatter) { fprintf(scatter,"showpage\n"); } /** Read in scores output by "scboth -v" which are in format: "seq# stock_score paircoil_score seq_code seq_name". Then output them as points in scatter plot as "stock_score paircoil_score ." **/ void scatter_points(FILE *scatter, FILE *scores) { char line[300]; double stock_score, paircoil_score; while (fgets(line,300,scores)) { sscanf(line,"%*d %lf %lf", &stock_score, &paircoil_score); fprintf(scatter,"%lf %lf .\n", stock_score, paircoil_score); } } void scatter_header(FILE *scatter, double a, double b, double lower_like, double upper_like) { fprintf(scatter,"%%!PS-Adobe-\n"); fprintf(scatter,"%%%%Center and scale the figure\n"); fprintf(scatter,"8.5 36 mul 11 36 mul translate\n"); fprintf(scatter,"0.8 0.8 scale\n"); fprintf(scatter,"-8.5 36 mul -11 36 mul translate\n"); fprintf(scatter,"\n%%%%The Figure Title.\n"); fprintf(scatter,"/name () def\n"); fprintf(scatter,"\n%%For printing on vertical and horizontal axis.\n"); fprintf(scatter,"/rshow {dup stringwidth pop neg 0 rmoveto show} def\n"); fprintf(scatter,"/cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n"); fprintf(scatter,"/label {/Times-Roman findfont 20 scalefont setfont cshow} def\n"); fprintf(scatter,"/ttlabel {/Courier findfont 20 scalefont setfont cshow} def\n"); fprintf(scatter,"%%%%Center and print the Title.\n"); fprintf(scatter,"4.5 72 mul 9.5 72 mul moveto name label\n"); fprintf(scatter,"\n/adjust {ymin sub 10 mul exch xmin sub 10 mul exch} def\n"); fprintf(scatter,"\n%%%%The following is used for printing out the scatter plot points.\n"); /** Unnecessary Stuff?? fprintf(scatter,"/choose { pop true } def %% choose will be redefined for selective plotting\n"); fprintf(scatter,"/. {choose {adjust .25 0 360 arc fill} {pop pop true} ifelse} def\n"); **/ fprintf(scatter,"/. {1 index 1 index ymin gt exch xmin gt and {adjust .9 0 360 arc fill} {pop pop}ifelse} def\n\n"); fprintf(scatter,"/xmin -2 def\n"); fprintf(scatter,"/ymin -26 def\n"); fprintf(scatter,"/xmax 30 def\n"); fprintf(scatter,"/ymax 28 def\n"); fprintf(scatter,"/vline 1.35 ln 28 mul def\n"); fprintf(scatter,"/e 2.7182818284590452353 def\n\n"); fprintf(scatter,"%%%%For computing down*gaussian(x) with input x mean sd\n"); fprintf(scatter,"/tgauss { \n"); fprintf(scatter," /stddev exch def sub stddev div dup mul -2 div %% -[(x-mean)/sd]^2/2\n"); fprintf(scatter," e exch exp %% e^\n"); fprintf(scatter," 2.5066282746310002 div stddev div down mul\n"); fprintf(scatter,"} def\n\n"); stock_likelihood_setup(scatter); define_paircoil_likelihoodline(scatter, a, b, lower_like, upper_like); fprintf(scatter,"\n/axes { %% Do the axes\n"); fprintf(scatter," 90 100 translate\n"); fprintf(scatter,"%%%% Horizontal\n"); fprintf(scatter," 0.5 setlinewidth\n"); fprintf(scatter," xmin ymin adjust moveto xmax ymin adjust lineto stroke\n"); fprintf(scatter," xmax ymin adjust moveto -8 4 rlineto stroke\n"); fprintf(scatter," xmax ymin adjust moveto -8 -4 rlineto stroke\n"); fprintf(scatter," vline ymin adjust moveto 0 -66 rmoveto\n"); fprintf(scatter," (NEWCOILS) ttlabel\n"); fprintf(scatter,"%%%% Vertical\n"); fprintf(scatter," xmin ymin adjust moveto xmin ymax adjust lineto stroke\n"); fprintf(scatter," xmin ymax adjust moveto 4 -8 rlineto stroke\n"); fprintf(scatter," xmin ymax adjust moveto -4 -8 rlineto stroke\n"); fprintf(scatter," xmin hline adjust moveto -80 0 rmoveto\n"); fprintf(scatter," gsave 90 rotate (PairCoil) ttlabel grestore\n"); fprintf(scatter,"%%%% Vertical Ticks\n"); fprintf(scatter," -24 2 26 {xmin exch adjust moveto -3 0 rmoveto 6 0 rlineto stroke} for\n"); fprintf(scatter,"%%%% Horizontal Ticks\n"); fprintf(scatter," [10 11 12 13 14 15 16 18 20 22 24 26 28]\n"); fprintf(scatter," {10 div ln 28 mul ymin adjust moveto 0 -3 rmoveto 0 6 rlineto stroke} forall\n"); fprintf(scatter,"%%%% Horizontal Numbers\n"); fprintf(scatter," /Times-Roman findfont 15 scalefont setfont\n"); fprintf(scatter," [10 11 12 13 14 16 18 20 24 28]\n"); fprintf(scatter," {10 div dup ln 28 mul ymin adjust -22 add moveto\n"); fprintf(scatter," 5 string cvs cshow} forall\n"); stock_percentages(scatter); fprintf(scatter,"%%%% Vertical Numbers (paircoil scores)\n"); fprintf(scatter," -24 2 26 {dup xmin exch adjust moveto -10 -4 rmoveto\n"); fprintf(scatter," 4 string cvs rshow} for\n"); fprintf(scatter,"%%%% Vertical Likelihoods (percentages)\n"); fprintf(scatter," -10 2 18 {dup xmin exch adjust moveto -10 -4 rmoveto -60 0 rmoveto sfit show (%%) show} for\n"); fprintf(scatter," newpath\n"); fprintf(scatter,"} def %%Ends the Definition of the Axes\n\n\n"); fprintf(scatter,"0 setgray\n"); fprintf(scatter,"axes\n"); fprintf(scatter,"newpath\n"); fprintf(scatter,"vline dup\n"); fprintf(scatter,"ymin adjust moveto ymax adjust lineto\n"); fprintf(scatter,"stroke\n"); fprintf(scatter,"hline dup\n"); fprintf(scatter,"xmin exch adjust moveto xmax exch adjust lineto\n"); fprintf(scatter,"stroke\n"); fprintf(scatter,"newpath\n"); } /***This prints out post script stuff for the least square line ***/ /***between 20% and 80% determined by the rest of the likelihood program ***/ void define_paircoil_likelihoodline(FILE *scatter, double a, double b, double lower_like, double upper_like) { /** x sfit should ouput the likelihood score for x in the linear **/ /** region of the likelihood graph (between lower_like and upper_like) **/ /** where the curve is like=ax+b. **/ fprintf(scatter,"\n%%%%x sfit should ouput the likelihood score for x in the linear region of the\n"); fprintf(scatter,"%%%%likelihood graph (between lower_like and upper_like) where the curve is like=ax+b.\n"); fprintf(scatter,"/a %lf def\n", a); fprintf(scatter,"/b %lf def\n", b); fprintf(scatter,"%%%%hline is used to draw the 50 percent likelihood line for PairCoil\n"); fprintf(scatter,"/hline .5 b sub a div def\n"); fprintf(scatter,"/lower_like %.0lf def\n",lower_like); fprintf(scatter,"/upper_like %.0lf def\n",upper_like); fprintf(scatter,"/sfit {\n"); fprintf(scatter," a mul b add 100 mul\n"); fprintf(scatter," .5 add cvi %%Rounding to nearest percentage.\n"); fprintf(scatter," dup upper_like gt %%If likelihood > upper_like\n"); fprintf(scatter," {(>) show upper_like 3 string cvs}\n"); fprintf(scatter," %%then print > and make string upper_like.\n"); fprintf(scatter," {dup lower_like lt %%else if < lower_like \n"); fprintf(scatter," {(<) show lower_like 3 string cvs}\n"); fprintf(scatter," %%then print < and make string lower_like.\n"); fprintf(scatter," {( ) show 3 string cvs} ifelse %%else print a space and\n"); fprintf(scatter," %%make likelihood a string and convert to text\n"); fprintf(scatter," } ifelse %%ends else for when upper_like>likelihood(x)\n"); fprintf(scatter,"} def\n\n"); } void stock_likelihood_setup(FILE *scatter) { fprintf(scatter,"%%%%To compute the likelihood of STOCK use:\n"); fprintf(scatter,"%%%% gauss_positives(x)/[gauss_pos(x)+30*gauss_genbnk(x)]\n"); fprintf(scatter,"%%%%NewCoil parameters for computing likelihoods.\n"); fprintf(scatter,"/mean-g 0.77 def\n"); fprintf(scatter,"/sd-g 0.20 def\n"); fprintf(scatter,"/mean-cc 1.63 def\n"); fprintf(scatter,"/sd-cc 0.24 def\n"); fprintf(scatter,"/down 1 def\n\n"); fprintf(scatter,"/prob { dup mean-g sd-g tgauss 30 mul exch\n"); fprintf(scatter," mean-cc sd-cc tgauss dup 3 1 roll add div\n"); fprintf(scatter,"} def\n"); fprintf(scatter,"/prob-log {28 div 2.7182818284590452353 exch exp prob} def\n\n"); } void stock_percentages(FILE *scatter) { fprintf(scatter,"%%%% Stock Percentages at 1.2, 1.3, 1.4.\n"); fprintf(scatter," 1.1 ln 28 mul ymin adjust -44 add moveto (< 1%) cshow\n"); fprintf(scatter," 1.6 ln 28 mul ymin adjust -44 add moveto (>99%) cshow\n"); fprintf(scatter," [12 13 14]\n"); fprintf(scatter," {10 div ln 28 mul dup ymin adjust -44 add moveto\n"); fprintf(scatter," prob-log 100 mul cvi 12 string cvs cshow} forall\n\n"); }