/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "rav.p" */ #include /* rav: running average information curve Tom Schneider NCI/FCRDC Bldg 469. Room 144 P.O. Box B Frederick, MD 21702-1201 (301) 846-5581 (-5532 for messages) permanent email: toms@alum.mit.edu toms@ncifcrf.gov http://www-lecb.ncifcrf.gov/~toms/ National Cancer Institute Laboratory of Experimental and Computational Biology */ /* end of program */ /* begin module version */ #define version 1.22 /* of rav.p 1998 January 27 1998 January 26: added color key origin 1998 January 21 */ #define updateversion 1.20 /* defines lowest acceptable current parameter file */ /* end module version */ /* begin module describe.rav */ /* name rav: running average information curve synopsis rav(symvec: in, ravp: in, xyin: out, marks: out, colorkey: out output: out) files symvec: output of rsdata or alpro, contains an information curve. ravp: parameters to control the program. The file must contain the following parameters, one per line: parameterversion: The version number of the program. This allows the user to be warned if an old parameter file is used. minwindowsize: (integer, bases) minimum size of windows to average. maxwindowsize: (integer, bases) maximum size of windows to average. stepwindowsize: (integer, bases) how to increment window size. colorX, colorY: (real, cm) x and y coordinates of lower left corner of color key for color to bit conversion. colorWidth, colorHeight: (real, cm) x and y size of color key (respectively) colorThick: (real, cm) thickness of lines in color key strip colorTransform: (integer) Defines how the PostScript colors are to be transformed. PostScript colors run from red through red, so the highest values have the same color as the lowest ones. This parameter determines how to transform the color scale. 0: no transform 1: transform as given by hue := colorAtransform*hue + colorBtransform This gives the user full control as a linear transformation. 2: cutoff: If the hue is less than colorAtransform/2 or larger than colorBtransform/3 then hue is set to 0.3; otherwise it is zero 0 (red). This gives a red line in the region between the two parameters. 3. Transform that gives yellow up to red: hue := 0.84*hue + 0.16 4: Roy G. Biv transform: the spectrum. (red, orange, yellow, green, blue, indigo, violet) where red is highest. The transform is hue := -magic*hue + magic; with magic = 0.85 colorAtransform: multiplicative parameter for colorTransform = 2 colorBtransform: additive parameter for colorTransform = 2 xyin: output of this program for plotting by xyplo. data columns: 1: position of window center 2: window size 3: average of information in window (bits) 4: standard deviation of information (bits) 5: color hue (0 to 1) 6: color saturation (0 to 1) 7: color brightness (0 to 1) 8: xsize 9: ysize marks: output of this program for inclusion with a sequence logo created by makelogo. This is the best way to display the data! colorkey: a PostScript file that creates a key showing the relationship between bits and colors. When one runs xyin through the xyplo program the color shows the degree of (averaged) conservation; this key can be inserted into the xyout graphic by including it in your xyplom file. output: messages to the user description The program rav creates a Running AVerage information curve starting from the symvec format produced by alpro or rseq. This is an improvement over the winfo program because winfo only used the rsdata format from rseq and therefore cannot handle sequences that come from alpro (which are in symvec format). In addition, because results are reported in the marks file, this program supercedes winfo. Windows coordinates are reported at the center of the window. The marks file is set up so that the the average value is shown by a rectangle, ranging one standard deviation above and below the mean. Note that for odd window sizes the boxes are placed surrounding a base in the sequence logo, while odd ones are offset between bases. You can check that the program is working correctly by trying window sizes of 1 (which should match the heights of the logo and have boxes that match the error bars); size 2 (which should be between bases and start to show averaging) etc. examples ravp: 1.20 version of rav that this parameter file is designed for. 1 minimum window size to average. 50 maximum window size to average. 5 increment of window size. 18.3 10 colorX, colorY (cm): coordinates of lower left corner of key 0.5 8 colorWidth, colorHeight (cm): size of key 0.2 colorThick (cm): thickness of lines in color key strip 4 colorTransform 0:none,1:control,2:cutoff,3:.84hue+.16,4:-.85hue+.85 0.84 colorAtransform multiplicative color transform parameter 0.16 colorBtransform additive color transform parameter xyplop: 2 6 zerox zeroy graph coordinate center zx 0 3000 zx min max (character, real, real) if zx='x' then set xaxis zy 0 100 zy min max (character, real, real) if zy='y' then set yaxis 10 10 1 1 xinterval yinterval xsubintervals ysubintervals: axis intervals 6 4 xwidth ywidth width of numbers in characters 0 0 xdecimal ydecimal number of decimal places 16 16 xsize ysize size of axes in cm position window size a zc 'c' crosshairs, axXyYnN n 2 zxl base if zxl='l' then make x axis log to the given base n 2 zyl base if zyl='l' then make y axis log to the given base ********************************************************************* 1 2 xcolumn ycolumn columns of xyin that determine plot location 0 symbol column the xyin column to read symbols from 8 9 xscolumn yscolumn columns of xyin that determine the symbol size 5 6 7 hue saturation brightness columns for color manipulation ********************************************************************* R symbol-to-plot c(circle)bd(dotted box)x+Ifgpr(rectangle) 0 symbol-flag character in xyin that indicates that this symbol -1.0 symbol sizex side in cm on the x axis of the symbol. -1.0 symbol sizey as for the x axis, get size from yscolumn n no connection (example for connection is c- 0.125 for dashed 0.125 cm) n 0.125 linetype size linetype l.-in and size of dashes or dots ********************************************************************* . **** version 8.50 of xyplop: uses cm for distances documentation see also winfo.p, rseq.p, xyplo.p, makelogo.p, struco.p, ravp, ravxyplop author Thomas Dana Schneider bugs technical notes Constant maxwin is the largest window size allowed. When the sequence logo runs over multiple lines, the marks from rav have to be triggered AFTER the letter. If this is not done, the box appears on the previous line, and is missing from the first stack of the next line. The trick is to write the box backwards - right to left. The box is therefore triggered after the stack has been written and is always associated with the stack. This has the consequence that the box is written on top of the stack, but it is completely clean. dedicated to Pete Rogan */ /* end module describe.rav */ /* begin module rav.const */ #define maxwin 2000 /* the largest window size allowed. */ #define infowid 10 /* width of reported information */ #define infodec 5 /* decimals for reported information */ /* end module rav.const */ /* begin module rav.type */ typedef struct infopoint { /* information about a point in an alignment */ long position; /* position in the alignment */ long samples; /* number of sequences at this position */ double information; /* information in bits */ double variance; /* variance of information */ } infopoint; typedef struct infoarray { /* store some infopoints */ long symbols; /* number of symbols in symvec */ infopoint data[maxwin]; /* array of infopoints */ long minwindowsize; /* minimum size of window requested */ long maxwindowsize; /* maximum size of window requested */ long stepwindowsize; /* step size of window requested */ long windowsize; /* current size of window */ double halfwindowsize; /* half of current windowsize */ long current; /* current window position */ long last; /* the last window position */ double Rtotal; /* current total of information R in the window */ double Vtotal; /* current total of variance in the window */ /* color key controls */ double colorX; /* x coordinate of lower left of color key strip */ double colorY; /* y coordinate of lower left of color key strip */ double colorWidth; /* horizontal width of color key strip */ double colorHeight; /* vertical width of color key strip */ double colorThick; /* thickness of lines in color key strip */ long colorTransform; /* kind of color transform to do */ double colorAtransform, colorBtransform; /* parameters of the transform */ } infoarray; /* end module rav.type */ Static _TEXT symvec; /* file used by this program */ Static _TEXT ravp; /* file used by this program */ Static _TEXT xyin, marks, colorkey; /* file used by this program */ Static jmp_buf _JL1; /* begin module halt */ Static Void halt() { /* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. */ printf(" program halt.\n"); longjmp(_JL1, 1); } #define magic 0.85 /* a nice cutoff for colors */ /* end module halt version = 'delmod 6.16 84 mar 12 tds/gds'; */ /* begin module docolortransform */ Static Void docolortransform(colorTransform, colorAtransform, colorBtransform, hue) long colorTransform; double colorAtransform, colorBtransform, *hue; { /* transformations of the hue, which runs from 0 to 1 into another range which gives different colors. 0: no transform 1: transform as given by a hue := colorAtransform*hue + colorBtransform 2: cutoff: hue less than colorAtransform/2 or larger than colorBtransform/3 is 0.3 otherwise 0 (red) 3. Transform that gives yellow up to red: hue := 0.84*hue + 0.16 4: Roy G. Biv transform (red, orange, yellow, green, blue, indigo, violet) where red is highest. hue := -magic*hue + magic; with magic = 0.85 */ switch (colorTransform) { case 0: /* blank case */ break; case 1: *hue = colorAtransform * *hue + colorBtransform; break; case 2: if (*hue < colorAtransform / 2) *hue = 0.3; else if (*hue > colorBtransform / 2) *hue = 0.3; else *hue = 0.0; break; case 3: *hue = 0.84 * *hue + 0.16; break; case 4: *hue = magic - magic * *hue; break; } } #undef magic #define ticjump 0.1 /* distance to jump to the right of the color strip to get to the tic mark, cm */ #define ticlength 0.2 /* length of tic mark, cm */ /* end module docolortransform */ /* begin module makecolorkey */ Static Void makecolorkey(psout, colorX, colorY, colorWidth, colorHeight, colorThick, colorTransform, colorAtransform, colorBtransform) _TEXT *psout; double colorX, colorY, colorWidth, colorHeight, colorThick; long colorTransform; double colorAtransform, colorBtransform; { /* make a color key bar with lower left hand coordinate at (colorX, colorY) of height colorHeight of width colorWidth, and with lines of thickness colorThick. The bit range is from 0 to 2 bits. */ double cmtopoints = 72 / 2.54; /* conversion factor from cm to points. */ long hueinteger; /* used to scan the hues */ double hue; /* a number between 0 and 1 that defines a color shade */ double y; /* current y coordinate of color strip in points */ if (*psout->name != '\0') { if (psout->f != NULL) psout->f = freopen(psout->name, "w", psout->f); else psout->f = fopen(psout->name, "w"); } else { if (psout->f != NULL) rewind(psout->f); else psout->f = tmpfile(); } if (psout->f == NULL) _EscIO2(FileNotFound, psout->name); SETUPBUF(psout->f, Char); fprintf(psout->f, "%%! colorkey\n"); /* 72 points/inch * 2.54 inch/cm */ /* define a new PostScript function for color: */ fprintf(psout->f, "/Times-Roman findfont 10 scalefont setfont\n"); fprintf(psout->f, "/c {1 1 sethsbcolor} def\n"); fprintf(psout->f, "%% define the colors\n"); fprintf(psout->f, "gsave\n"); fprintf(psout->f, "%8.5f setlinewidth\n", colorThick * cmtopoints); for (hueinteger = 0; hueinteger <= 100; hueinteger++) { hue = hueinteger / 100.0; /* at this point hue is in cm */ y = (colorHeight * hue + colorY) * cmtopoints; docolortransform(colorTransform, colorAtransform, colorBtransform, &hue); fprintf(psout->f, "%8.5f c", hue); fprintf(psout->f, " %8.5f %8.5f moveto %8.5f 0 rlineto stroke\n", colorX * cmtopoints, y, colorWidth * cmtopoints); if (hueinteger % 5 == 0) { /* tic marks */ fprintf(psout->f, "gsave\n"); fprintf(psout->f, "1 setlinewidth\n"); fprintf(psout->f, "1 0 0 sethsbcolor\n"); /* black */ fprintf(psout->f, " %8.5f %8.5f moveto\n", (colorX + colorWidth + ticjump) * cmtopoints, y); fprintf(psout->f, "currentpoint translate 0 0 moveto\n"); fprintf(psout->f, "%8.5f 0 rlineto currentpoint stroke\n", ticlength * cmtopoints); fprintf(psout->f, "translate 0 0 moveto\n"); if (hueinteger % 50 == 0) { /* numbers */ fprintf(psout->f, "(%2ld\n", (long)floor(hueinteger / 50.0 + 0.5)); if (hueinteger == 100) fprintf(psout->f, " bits\n"); fprintf(psout->f, ") show\n"); } /* p2c: rav.p, line 359: * Note: Using % for possibly-negative arguments [317] */ fprintf(psout->f, "grestore\n"); fprintf(psout->f, "%8.5f setlinewidth\n", colorThick * cmtopoints); } /* p2c: rav.p, line 350: * Note: Using % for possibly-negative arguments [317] */ } fprintf(psout->f, "grestore\n"); } #undef ticjump #undef ticlength /* end module makecolorkey */ /* begin module copyaline */ Static Void copyaline(fin, fout) _TEXT *fin, *fout; { /* copy a line from file fin to file fout */ while (!P_eoln(fin->f)) { putc(P_peek(fin->f), fout->f); getc(fin->f); } fscanf(fin->f, "%*[^\n]"); getc(fin->f); putc('\n', fout->f); } /* copyaline */ /* end module copyaline */ /* begin module rav.setheaders */ Static Void setheaders(symvec, xyin, marks, store) _TEXT *symvec, *xyin, *marks; infoarray *store; { /* read header information from symvec and transfer it to xyin and marks. Obtain number of symbols from symvec and report it too. */ if (*symvec->name != '\0') { if (symvec->f != NULL) symvec->f = freopen(symvec->name, "r", symvec->f); else symvec->f = fopen(symvec->name, "r"); } else rewind(symvec->f); if (symvec->f == NULL) _EscIO2(FileNotFound, symvec->name); RESETBUF(symvec->f, Char); if (*xyin->name != '\0') { if (xyin->f != NULL) xyin->f = freopen(xyin->name, "w", xyin->f); else xyin->f = fopen(xyin->name, "w"); } else { if (xyin->f != NULL) rewind(xyin->f); else xyin->f = tmpfile(); } if (xyin->f == NULL) _EscIO2(FileNotFound, xyin->name); SETUPBUF(xyin->f, Char); fprintf(xyin->f, "* rav %4.2f\n", version); while (P_peek(symvec->f) == '*') copyaline(symvec, xyin); if (*symvec->name != '\0') { if (symvec->f != NULL) symvec->f = freopen(symvec->name, "r", symvec->f); else symvec->f = fopen(symvec->name, "r"); } else rewind(symvec->f); if (symvec->f == NULL) _EscIO2(FileNotFound, symvec->name); RESETBUF(symvec->f, Char); if (*marks->name != '\0') { if (marks->f != NULL) marks->f = freopen(marks->name, "w", marks->f); else marks->f = fopen(marks->name, "w"); } else { if (marks->f != NULL) rewind(marks->f); else marks->f = tmpfile(); } if (marks->f == NULL) _EscIO2(FileNotFound, marks->name); SETUPBUF(marks->f, Char); fprintf(marks->f, "* rav %4.2f\n", version); while (P_peek(symvec->f) == '*') copyaline(symvec, marks); fscanf(symvec->f, "%ld%*[^\n]", &store->symbols); getc(symvec->f); fprintf(xyin->f, "* number of symbols = %ld\n", store->symbols); fprintf(marks->f, "* number of symbols = %ld\n", store->symbols); fprintf(xyin->f, "* minimum window size = %ld\n", store->minwindowsize); fprintf(xyin->f, "* maximum window size = %ld\n", store->maxwindowsize); fprintf(marks->f, "* minimum window size = %ld\n", store->minwindowsize); fprintf(marks->f, "* maximum window size = %ld\n", store->maxwindowsize); fprintf(xyin->f, "*\n"); fprintf(xyin->f, "* data columns:\n"); fprintf(xyin->f, "* 1: position of window center\n"); fprintf(xyin->f, "* 2: window size\n"); fprintf(xyin->f, "* 3: average of information in window (bits)\n"); fprintf(xyin->f, "* 4: standard deviation of information (bits)\n"); fprintf(xyin->f, "* 5: color hue (0 to 1)\n"); fprintf(xyin->f, "* 6: color saturation (0 to 1)\n"); fprintf(xyin->f, "* 7: color brightness (0 to 1)\n"); fprintf(xyin->f, "* 8: xsize\n"); fprintf(xyin->f, "* 9: ysize\n"); } Local Void die() { /* die at unexpected end of file */ printf("unexpected end of symvec file\n"); halt(); } /* end module rav.setheaders */ /* begin module rav.readoneposition */ Static Void readoneposition(symvec, i, symbols) _TEXT *symvec; infopoint *i; long symbols; { /* read the symvec information for one position into i, skip the symbols */ long s; /* index to symbols */ if (BUFEOF(symvec->f)) die(); fscanf(symvec->f, "%ld%ld%lg%lg%*[^\n]", &i->position, &i->samples, &i->information, &i->variance); getc(symvec->f); /* skip symbol data */ for (s = 1; s <= symbols; s++) { if (BUFEOF(symvec->f)) die(); fscanf(symvec->f, "%*[^\n]"); getc(symvec->f); } } /* end module rav.readoneposition */ /* begin module rav.reportwindow */ Static Void reportwindow(xyin, marks, store) _TEXT *xyin, *marks; infoarray *store; { /* report the current window data */ double hue; /* a number between 0 and 1 that defines a color shade */ double Rav; /* information average */ double Vav; /* information variance average */ double stdev; /* information standard deviation */ double windowcenter; /* center of the current window */ double saturation = 1.0; /* color saturation */ double brightness = 1.0; /* color brightness */ double xsize = 1.0; /* x axis box size */ double ysize; /* y axis box size */ ysize = store->stepwindowsize; windowcenter = store->data[store->current - 1].position - store->halfwindowsize; windowcenter += 0.5; /* writeln(output,'halfwindowsize: ', halfwindowsize:10:5); writeln(output,'data[current].position: ', data[current].position:10:5); writeln(output,'windowcenter: ', windowcenter:10:5); */ Rav = store->Rtotal / store->windowsize; Vav = store->Vtotal / store->windowsize; /* unfortunately Vtotal can become 'negative zero' - protect against this stupidity */ if (store->Vtotal <= 0.0) stdev = 0.0; else stdev = sqrt(Vav); fprintf(xyin->f, "%*.*f", infowid, infodec, windowcenter); fprintf(xyin->f, "%*ld", infowid, store->windowsize); fprintf(xyin->f, " %*.*f", infowid, infodec, Rav); fprintf(xyin->f, " %*.*f", infowid, infodec, stdev); hue = Rav / 2.0; docolortransform(store->colorTransform, store->colorAtransform, store->colorBtransform, &hue); fprintf(xyin->f, " %*.*f", infowid, infodec, hue); fprintf(xyin->f, " %*.*f", infowid, infodec, saturation); fprintf(xyin->f, " %*.*f", infowid, infodec, brightness); fprintf(xyin->f, " %*.*f", infowid, infodec, xsize); fprintf(xyin->f, " %*.*f\n", infowid, infodec, ysize); /* This is the point where the mark is put in backwards */ fprintf(marks->f, "bs %*.*f %*.*f %*.*f %*.*f\n", infowid, infodec, windowcenter + 0.5, infowid, infodec, Rav - stdev, infowid, infodec, windowcenter - 0.5, infowid, infodec, Rav + stdev); } /* end module rav.reportwindow */ /* begin module rav.initialread */ Static Void initialread(symvec, store) _TEXT *symvec; infoarray *store; { /* read the first windowsize data points */ long w; /* index to window */ long FORLIM; if (*symvec->name != '\0') { if (symvec->f != NULL) symvec->f = freopen(symvec->name, "r", symvec->f); else symvec->f = fopen(symvec->name, "r"); } else rewind(symvec->f); if (symvec->f == NULL) _EscIO2(FileNotFound, symvec->name); RESETBUF(symvec->f, Char); while (P_peek(symvec->f) == '*') { /* skip header */ fscanf(symvec->f, "%*[^\n]"); getc(symvec->f); } fscanf(symvec->f, "%*[^\n]"); getc(symvec->f); /* skip digits */ store->halfwindowsize = store->windowsize / 2.0; store->current = 0; store->Rtotal = 0.0; store->Vtotal = 0.0; FORLIM = store->windowsize; for (w = 0; w < FORLIM; w++) { readoneposition(symvec, &store->data[w], store->symbols); store->current++; store->Rtotal += store->data[store->current - 1].information; store->Vtotal += store->data[store->current - 1].variance; } store->last = 1; /* last tracks the hind end of the window */ } /* end module rav.initialread */ /* begin module rav.doaposition */ Static Void doaposition(symvec, xyin, marks, store) _TEXT *symvec, *xyin, *marks; infoarray *store; { /* read a symvec position and compute the window */ store->current++; if (store->current > store->windowsize) store->current = 1; store->Rtotal -= store->data[store->current - 1].information; store->Vtotal -= store->data[store->current - 1].variance; readoneposition(symvec, &store->data[store->current - 1], store->symbols); store->Rtotal += store->data[store->current - 1].information; store->Vtotal += store->data[store->current - 1].variance; reportwindow(xyin, marks, store); } /* end module rav.doaposition */ /* begin module rav.readwindowsize */ Static Void readwindowsize(f, windowsize) _TEXT *f; long *windowsize; { /* read the window size from file f */ fscanf(f->f, "%ld%*[^\n]", windowsize); getc(f->f); if (*windowsize > maxwin) { printf("windowsize requested was %ld; ", *windowsize); printf("windowsize must not exceed %ld\n", (long)maxwin); halt(); } if (*windowsize >= 1) return; printf("windowsize requested was %ld; ", *windowsize); printf("windowsize must be positive\n"); halt(); } /* end module rav.readwindowsize */ /* begin module rav.themain */ Static Void themain(symvec, ravp, xyin, marks, colorkey) _TEXT *symvec, *ravp, *xyin, *marks, *colorkey; { /* the main procedure of the program */ double parameterversion; /* parameter version number */ infoarray store; /* storage of current window */ long w; /* current window size */ printf("rav %4.2f\n", version); if (*ravp->name != '\0') { if (ravp->f != NULL) ravp->f = freopen(ravp->name, "r", ravp->f); else ravp->f = fopen(ravp->name, "r"); } else rewind(ravp->f); if (ravp->f == NULL) _EscIO2(FileNotFound, ravp->name); RESETBUF(ravp->f, Char); fscanf(ravp->f, "%lg%*[^\n]", ¶meterversion); getc(ravp->f); if (parameterversion < updateversion) { printf("You have an old parameter file!\n"); halt(); } readwindowsize(ravp, &store.minwindowsize); readwindowsize(ravp, &store.maxwindowsize); fscanf(ravp->f, "%ld%*[^\n]", &store.stepwindowsize); getc(ravp->f); fscanf(ravp->f, "%lg%lg%*[^\n]", &store.colorX, &store.colorY); getc(ravp->f); fscanf(ravp->f, "%lg%lg%*[^\n]", &store.colorHeight, &store.colorWidth); getc(ravp->f); fscanf(ravp->f, "%lg%*[^\n]", &store.colorThick); getc(ravp->f); fscanf(ravp->f, "%ld%*[^\n]", &store.colorTransform); getc(ravp->f); fscanf(ravp->f, "%lg%*[^\n]", &store.colorAtransform); getc(ravp->f); fscanf(ravp->f, "%lg%*[^\n]", &store.colorBtransform); getc(ravp->f); makecolorkey(colorkey, store.colorX, store.colorY, store.colorHeight, store.colorWidth, store.colorThick, store.colorTransform, store.colorAtransform, store.colorBtransform); setheaders(symvec, xyin, marks, &store); w = store.minwindowsize; while (w <= store.maxwindowsize) { printf("working on window size %ld\n", w); store.windowsize = w; initialread(symvec, &store); reportwindow(xyin, marks, &store); while (!BUFEOF(symvec->f)) doaposition(symvec, xyin, marks, &store); w += store.stepwindowsize; } } /* end module rav.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; colorkey.f = NULL; strcpy(colorkey.name, "colorkey"); marks.f = NULL; strcpy(marks.name, "marks"); xyin.f = NULL; strcpy(xyin.name, "xyin"); ravp.f = NULL; strcpy(ravp.name, "ravp"); symvec.f = NULL; strcpy(symvec.name, "symvec"); themain(&symvec, &ravp, &xyin, &marks, &colorkey); _L1: if (symvec.f != NULL) fclose(symvec.f); if (ravp.f != NULL) fclose(ravp.f); if (xyin.f != NULL) fclose(xyin.f); if (marks.f != NULL) fclose(marks.f); if (colorkey.f != NULL) fclose(colorkey.f); exit(EXIT_SUCCESS); } /* End. */