By downloading this code you agree to the
Source Code Use License (PDF). |
{ version = 2.21; (* of diana.p 2013 Apr 10} (* begin module describe.diana *) (* name diana: dinucleotide analysis of an aligned book synopsis diana(book: in, inst: in, protseq: in, dianap: in, da: out, output: out); files book: the standard delila book that is to be analyzed inst: the instruction set that was used to make the book protseq: aligned sequences in the protseq format (see alpro.p). The sequences can only contain bases a, c, g, t. This file may be empty, in which case the book and inst are used. Likewise, if the book is empty, the protseq will be used. One of the two must be empty. dianap: Parameters to control the program, one per line: 0. The version of diana that this parameter file is designed for. If the program finds an old version, it will *upgrade* the dianap file. 1. fromparam and toparam (two integers), determine the from-to range over which to do the analysis. If dianap is empty, the range from book and inst is used. 2. fromprotseq and toprotseq (two integers) define the range of the protseq. They are only used if the protseq file is used. 3. dicontrol: The first character of the third line controls printing to da: 'd': just sets of 16 correlation data are printed (ie, only d in column 1 of the da file). 'i': just information about the correlation is printed (ie, only i in column 1 of the da file). 'b': both data and information printed 4. colorslope, colorconstant: two real numbers on the fourth line control the reported value of the normalized information content (column 13 below). This allows one to adjust the range of colors produced in a color plot. For PostScript colors, 0.84 and 0.16 respecitvely gives yellow ranging up to red. Once the Rnormalized has been calculated, it is replaced by the value: (Rnormalized times colorslope) added to colorconstant. 5. alignmenttype: the first character on the first line, defines how to align the pieces. See the alist program for the detailed logic. There are three choices, as in the alist program: 'f' (for 'first') then the sequences are always aligned by their first base. 'i' then the sequences are aligned by the delila instructions. If the inst file is empty, alignment is forced to the 'b' mode. 'b' (for 'internal') then the alignment is on the internal zero of the book's sequence. This option is to be used when "default coordinate zero" is used in the Delila instructions. 6. rangecontrol: if the first character on the first line is an 'r' then the next two tokens must be two real numbers rangelower and rangeupper that define the range of Rnormalized values to use. When Rnormalized is between (inclusive) rangelower and rangeupper, column 17 is set to 'INN' otherwise it is 'OUT'. This allows one to grep the range. 7. maxsequences: (integer) maximum number of sequences to read. 8. maxdistance: (integer) maximum number distance between positions to compute. By making this small the total computation is faster. da: the di-analysis file that contains the output triangular array. column 1: Tells what that row of data is. There are three choices: n : the column is a normal data element in the triangle. d : the column is an element on the diagonal of the triangle, where the two coordinates are equal. i : the column contains the information content correlating the two positions defined by columns 3 and 4. column 2: Tells what the dinucleotide is. There are 16 dinucleotides aa, ac, ag, at, ... , tt as well as `in' or `id', which denote columns that contain information content. `id' means that the information is for the diagonal, while `in' are off-diagonal. Combining "in" with "id" gives the entire information triangle. column 3: The position on the sequence that corresponds to the x coordinate on a Cartesian graph. column 4: The position on the sequence that corresponds to the y coordinate on a Cartesian graph. column 5: A column of constant 1's usable by xyplo (usually xscolumn). column 6: A column of constant 1's usable by xyplo (usually xscolumn). column 7: The number of data points at a position column 8: Rcorrelation (with small sample correction) (bits) The frequency of the dinucleotide in column 2 at position (column 3, column 4). If the column is an information column, then this is the information at that position. When the position is off diagonal, this is the correlation information Rcorrelation with small sample correction (ssc16). When the position is on diagonal it is Rsequence. (2008 Jul 4: "without small sample correction" was here, but it appears to be in the code as ssc4.) column 9: One minus the frequency (or 1 - column 8). If this is an information column, then this is the chi-square value. column 10: A column of constants usable by xyplo. If this is in an information row, then this column is the number of degrees of freedom at that position column 11: A column of constant 1's usable by xyplo (usually hue). column 12: A column of constant 0's usable by xyplo (usually saturation). column 13: Rnormalized = Rcorrelation/4 or 1-(Rsequence/2) Information from column 8 normalized by dividing by the maximum possible information (4 for non-diagonal, 2 for diagonal). This gives a number between 0 and 1. This number is then multiplied by the parameter colorslope and then added to colorconstant. column 14: 1-Rnormalized; 1 minus column 13 column 15: correlation coefficient for x to y on information (in or id) rows column 16: A column of constant 1's usable by xyplo (saturation alternative). column 17: INN if rangelower <= Rcorrelation <= rangeupper OUT otherwise output: error messages from the program description Diana goes through a book and looks at relationships of dinucleotide pairs within the sequences. The output of the program is in the form of a triangular array rather than the entire square of x to y relationships. This saves half of the calculations. For every pair of coordinates, the frequency of the dinucleotide pair is tabulated. The program also calculates the chi-square for the dinucleotide given the expected mono-nucleotide frequencies. The correlation information is calculated as the information in the dinucleotide frequencies less the information in each of the two mononucleotides. The output of the program may be sent into xyplo for graphics. Note the distinction between the `in' and `id' information columns. Information in a binding site is usually going to appear on the diagonal, so by making this distinction, one can eliminate the diagonal information peaks for statstical analysis with genhis. The program now has an error correction for small sample size. See Schneider1986 p. 427, eqn A6. s = 4 or s = 16. examples See xyplop.diana and xyplop.diana.color for example xyplop files. To create a xyin file when the di control is 'b', extract all lines from the da file that contain ' in '. To do this in Unix use: grep ' in ' da > xyin That will not include the diagonal, which contains Rsequence and so normally should not be compared to the correlations. To include the diagonal use: grep i da > xyin To sort the output files to find the most significant correlations, one can use the Unix commands: grep in da | sort +8 -9 -n -r | cat > da.sorted The -r puts the most significant ones first. documentation @article{Stephens.Schneider.Splice, author = "R. M. Stephens and T. D. Schneider", title = "Features of spliceosome evolution and function inferred from an analysis of the information at human splice sites", journal = "J. Mol. Biol.", volume = "228", pages = "1124--1136", year = "1992"} see also Explanation of Diana: http://alum.mit.edu/www/toms/diana.html Example dianap file: dianap Example xyplop files used in conjunction with this program: xyplop.diana, xyplop.diana.color The program to plot the output: xyplo.p Program to list aligned sequences: alist.p Program to graph histograms: genhis.p The alpro program defines the protseq format: alpro.p Related programs for plotting the data in three dimensions: da3d.p, da3drotate.p author R. Michael Stephens Modified by Tom Schneider to have user requested range and: faster calculation in main loop [1995 Jan 9], new input routine that reads protseq format [1995 Jan 9]. bugs It would be nice to allow a restricted list of pairs to be done. protseq has no definition of coordinates, and this forces the use of the fromprotseq and toprotseq. A definition for multiply aligned sequences with coordinates should be created. technical notes *) (* end module describe.diana *) {This manual page was created by makman 1.45}{created by htmlink 1.62}