JP4069208B2 - Gene interaction estimation method, gene interaction estimation program, gene interaction estimation device, binding site estimation method, binding site estimation program, and binding site estimation device - Google Patents

Description

The present invention relates to a gene interaction estimation method, a gene interaction estimation program, a gene interaction estimation apparatus , and a The present invention relates to a binding site estimation method, a binding site estimation program, and a binding site estimation apparatus suitable for use in a gene interaction estimation method, a gene interaction estimation program, and a gene interaction estimation apparatus .

  Along with the progress of genomic science in recent years, it has become possible to simultaneously acquire a large number of gene expression data by a DNA microarray technology which is a technique for quantitatively measuring the expression intensity of a large number of genes in a certain environment. Research to elucidate the function of genes is also moving from the stage of analyzing individual genes to the stage of knowing how intracellular activities are performed by gene interactions.

  Here, for example, as shown in FIG. 27, the gene has a control region in the upstream portion, and binding of the control factor thereto promotes or suppresses gene expression. In addition, when a factor (protein) is generated by the expression of a gene, the factor affects the expression of other genes, but it is estimated how each gene affects other genes including itself. The problem is called a gene interaction estimation problem.

  As a method for estimating the gene interaction, a method based on model estimation using only expression data has been proposed. That is, a model that can best explain the expression data is estimated by using only the expression data obtained by the DNA microarray.

  First, as a method for estimating this kind of gene interaction, a Boolean network using Boolean algebra has been proposed (see, for example, Non-Patent Document 1). A Boolean network is a model like a logic circuit that operates synchronously. For example, in the case of the model shown in FIG. 28, one node (vertex) corresponds to one gene, and each vertex is “0 (not expressed)” or “1 (expressed)”. Take the state. Further, the state of each vertex changes synchronously at each unit time, and the gene expression control rule is expressed in the form of a Boolean function. In such a method, time series data of expression is given, and the problem of estimating a Boolean network that minimizes the mismatch between the data and the model is solved. However, in this method, there is a problem that quantitative expression intensity cannot be handled.

  In addition, as a method for estimating this kind of gene interaction, a method using a Bayesian network has been proposed (see, for example, Non-Patent Document 2). A probability model based on a graph structure in which random variables are represented by vertices and linked between variables having a dependent relationship such as a causal relationship or correlation is called a graphical model. Among them, for example, as shown in FIG. A direction in which the direction is the direction of causality and the path following this link does not circulate is called a Bayesian network. A Bayesian network represents a qualitative dependency relationship between random variables by a graph structure, and represents a quantitative dependency relationship between random variables by a conditional probability defined between the random variables. In such a technique, the expression intensity of each gene is used as a random variable, and the dependence between them is estimated as a Bayesian network. However, such a method has a problem that a plurality of models that can explain the same data exist.

  Furthermore, as a method for estimating this kind of gene interaction, a method using graphical Gaussian modeling has been proposed (see, for example, Non-Patent Document 3). The problem with clustering by correlation coefficient is that there is a correlation between direct interaction and a correlation through the expression of several other genes, or the expression of a common gene. It is indistinguishable from the case of the pseudo-correlation seen in the two affected genes. Here, assuming that the expression pattern in an environment with n genes is a sample from an n-dimensional normal distribution, the pseudo correlation is eliminated depending on whether or not the partial correlation coefficient is “0”. Can do. This is a technique using graphical Gaussian modeling. However, in such a method, when there is a gene whose expression behavior is very similar, the correlation coefficient matrix is not regular, and the partial correlation coefficient cannot be calculated. Therefore, in this method, clustering needs to be performed as preprocessing, and there is a problem that the dependency relationship can be estimated only in gene group units.

  Furthermore, as a method for estimating this kind of gene interaction, a method using a dynamic system model has been proposed (see, for example, Non-Patent Document 4). This method assumes a system model based on a differential equation or the like, formulates the parameter estimation as an optimization problem with the degree of coincidence with expression data as an objective function, and solves it by a genetic algorithm. However, this method has a problem that the calculation cost is high, thousands of gene networks cannot be solved at once, and processing such as division and hierarchization is required. In addition, in this method, there is a possibility of outputting a model that also explains noise, and there is a problem of how to eliminate noise.

  As described above, although the DNA microarray is an effective means for analyzing gene interactions, all of the methods described above use only expression data. There is a problem that only a few hundreds of experimental data can be obtained compared to the number of thousands of genes from the viewpoint of the problem and the experimental cost, and there is a problem that the relationship between the transcription factor and the transcription factor binding site cannot be obtained, The estimation accuracy is not high.

  In order to solve this problem, there is a method for deriving the relationship between transcription factors and genes controlled by these by combining expression data from DNA microarrays and analysis of patterns contained in control regions in the upstream part of genes. (For example, refer nonpatent literature 5, nonpatent literature 6, and patent literature 1 etc.).

  Specifically, Non-Patent Document 5 uses a series of expression data to examine the correlation between the expression pattern of transcription factors and the synthesis of the expression pattern of a gene having a common motif in the control region. A method for excluding the influence of a combination of transcription factors and estimating the relationship between the transcription factor and the binding site is disclosed. Non-Patent Document 6 discloses a technique of clustering genes based on the correlation between the expression pattern of transcription factors and the expression pattern of each gene, and then finding a motif. Further, Patent Document 1 discloses a method for estimating a control relationship between genes by a computer based on gene expression level data and gene sequence data. In other words, in this method, an early expression gene and a late expression gene are arbitrarily selected from time-series data of gene expression levels, a cross-correlation function for the two genes is obtained, and an uncorrelated test is performed to determine a null hypothesis. If is rejected, record the two genes. In this method, this procedure is repeated until the early expression gene and the late expression gene are completely selected.After that, a gene having a large number of consensus control sequences for the transcription factor is included in the upstream sequence of the transcription factor. The recorded early expression gene and late expression gene are collated, and when they match, the two genes are recorded, and it is estimated that these two genes have a regulatory relationship.

Takuya Akutsu, “Gene Network Estimation Algorithm”, Mathematical Sciences June 1999 issue, Science Inc., 1999, No. 432, p. 40-46 N. Friedman, M. Linial, I. Nachman and D. Pe’er, “Using Bayesian networks to analyze expression data”, Journal of Computational Biology, vol. 7, no. 3/4, 2000, p. 601-620 H. Toh, K. Horimoto, "Inference of a genetic network by a conbined approach of clusteranalysis and graphical Gaussian modeling", Bioinformatics, vol. 18, no. 2, 2002, p. 287-297 Masahiro Okamoto, “Estimation of gene interaction by S-system”, Genome Information Biology, Nakayama Shoten Co., Ltd., 2000, p. 165-188 K. Birnbaum, PN Benfey and DE Shasha, “cis Element / Transcription Factor Analysis (cis / TF): AMethod for Discovering Transcription Factor / cis Element Relationships”, Genome Res. 2001 September, vol. 11, no. 9, 2001 Year, p. 1567-1573 Z. Zhu, Y. Pilpel and G. M. Church, “Computational Identification of Transcription Factor Binding Sitesvia a Transcription-factor-centric Clustering (TFCC) Algorithm”, J. Mol. Biol. Vol. 318, 2002, p. 7181 JP 2003-141123 A

  By the way, although the expression level of each gene changes depending on various factors, the change in the expression intensity of the transcription factor is extremely small with respect to the background change such as noise. On the other hand, since the conventional methods described in Non-Patent Document 5, Non-Patent Document 6 and Patent Document 1 described above are all based on obtaining correlation between expression data, the influence of noise is reduced. There was a problem of receiving large.

  In addition, many transcription factors are known to be constantly expressed and controlled by the subsequent action, but there is a problem that conventional methods cannot be applied to such transcription factors. there were.

  Furthermore, Japanese Patent Application Laid-Open No. 2000-342257 discloses a technique that can clarify the transcriptional control of a gene by exploring the internal structure of an enhancer or promoter. That is, in this method, in estimating the binding site structure of a gene regulatory factor, an enhancer or promoter region to which a protein serving as a transcription factor in the upstream region or downstream region of the coding region encoding the gene binds. Set the gene for which the regulatory structure of the binding site of the regulatory factor is to be estimated, and the locus and other genes of the regulatory factor involved in the gene in the enhancer or promoter region or the binding site of the hypothetical regulatory factor A calculation model with the expression factor of the gene as a parameter is constructed, the transcription amount of the gene is calculated by the constructed calculation model, and the expression of the set gene obtained in the experiment is calculated using a parameter search algorithm. Search the parameters of the calculation model to obtain It has estimated the micro-structure of over or promoter. However, this method does not estimate the gene regulatory relationship itself.

The present invention has been made in view of such circumstances, and a gene interaction estimation method, a gene interaction estimation program capable of estimating a gene interaction with high accuracy while using a DNA microarray, and Provided are a gene interaction estimation device , a gene interaction estimation method, a gene interaction estimation program, and a binding site estimation method, a binding site estimation program, and a binding site estimation device suitable for use in the gene interaction estimation device. For the purpose.

  The applicant of the present application, as a result of earnest research on the gene interaction estimation problem to estimate the effect on other genes that are controlled by the factor when the gene is controlled by the factor, Fusing three independent estimation methods by filtering processing of 1, second filtering processing which is a completely new technique for estimating binding sites, and third filtering processing using gene expression data to which an experimental operation is added, Furthermore, the present invention has been made by devising an extremely novel technique as a specific process to be performed as the second filtering process.

That is, the gene interaction estimation method according to the present invention that achieves the above-described object is such that when gene a is controlled by factor F, gene interaction that estimates the effect on other gene b that is controlled by factor F is estimated. In the estimation method, an input step for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation and DNA base sequence data, which is data representing the base sequence, to a computer, and the computer The processor in (1) reads the DNA microarray disrupted strain data input in the input step, and the predetermined number of the upper predetermined number having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray disrupted strain data. and genes selection step of selecting a gene group B consisting of the genes b, the process Sa is the data of the gene a, for a pair of each gene b data of said gene group B selected in the genes selection step, based on the DNA base sequence data input in said input step When the control area of the gene a is taken out as a character string, a continuous area of a predetermined length l from each position i in the control area is defined as a window W _a [ i], and when the control region of each gene b of the gene group B is extracted as a character string, the position of the j-th character from the beginning is defined as a position j. cut out the continuous region length l as a window W _b [j], these windows W _{a [i],} the character string data contained in W b [j], these c When all the similarities between character strings of a predetermined length within a range deviated by a predetermined number at each position in the window are calculated, and the maximum value of the calculated similarity is defined as the similarity at that position . The window similarity that is the sum of the similarities for each position of the gene group B is calculated, and each window in the control region of the gene a is calculated based on the window similarity data including the window similarities of the genes b of the gene group B. A window similarity calculation step of calculating a maximum value Max _aB [i] and an average value Avg _aB [i] of the similarity of the gene group B to the position i, and the processor selected in the gene group selection step and data of each gene b of the genes B, the maximum value data and the average value a consisting of said window the maximum value Max aB calculated in the similarity calculation step _[i] reads the average value data consisting of g _{aB [i],} of the control region of the gene a, statistically specifically greater in the distribution of the windows similarities at the maximum value Max aB _[i] is the position i The gene with respect to a change in the position i of the control region of the gene a based on the window similarity data of each gene b of the gene group B calculated in the window similarity calculation step while obtaining a first region A second region that includes the peak position when the position that gives the maximum value when the maximum value Max _ab [i] of the window similarity of each gene b of group B is viewed is determined as the peak position; A region searching step for searching for a region R that is one region and the second region, and the gene selected by the processor in the gene group selection step Region data consisting of the data of each gene b of B, the window similarity data of each gene b of the gene group B calculated in the window similarity calculation step, and the region R obtained in the region search step DOO reads, among the genes b of the genes B, and each region R, there is a peak located in the area R, and obtains the genes B _R of similarity is a high gene b than the predetermined value The gene group extraction step to extract the expression intensity, and the expression intensity of the wild strain not added to the experimental operation based on the DNA microarray disrupted strain data input in the input step and the disrupted strain added the experimental operation using the change ratio, for each region R, among the genes of the genes B _R extracted by said genes extraction step largely expression intensity in the experimental manipulation is It is characterized in that it comprises only what is turned into a second group of genes extraction step of extracting as a gene group B _R ^F.

In addition, the gene interaction estimation program according to the present invention that achieves the above-described object can be executed by a computer that estimates the influence on the other gene b controlled by the factor F when the gene a is controlled by the factor F. In a gene interaction estimation program, the computer inputs DNA microarray disrupted strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data that is data representing the base sequence The input means reads the DNA microarray disrupted strain data input by the input means, and the upper predetermined number of the other genes b having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray disrupted strain data A gene group selection means for selecting a gene group B comprising the gene a For a pair of data and data of each gene b of the gene group B selected by the gene group selection means, a control region of the gene a is determined based on the DNA base sequence data input by the input means. When the position of the i-th character from the beginning is taken as position i when taken out as a character string, a continuous area of a predetermined length l is cut out from each position i in the control area as _a window W _a [i], and the gene When the control region of each gene b in group B is taken out as a character string, a continuous region having a predetermined length l from each position j in the control region is defined as a window W when the position of the jth character from the beginning is defined as a position j. excised as _b [j], these windows W _{a [i],} the character string data contained in W b [j], in these windows All the similarities between character strings of a predetermined length in a range that is shifted by a predetermined number at the front and back at each position, and each value in the window when the calculated maximum similarity is defined as the similarity of the position A window similarity that is a sum of similarities with respect to a position is calculated, and each position i in the control region of the gene a is calculated based on the window similarity data including the window similarity of each gene b of the gene group B. Window similarity calculation means for calculating the maximum value Max _aB [i] and the average value Avg _aB [i] of the similarity of the gene group B to each gene of the gene group B selected by the gene group selection means and b data the window maximum value data and the consisting calculated by the similarity calculating unit that said maximum value _Max aB [i] mean _Avg aB [i] or Becomes the average value is read and data of the control region of the gene a, obtaining a statistically specifically large first region the maximum value Max aB _[i] is in the distribution of the window similarity at the position i In addition, based on the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, each gene b of the gene group B with respect to a change in the position i of the control region of the gene a A second region that includes the peak position when a position that gives a maximum value when the maximum value Max _ab [i] of the window similarity is viewed is a peak position, and is the first region; Region search means for searching for the region R as the second region, data of each gene b of the gene group B selected by the gene group selection means; Reading the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means and the area data consisting of the area R obtained by the area search means, among the genes b, for each region R, there is a peak located in the area R, and, the gene group extracting means similarity extracted seeking genes B _R consisting of high gene b than the predetermined value, and Using the change ratio of the expression intensity of the wild strain not subjected to the experimental operation based on the DNA microarray disrupted strain data input by the input means and the disrupted strain subjected to the experimental operation, for each region R, among the genes of the genes B _R extracted by the genetic group extracting means, which expression intensity greatly changes in the experimental procedure It is characterized in that to function only as a second gene group extracting means for extracting a group of genes B _R ^F.

Furthermore, the gene interaction estimation apparatus according to the present invention that achieves the above-described object provides a gene interaction that estimates an influence on another gene b controlled by the factor F when the gene a is controlled by the factor F. In the estimation apparatus, input means for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data which is data representing the base sequence, and the input means Read the inputted DNA microarray disruption strain data, and select the gene group B consisting of the predetermined number of the other genes b having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray disruption strain data Selected by the gene group selection means, the data of the gene a, and the gene group selection means When the control region of the gene a is extracted as a character string based on the DNA base sequence data input by the input means for the pair of the data of each gene b of the gene group B, the i th from the head When the position of the character is defined as position i, a continuous region of a predetermined length l is cut out from each position i in the control region as _a window W _a [i], and the control region of each gene b of the gene group B is a character string. When the position of the j-th character from the beginning is taken as position j, a continuous area of a predetermined length l is cut out as window W _b [j] from each position j in the control area , and these windows W _a With respect to the character string data included in [i], W _b [j] , a character having a predetermined length in a range shifted by a predetermined number at each position in these windows . Calculate all the similarities between columns , calculate the window similarity that is the sum of the similarities for each position in the window when the calculated maximum similarity is defined as the similarity of that position , and calculate Based on window similarity data composed of the window similarity of each gene b of the gene group B, the maximum value Max _aB [i] of the similarity of the gene group B to each position i in the control region of the gene a and Window similarity calculating means for calculating an average value Avg _aB [i], data of each gene b of the gene group B selected by the gene group selecting means, and the window similarity calculating means calculated by the window similarity calculating means It reads the average value data consisting of a maximum value data and the average value _Avg aB consisting maximum _{Max aB [i] [i]} , the gene a Of your area portions to determine statistically specifically large first area in the distribution of the window similarity the maximum value Max aB _[i] is in the position i, calculated by the window similarity calculation means Based on the window similarity data of each gene b of the gene group B, the maximum value Max _ab [i] of the window similarity of each gene b of the gene group B with respect to the change in the position i of the control region of the gene a. ], The second region that includes the peak position when the position that gives the maximum value is determined as the peak position, and the region R that is the first region and the second region is searched for. Data of each gene b of the gene group B selected by the region search means, the gene group selection means, and the window similarity calculation means. In addition, the window similarity data of each gene b of the gene group B and the area data composed of the area R obtained by the area searching means are read, and among the genes b of the gene group B, each area R for, there is a peak located in the area R, and the DNA similarity entered and genetic group extracting means for extracting seeking genes B _R consisting of high gene b than the predetermined value, by said input means The gene extracted by the gene group extraction means for each region R using a change ratio of expression intensity between a wild strain not subjected to the experimental operation based on microarray disrupted strain data and a disrupted strain subjected to the experimental operation among the genes of group B _R, the second gene to extract only those expression intensity greatly changes in the experimental procedure as genes B _R ^F It is characterized in that it comprises an extraction unit.

In such a gene interaction estimation method , a gene interaction estimation program, and a gene interaction estimation apparatus according to the present invention, by introducing the concept of window similarity, it is not the entire control region of genes a and b that is localized. It is possible to find a similar degree of similarity.

Here, the gene interaction estimation method according to the present invention, the gene interaction estimation program, and in gene interaction estimation apparatus when calculating the window similarity statistically specifically frequency is higher string for, you a similarity to it with 0. That is, the applicant of the present application analyzed the frequency of occurrence of oligonucleotides in the control region, and newly discovered that oligonucleotides with a statistically high frequency of appearance are less likely to be binding strings. In the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention, it is possible to find a similar region other than the binding site by reflecting this fact in the region similarity calculation. Can be reduced.

Moreover, the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention are based on window similarity data composed of the window similarity of each gene b of the gene group B. The maximum value Max _aB [i] and the average value Avg _aB [i] of the similarity degree of the gene group B for each position i in the control region of the gene a are calculated. A gene interaction estimation method, a gene interaction estimation program, and a gene interaction estimation apparatus according to the present invention include a maximum value consisting of data of each gene b of the gene group B and the maximum value Max _aB [i]. reads the average value data consisting of data and the average value _Avg aB [i], of the control region of the gene a, statistically the maximum value _Max aB [i] is in the distribution of the window similarity at the position i together determine the specific size has a first region, on the basis of the window similarity data for each gene b of the genes B, each of the genes B with respect to a change in position i of the control region of the gene a the peak position when the position giving the maximum value and the peak position when viewed the maximum value Max ab window similarity _[i] of the gene b It obtains a second region that includes a, is characterized by searching for a region R is the a first region and the second region. That is, in the gene interaction estimation method , gene interaction estimation program, and gene interaction estimation apparatus according to the present invention, when a common binding site exists in genes a and b, a window from this peak position is displayed. Therefore, it is possible to narrow down the binding site candidates by searching for a region R that includes such a peak position.

Furthermore, the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention include the data of each gene b in the gene group B and the window similarity of each gene b in the gene group B. Degree data and region data composed of the region R, and for each region R among the genes b of the gene group B , a peak position exists in the region R, and the similarity is lower than a predetermined value. It is characterized by extracting seeking genes B _R consisting of high gene b. Thus, a gene interaction estimation method according to the present invention, the gene interaction estimation program, and in gene interaction estimation apparatus of genes B, genes B _R having a predetermined value or more similarity to the control region Can be extracted.

Thus, in the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention, these processes are performed as the second filtering process, without assuming a specific pattern, A gene having a control region similar to the designated gene a is exhaustively searched. Thus, it can be estimated that the area | region where the similarity degree calculated | required is a site which a common factor couple | bonds.

Furthermore, the gene interaction estimation method , gene interaction estimation program, and gene interaction estimation apparatus according to the present invention can be obtained by performing gene disruption as a specific experimental operation prior to calculating the window similarity. Reading DNA microarray-disrupted strain data, and selecting a gene group B composed of a predetermined number of the other genes b having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray-disrupted strain data It is characterized by. That is, in the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention, the correlation coefficient of the DNA microarray disrupted strain data is used only for filtering genes that are subject to common control. As described above, in the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention, the control region analysis process having a high calculation cost, that is, the second filtering process is performed before the second filtering process. By performing the process using the correlation of the expression data as the first filtering process, the calculation time of the second filtering process having a high calculation cost can be significantly reduced.

Further, the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention include a wild strain that does not add the experimental operation based on the DNA microarray disrupted strain data, and a disruption that adds the experimental operation. using variation ratio of the expression intensities of the line extract, for each region R, among the genes of the genes B _R, only those expression intensity greatly changes in the experimental procedure as genes B _R ^F It is characterized by doing. Specifically, a gene interaction estimation method according to the present invention, the gene interaction estimation program, and genetic interactions estimation apparatus factor group of strongly influenced the DNA microarray disrupted strain data to be had based Dzu the gene a the calculated, for each factor F was determined, based on a change ratio of the expression intensity in the disrupted strain of the factor F, among the genes of the genes B _R, a gene consisting of strong received genes influence of said factor F you extract the group as the group of genes _B ^{R F.} Thus, a gene interaction estimation method according to the present invention, the gene interaction estimation program, and in gene interaction estimation apparatus as a third filtering process, the genes B _R extracted by the second filtering process Among them, by selecting only those whose expression intensity is greatly changed in a specific experimental operation, a factor generated by the experimental operation or an unknown factor that is strongly influenced by the factor can be estimated as a transcription factor.

Furthermore, the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention include the data of each gene in the gene group B _R ^F for each region R and each factor F, and the base sequence Character string data included in a window from the peak position in each region R, for each gene of the gene a and the gene group B _R ^F , based on DNA base sequence data that is data representing It is characterized by performing multiple alignment for the peak position . Thereby, in the gene interaction estimation method , the gene interaction estimation program, and the gene interaction estimation apparatus according to the present invention, when the existence of the common pattern can be confirmed, it can be determined that the estimation is successful. It becomes possible.

Furthermore, the binding site estimation method according to the present invention that achieves the above-described object is characterized in that the binding site estimation method estimates a binding site to which a common factor that controls gene a and any other gene b binds. An input process for inputting DNA microarray disruption strain data obtained by performing gene disruption as an experimental operation and DNA base sequence data, which is data representing the base sequence, to a computer, and a processor in the computer is designated by A pair of the data of the gene a and the data of each gene b of the gene group B selected based on the DNA microarray disruption strain data input in the input step was input in the input step Based on the DNA base sequence data, the control region of the gene a is extracted as a character string. The continuous region of a predetermined length l from each position i in the control region when from the head to the location of the i-th character is position i with cut out as the window W _{a [i]} when, each gene of the genes B When the control area of b is extracted as a character string, the position of the jth character from the beginning is taken as position j, and a continuous area of a predetermined length l from each position j in the control area is defined as a window W _b [j]. For the character string data cut out and included in these windows W _a [i], W _b [j] , the similarity between character strings of a predetermined length in a range that is shifted by a predetermined number at each position in these windows the calculated all the maximum value of the calculated similarity to calculate the window similarity is the sum of similarity to each position in the window when defining the degree of similarity of its position, calculates Based on the window similarity data composed of the window similarity of each gene b of the gene group B, the maximum value Max _aB [i] of the gene group B for each position i in the control region of the gene a Window similarity calculation step for calculating the average value Avg _aB [i], the processor calculates the data of each gene b of the selected gene group B, and the maximum value calculated in the window similarity calculation step The maximum value data consisting of Max _aB [i] and the average value data consisting of the average value Avg _aB [i] are read, and the maximum value Max _aB [i] in the control region of the gene a is at the position i A first region that is statistically specifically large in the window similarity distribution is obtained and calculated in the window similarity calculation step. Based on the window similarity data of each gene b of the gene group B, the maximum value Max _ab [ max _{ab of the} window similarity of each gene b of the gene group B with respect to the change in the position i of the control region of the gene a When the position that gives the maximum value when i] is taken as the peak position, a second area including the peak position is obtained, and the area R that is the first area and the second area is searched. And an output step in which the processor outputs information indicating a binding site to which a factor common to the genes a and b binds based on the region data including the region R obtained in the region searching step. It is characterized by comprising.

The binding site estimation program according to the present invention that achieves the above-described object is a computer-executable binding site estimation program that estimates a binding site to which a common factor that controls gene a and any other gene b binds. In the above, an input means for inputting DNA microarray disrupted strain data obtained by performing gene disruption as a specific experimental operation and DNA base sequence data that is data representing the base sequence is designated. The DNA base sequence inputted by the input means for a pair of the data of the gene a and the data of each gene b of the gene group B selected based on the DNA microarray disrupted strain data inputted by the input means Based on the data, the control region of gene a is extracted as a character string The continuous region of a predetermined length l from each position i in the control region when from the head to the location of the i-th character is position i with cut out as the window W _{a [i]} when a respective of said gene group B When the control region of gene b is taken out as a character string, the position of the jth character from the beginning is taken as position j, and a continuous region of a predetermined length l from each position j in the control region is represented by window W _b [j] As for character string data included in these windows W _a [i] and W _b [j] , similarities between character strings of a predetermined length within a range shifted by a predetermined number at each position in these windows All the degrees are calculated, and the window similarity that is the sum of the similarities for each position in the window when the calculated maximum similarity is defined as the similarity of the position is calculated. Based on the window similarity data composed of the window similarity of each gene b of the gene group B, the maximum value Max _{aB of the} gene group B for each position i in the control region of the gene a Max _iB [i ] And the average value Avg _aB [i], the data of each gene b in the selected gene group B, and the maximum value Max _aB [i calculated by the window similarity calculation unit ] And the average value data consisting of the average value Avg _aB [i], and among the control region of the gene a, the maximum value Max _aB [i] is the window similarity at the position i. A first region that is statistically specifically large in the distribution of, and is calculated by the window similarity calculation unit Based on the window similarity data of each gene b of the gene group B, the maximum value Max _ab [i of the window similarity of each gene b of the gene group B with respect to the change of the position i of the control region of the gene a. ], The second region that includes the peak position when the position that gives the maximum value is determined as the peak position, and the region R that is the first region and the second region is searched for. Functioning as region search means and output means for outputting information indicating a binding site to which a factor common to each gene a and b binds based on region data consisting of the region R obtained by the region search means. It is characterized by.

Furthermore, a binding site estimation apparatus according to the present invention that achieves the above-described object is a binding site estimation apparatus that estimates a binding site to which a common factor that controls gene a and any other gene b binds. DNA microarray disruption strain data obtained by performing gene disruption as an experimental operation, input means for inputting DNA base sequence data that is data representing the base sequence, data of the designated gene a, Based on the DNA base sequence data input by the input unit, the gene for the pair of the data of each gene b of the gene group B selected based on the DNA microarray disrupted strain data input by the input unit When the control area of a is extracted as a character string, the position of the i-th character from the beginning is set as the position i. The continuous region of a predetermined length l from each position i in the control region with cut out as the window W _{a [i]} when the, from the beginning when the removed control region of each gene b of the genes B as a string When the position of the j-th character is a position j, a continuous area of a predetermined length l is cut out from each position j in the control area as a window W _b [j], and these windows W _a [i], W _b [j the character string data included in, the position the maximum value of the calculated all the calculated similarity similarity string between the predetermined length in the range of a predetermined number shifted back and forth at each position within the windows Window similarity, which is the sum of the similarities for each position in the window when defined as the similarity of the window, and the windows of each gene b of the gene group B calculated The maximum value Max _aB [i] and the average value Avg _aB [i] of the gene group B with respect to each position i in the control region of the gene a are calculated based on the window similarity data composed of the similarity. Maximum value data consisting of window similarity calculation means, data of each gene b of the selected gene group B, and the maximum value Max _aB [i] calculated by the window similarity calculation means, and the average value Avg _aB Average value data consisting of [i] is read, and the maximum value Max _aB [i] of the control region of the gene a is statistically specifically large in the distribution of the window similarity at the position i. And calculating the window similarity degree of each gene b of the gene group B calculated by the window similarity degree calculating means. Based on the data, the position giving the maximum value when the maximum value Max _ab [i] of the window similarity of each gene b of the gene group B with respect to the change of the position i of the control region of the gene a is peaked A second region that includes the peak position when it is set as a position, and a region search unit that searches for the region R that is the first region and the second region; and the region search unit And output means for outputting information indicating a binding site to which a factor common to the genes a and b binds based on the region data consisting of the region R.

In such a binding site estimation method, binding site estimation program, and binding site estimation apparatus according to the present invention, by introducing the concept of window similarity, local similarity rather than the entire control region of genes a and b is adopted. It becomes possible to discover. Then, the binding site estimation method, the binding site estimation program, and the binding site estimation apparatus according to the present invention are based on window similarity data including the window similarity of each gene b of the gene group B. Maximum value Max of similarity of gene group B for each position i in the control region _aB [I] and average value Avg _aB [I] is calculated, the data of each gene b of the gene group B, and the maximum value Max. _aB Maximum value data consisting of [i] and the average value Avg _aB Average value data consisting of [i] and reading the maximum value Max of the control region of the gene a _aB [I] obtains a first region that is statistically specifically large in the distribution of window similarity at the position i, and based on the window similarity data of each gene b of the gene group B, Maximum value Max of the window similarity of each gene b of the gene group B with respect to a change in the position i of the control region of a _ab When the position that gives the maximum value when viewing [i] is the peak position, a second region that includes the peak position is obtained, and the region R that is the first region and the second region is determined. Explore. Thereby, in the binding site estimation method, the binding site estimation program, and the binding site estimation apparatus according to the present invention, when there is a common binding site in the genes a and b, the window from this peak position is within the window. Since there is a high possibility of existence, it is possible to narrow down candidates for binding sites by searching for a region R including such a peak position.

According to the present invention, instead of examining how common each character string is included in the control area, a partial area (window) of a certain length is regarded as a collection of character strings, and the window as a whole By comprehensively calculating whether or not they are similar, it is possible to cope with expression control by a combination of character strings. According to the present invention, the influence of DNA microarray data and the influence of ambiguity included in the binding pattern can be reduced, and the gene interaction and the binding site are estimated with high accuracy while using the DNA microarray. be able to.

  Hereinafter, specific embodiments to which the present invention is applied will be described in detail with reference to the drawings.

  In this embodiment, when the gene a is controlled by the factor F (positive control or negative control), it is estimated how the gene interaction is estimated to affect the other gene b controlled by the factor F. It is. In particular, this gene interaction estimation method can estimate the gene interaction with high accuracy while performing DNA microarray by applying three independent filtering processes to the input data. is there.

  In the following description, actual data using gene expression data of Bacillus subtillis will be appropriately described as gene expression data obtained by performing an experimental operation on the organism to be estimated.

  A gene interaction estimation apparatus that executes this gene interaction estimation method is generally realized using a computer that executes a predetermined gene interaction estimation program by a processor. FIG. 1 shows a series of processes in a gene interaction estimation method to which the present invention is applied.

  First, in this gene interaction estimation method, DNA microarray disruption strain data, DNA base sequence data, and gene a data to be estimated for interaction are input to a gene interaction estimation device. The DNA microarray data is the expression data of the gene a obtained by the DNA microarray technology. The DNA microarray disrupted strain data is one type of DNA microarray data, and the expression level of the disrupted strain and the wild strain It is the data which observed ratio. Here, the disrupted strain is a strain in which the target gene a is destroyed by an experimental method, and the wild strain is a strain that has not destroyed the gene a. A strain is another individual having the same genetic character. Furthermore, DNA is a substance in which one base sequence composed of four base letters (a, t, g, c) and its complementary sequence are bonded to each other to form a double helix structure. The base sequence data is data representing the base sequence on one side of the sequence.

In this gene interaction estimation method, the first filtering process using the correlation of the expression data is performed on the input in order to reduce the influence of the background noise of the expression data by the DNA microarray. That is, in this gene interaction estimation method, in step S1, the input DNA microarray disruption strain data is read by the gene group selection means in the gene interaction estimation device , and the designated gene a and DNA microarray disruption strain data are read. A gene group B consisting of a predetermined number of other genes b having a large absolute value of the correlation coefficient is selected. Specifically, in this gene interaction estimation method, the top 300 gene groups B having a large absolute value of the correlation coefficient between the gene a and the DNA microarray disrupted strain data are selected. For example, if the gene a is "ahpC", as shown in FIG. 2, in order large absolute value of the correlation coefficient of the DNA microarray disrupted strain data, "ahpC", "ahpF" , "yurV", " The top 300 genes such as “yurW”, “yurU”,... are selected and set as gene group B.

That is, in this gene interaction estimation method, the correlation coefficient of DNA microarray disruption strain data is used only for filtering genes that are subject to common control. In other words, in this gene interaction estimation method, in order to narrow down the candidate genes that are subject to common control, the correlation coefficient value of the DNA microarray disruption strain data itself is not limited to a certain value or more. Processing is performed so that all those having a correlation coefficient are left as candidate genes that are subject to common control.

  Subsequently, in this gene interaction estimation method, a second filtering process is performed on the gene group B extracted by the first filtering process, and a binding site (binding site) by searching for a similar partial region in the control region. ).

First, in this gene interaction estimation method, as part of the second filtering process, in step S2 in FIG. 1, each of the data of gene a and each of gene group B is calculated by the window similarity calculation means in the gene interaction estimation device. For the pair with the gene b data , the window similarity of the control region is calculated based on the DNA base sequence data .

  Here, for example, as shown in FIG. 3, the window similarity is a concept newly devised in order to discover local similarity while allowing a slight shift of similar parts, not the entire control region of a gene. It is. At this time, in this gene interaction estimation method, a character string having an appearance frequency that is specifically higher than the expected value does not become a binding site, and the character string is commonly included in the control region of two genes to be compared. Based on the applicant's discovery that common factors do not bind even if they are included, only the similarity of strings that are likely to become binding sites is considered, and the characters are not related to gene regulation. Ignore if the columns are similar.

Specifically, for example, as shown in FIG. 4, the window similarity is obtained by cutting out continuous regions of length l (for example, about 30) from the control regions of two genes a and b to be compared as windows. The maximum value of the similarities between character strings (oligonucleotides) of length k (for example, about 6) having close start positions included in is defined as the similarity of the positions. Is calculated as The similarity between oligonucleotides is defined as a monotonically increasing function related to the number of matching characters. At this time, all the similarities in a range shifted by several (for example, about 3) before and after are calculated, and the maximum value is defined as the similarity of the position. As a result, this gene interaction estimation method approximates the motif alignment, and can correspond to a pattern including an arbitrary character string with a variable length in between, which is frequently seen in a combined character string. It becomes possible. Further, for oligonucleotides that are statistically specific and have a high appearance frequency, based on the discovery that they do not become binding sites, the similarity to them is set to “0” and excluded from the calculation target. Actually, for example, as shown in FIG. 4, if the window starting from the i-th character in the control region of gene a is W _a [i], then the window W _{b of} gene b for this window W _a [i]. The similarity wsim _{a, b} (i, j) of [j] is assumed to be “6” as the length of the character string.
wsim _{a, b} (i, j): = Σ _{k = 0, 1−6} sim _{a, b} (i + k, j + k)
Is calculated as

Here, the calculation amount of the similarity for one window pair is considered to be proportional to the window length l, and the total number of window pairs is n, where n is the length of the entire control region of the genes a and b. at most ^{two n.} Therefore, the amount of calculation for all window pairs should be O (n ² · l) because there are n ² pairs of windows with the amount of calculation being l. In contrast, in this gene interaction estimation method, sim _{a, b} (i, j) for all positions i, j is calculated in advance, where n is the length of the entire control region of genes a, b. The amount O (n ² ) can be calculated, and the similarity for all window pairs is wsim _{a, b} (i + 1, j + 1) = wsim _{a, b} (i, j) −sim _{a, b} ( If the relationship of i, j) + sim _{a, b} (i + 1-5, j + 1-5) is used, the calculation amount O (n ² ) can be used. Therefore, in the gene interaction estimation method, the similarity to all window pairs can be calculated with the calculation amount O (n ² ) without depending on the window length l, and the calculation amount is very small. You can do it.

Thus, in this gene interaction estimation method, the local similarity is found instead of the entire control region of genes a and b by incorporating the concept of window similarity. Further, in this gene interaction estimation method, when calculating the window similarity, each position i in the control region of the gene a is based on the window similarity data including the window similarity of each gene b of the gene group B. The maximum value Max _aB [i] and the average value Avg _aB [i] of the gene group B and the maximum value Max _ab [i] of the similarity of each gene b of the gene group B are calculated. Note that the position i means the position of the i-th character from the beginning when the control area is extracted as a character string. Further, the window of length l from the position i means partial character strings from the i-th to the (i + 1) -1th.

Subsequently, in this gene interaction estimation method, as part of the second filtering process, in step S3 in FIG. 1, the region search means in the gene interaction estimation device uses the data of each gene b of the gene group B, The maximum value data consisting of the maximum value Max _aB [i] and the average value data consisting of the average value Avg _aB [i] are read, and the maximum value Max _aB [i] in the control region of the gene a is similar to the window at the position i. A region R that is statistically specifically large in the degree distribution and includes the peak position is searched. More specifically, in this gene interaction estimation method, the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i] is calculated, and the gene a Among the control regions, a region R in which the calculated difference Max _aB [i] −Avg _aB [i] is larger than a predetermined value and includes the peak position is searched. Note that the peak position refers to the position of the control region of the gene a by shifting the position of the gene a back and forth based on the window similarity data of each gene b of the gene group B , for example, as indicated by an arrow in FIG. This is the position that gives the maximum value when _{looking at} the maximum value Max _ab [i] of the window similarity of each gene b of the gene group B with respect to the change of i. At this time, in this gene interaction estimation method, only the peak position having a similarity equal to or higher than a predetermined threshold is targeted, and therefore the difference Max between the maximum value Max _aB [i] and the average value Avg _aB [i]. _{aB [i] -Avg aB [i} ] is directed to a peak position included in the larger continuous area than a predetermined value. And in this gene interaction estimation method, when there is a common binding site for genes a and b, it is highly possible that they exist within the window from this peak position, so this peak position is included. Such a region R is searched.

Specifically, in this gene interaction estimation method, for example, as shown in FIG. 6, for example, the maximum value Max _aB [i] and the average value Avg _aB [i] with respect to a change in the position i of the control region of the gene a ] And a peak position indicated by a solid plot in the same figure is obtained. In this gene interaction estimation method, the peak position having a similarity of 18000 or more is targeted, and the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i]. Assuming that an area where the value of i] continuously has a value of 15000 or more and includes a peak position is searched, for example, as shown by the hatched portion in FIG. Asking.

In order to search the region R, the difference Max _aB [i] −Avg _aB [i] is used instead of the maximum value Max _aB [i] for the following reason. As described above, since an oligonucleotide with a statistically specific high appearance frequency is set to “0” and excluded from the calculation target, it is used for calculating the similarity for each window. The total number of strings will be different. Therefore, there is a possibility that the same degree of window similarity is calculated when there are only a few character strings with high similarity in the window and when there are many character strings with low similarity in the window. . In the gene interaction estimation method, it is necessary to select only the former because the purpose is to select one having a high degree of similarity. Here, in the former and the latter, since the latter has a higher average value, the former and the latter can be simplified by looking at the difference between the maximum value Max _aB [i] and the average value Avg _aB [i]. Can be distinguished. This is the reason why the difference Max _aB [i] −Avg _aB [i] is used instead of the maximum value Max _aB [i] to search the region R. However, in the gene interaction estimation method, if the maximum value Max _aB [i] is determined to be statistically specifically large in the distribution of window similarity at the position i, the difference Max _aB [I] -Avg _{aB Even a} method other than the method using [i] can be applied.

In gene interaction estimation method, this way obtaining the region R, as a part of the second filtering process, in step S4 in FIG. 1, the gene group extracting means in gene interaction estimation apparatus, the genes B The data of each gene b, the window similarity data of each gene b of gene group B, and the area data consisting of the obtained area R are read, and for each area R of each gene b of gene group B, the corresponding area there is a peak located in R, and the similarity to seek genes B _R consisting of high gene b than a predetermined value. That is, in this gene interaction estimation method, by performing this step S4, among the genes B extracted by the first filtering process, genes B _R having a predetermined value or more similarity to the control region Can be extracted.

  Thus, in the gene interaction estimation method, as the second filtering process, genes having a control region similar to the designated gene a are exhaustively searched without assuming a specific pattern. At this time, in this gene interaction estimation method, when calculating the similarity, a portion that is unlikely to be included in the binding site is excluded based on the statistical specificity regarding the appearance frequency of the partial character string. Thus, it can be estimated that the area | region where the similarity degree calculated | required is a site which a common factor couple | bonds.

In the gene interaction estimation method, performed in step S4 in FIG. 1, when obtaining the genes B _R, the process of estimating the transcription factor by using the gene expression data obtained by adding experimental procedure to be estimated biological. That is, in this gene interaction estimation method, the direct influence of the experimental operation can be known by comparing the change in gene expression level between the individual to which the experimental operation is not applied and the individual to which the experimental operation is added. Specifically, in the gene interaction estimation method, the gene extracted by the second filtering process for each region R as the third filtering process by the second gene group extracting means in the gene interaction estimating apparatus. From the group B _R, it selects only those large expression intensities in the experimental operation is changing. Here, the factor group which strongly influenced the gene a in the DNA microarray disruption strain data obtained by performing gene disruption as an experimental operation is obtained. Then, in this gene interaction estimation method, for each factor F obtained, among the genes B _R, extracts strongly received genes influence factors F, obtains it as B _R ^F. Incidentally, genes B _R ^F has a large change in the expression level of disrupted strain compared with the expression intensities of the wild strain by factor F, and, those that order is also higher. The change ratio of the expression intensity is represented by the absolute value of log ((expression intensity of the disrupted strain) / (expression intensity of the wild strain)).

For example, it is assumed that the genes B _R as shown in FIG. 8 are extracted. In this case, in the gene interaction estimation method, as shown in FIG. 9, the factor generated by the disrupted gene with respect to the disrupted strain data in which the expression intensity of “ahpC” as gene a has changed has a strong influence on the gene a. As a factor group. In this gene interaction estimation method, for example, as shown in FIG. 10, the change ratio of the expression intensity in the disrupted strain of the gene “perR” that produces “PerR” as the obtained factor F is obtained. In the figure, the numbers in parentheses indicate the rank order of the change in the expression intensity within the gene. Thus, in this gene interaction estimation method, among the genes B _R, influence strongly received genes of "ahpC" one of the strongly influential factor F is in "PerR" as a gene a Extracting “katA”, “yfmJ”, “mrgA”, “hemA”, “ponA”, “ykvW”, and “ydjL” as B _R ^F as shown in bold in FIGS. it can.

Thus, in the gene interaction estimation method, as a third filtering process, among the genes B _R extracted by the second filtering process, only those that expression intensity greatly changes in certain experimental manipulation By selecting, a factor generated by an experimental operation or an unknown factor strongly influenced by the factor can be estimated as a transcription factor.

In this gene interaction estimation method, in step S5 in FIG. 1, the data of each gene of the gene group B _R ^F for each region R and each factor F is input and input by the alignment means in the gene interaction estimation device. Based on the obtained DNA base sequence data, for each gene of gene a and gene group B _R ^F , character string data included in the window from the matched position (peak position in each region R) is extracted, for example, FIG. As shown in FIG. 11, multiple alignment is performed with respect to the peak position . In this gene interaction estimation method, if the existence of the common pattern can be confirmed, it is determined that the estimation is successful. Thus, in this gene interaction estimation method, for example, as shown in FIG. 12, genetic elements to control a F, control direction (positive or negative), genes other than the gene a controlled by factors F B _R Information indicating the binding site of each gene of ^F , gene a, and gene group B _R ^F is output, and the series of processing ends.

Thus, in the gene interaction estimation method, a plurality of results obtained by performing three filtering processes are used in combination, and gene a receives control of factor F (positive control or negative control). At this time, it is estimated what kind of influence the other gene b under the control of factor F has. At this time, in this gene interaction estimation method, data with low reliability such as correlation coefficient of DNA microarray disruption strain data is used only for rough selection, and data with high reliability such as experimental operation data is used. By using the value itself as a selection criterion, the reliability of the result can be improved.

  In the above description, specific results are shown for the case where “ahpC” as gene a is controlled by “PerR” as factor F, using FIG. 2 and FIGS. Hereinafter, other experimental results obtained by concretely executing such a method for estimating gene interaction are also shown for reference.

  First, specific experimental results when “purA” is used as gene a will be described.

First, in the experiment, similar to the processing in step S1 in FIG. 1, the top 300 gene group B having the large absolute value of the correlation coefficient between the gene a and the DNA microarray disrupted strain data was selected. Thus, for example, as shown in FIG. 13, the genes “purA”, “yumD”, “glyA”, “ykbA”, “purQ”,... In descending order of the absolute value of the correlation coefficient of the DNA microarray disrupted strain data. -Was selected as gene group B.

Subsequently, in the experiment, similarly to the processing in step S2 in FIG. 1, the window similarity of the control region is calculated for the pair of the data of gene a and the data of each gene b of gene group B, and in step S3 Similar to the processing, in the control region of gene a, the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i] is larger than a predetermined value, and The region R including the peak position was searched. Here, the peak position with a similarity of 18000 or more is targeted, and the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i] is 15000 or more. Since the condition was set so as to search for a region having values continuously and including the peak position, positions i = 4156019 to 4156029 were obtained as the region R as indicated by the hatched portion in FIG.

Subsequently, in the experiment, as in the processing in the step S4 1, there is a peak located in the determined area R, and the degree of similarity is extracted genes B _R consisting of high gene b than the predetermined value . Thus, genes B _R as shown in FIG. 15 is extracted. Then, in the experiment, determine the factor group gave strongly influence gene a in DNA microarray disrupted strain data, for each factor F obtained, among the genes of the genes B _R, was strongly influenced by factors F extract the genes consisting of the genes was determined it as B _R ^F. Thus, as shown in FIG. 16, broken lines expressing the strength of "purA" as genes a has changed, obtained as factor group of strongly influenced in the gene "purA", each of a group of genes B _R Among the genes , as a gene group B _R ^F strongly influenced by “PurR”, which is one of the factors F that strongly influenced the gene purA ”, as shown in bold letters in FIG. 15 and FIG. ydeQ "and" glyA "were extracted.

Then, in the experiment, as in the processing in the step S5 1, for each gene in the gene a and gene group B _R ^F, takes out character string data included in the window from the matched positions, performs multiple alignment, Results as shown in FIG. 18 were obtained.

Thereby, the factor F controlling the gene a, the control direction (positive or negative), the gene group B _R ^F other than the gene a controlled by the factor ^F , and the binding site of each gene of the gene a and the gene group B _R ^F As information indicating the above, a result as shown in FIG. 19 was obtained as an output.

  Next, specific experimental results when “phoD” is used as gene a will be described.

First, in the experiment, similar to the processing in step S1 in FIG. 1, the top 300 gene group B having the large absolute value of the correlation coefficient between the gene a and the DNA microarray disrupted strain data was selected. Accordingly, for example, as shown in FIG. 20, the genes “phoD”, “yqgG”, “tuaD”, “tuaB”, “yqgK”,... In descending order of the absolute value of the correlation coefficient of the DNA microarray disrupted strain data. -Was selected as gene group B.

Subsequently, in the experiment, similarly to the processing in step S2 in FIG. 1, the window similarity of the control region is calculated for the pair of the data of gene a and the data of each gene b of gene group B, and in step S3 Similar to the processing, in the control region of gene a, the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i] is larger than a predetermined value, and The region R including the peak position was searched. Here, the peak position with a similarity of 18000 or more is targeted, and the difference Max _aB [i] −Avg _aB [i] between the maximum value Max _aB [i] and the average value Avg _aB [i] is 15000 or more. Since the condition was set so as to search for an area having values continuously and including the peak position, positions i = 283493 to 283502 were obtained as the area R as indicated by the hatched portion in FIG.

Subsequently, in the experiment, as in the processing in the step S4 1, there is a peak located in the determined area R, and the degree of similarity is extracted genes B _R consisting of high gene b than the predetermined value . Thus, genes B _R as shown in FIG. 22 is extracted. Then, in the experiment, determine the factor group gave strongly influence gene a in DNA microarray disrupted strain data, for each factor F obtained, among the genes of the genes B _R, was strongly influenced by factors F extract the genes consisting of the genes was determined it as B _R ^F. Thus, as shown in FIG. 23, broken lines expressing the strength of "phoD" as genes a has changed, obtained as factor group of strongly influenced in the gene "phoD", each of a group of genes B _R of the genes, as strongly received genes _B ^{R F} the effect of the gene "phoD" one of the strongly influential factor F is in "PhoP", as shown in bold in FIGS. 22 and 24, “PhoA” and “phoB” were extracted.

Then, in the experiment, as in the processing in the step S5 1, for each gene in the gene a and gene group B _R ^F, takes out character string data included in the window from the matched positions, performs multiple alignment, Results as shown in FIG. 25 were obtained.

Thereby, the factor F controlling the gene a, the control direction (positive or negative), the gene group B _R ^F other than the gene a controlled by the factor ^F , and the binding site of each gene of the gene a and the gene group B _R ^F As information indicating the above, a result as shown in FIG. 26 was obtained as an output.

  As described above, in the gene interaction estimation method shown as the embodiment of the present invention, the first filtering process using the correlation of the expression data, the second method that is a completely new technique for estimating the binding site. The influence of DNA microarray data and the ambiguity included in the binding pattern are reduced by fusing three independent estimation methods based on the filtering process and the third filtering process using gene expression data with added experimental operations. It is possible to estimate a gene interaction with high accuracy while using a DNA microarray.

  In particular, in this gene interaction estimation method, by performing processing that uses the correlation of expression data as the first filtering process prior to the analysis process of the control region with high calculation cost, that is, the second filtering process, It is possible to greatly reduce the calculation time of the second filtering process with high calculation cost. The applicant of the present application uses an AMD CPU “Opteron” (registered trademark) having a speed about five times faster than Intel CPU “Pentium (registered trademark) 4, clock 3 GHz”. (Trademark) When 2.4 GHz "is executed using a computer having a 4-parallel configuration, it is confirmed that it takes about 10 days to calculate all 4000 genes. On the other hand, the applicant of this application has also confirmed that the time of 10 times or more of this number of days is required when the first filtering process is not performed.

  In addition, in this gene interaction estimation method, as the second filtering process, as is frequently performed in the discovery in the promoter region, it is not possible to examine how common each oligonucleotide is included in the control region. Rather, a partial region (window) of a certain length is regarded as an aggregate of oligonucleotides, and it is comprehensively calculated whether or not the windows are similar as a whole. Thereby, in this gene interaction estimation method, it is possible to cope with expression control by a combination of oligonucleotides. Furthermore, the applicant of the present application has analyzed the appearance frequency of oligonucleotides in the control region, and newly discovered that an oligonucleotide having a statistically high frequency of occurrence is less likely to be a binding character string. In this gene interaction estimation method, by reflecting this fact on the calculation of the similarity between regions, the possibility that a similar region is found outside the binding site can be reduced. Furthermore, in this gene interaction estimation method, since the binding site can be expected to exist at the peak position, it is also possible to narrow down the binding site candidates.

  In this gene interaction estimation method, the third filtering process depends on only single DNA microarray data and is susceptible to noise. Therefore, the third filtering process is applied to the final stage. By performing the filtering process, it is possible to give a sufficiently large change in the expression intensity with respect to the background noise of the expression data by the DNA microarray, and it is possible to derive the direct influence of the factor on each gene. Furthermore, in this gene interaction estimation method, when performing the third filtering process, it is possible to create a situation equivalent to performing a relative evaluation by considering not only the value but also the rank.

  Note that the present invention is not limited to the above-described embodiment, and it is needless to say that modifications can be made as appropriate without departing from the spirit of the present invention.

It is a flowchart of a series of processes in the gene interaction estimation method shown as an embodiment of the present invention. It is a figure explaining the specific example of the gene group B of the upper predetermined number with a large absolute value of the correlation coefficient of the gene a and DNA microarray destruction strain data . It is a figure for demonstrating window similarity. It is a figure for demonstrating the specific calculation method of a window similarity. It is a figure for demonstrating a peak position. It is a figure for demonstrating a mode that a peak position is calculated | required concretely. It is a figure for demonstrating a mode that the area | region R is calculated | required concretely. Specific examples of the extracted genes B _R is a diagram illustrating a. It is a figure explaining the specific example of the destruction strain from which the expression intensity of "ahpC" as gene a changed. FIG. 4 is a diagram for explaining a specific example of a change ratio of expression intensity in a disrupted strain of gene “perR” as factor F. It is a figure for demonstrating the specific method of multiple alignment. It is a figure explaining the specific example of the information output through a series of processes in the gene interaction estimation method shown as embodiment of this invention. It is a figure explaining the specific example of the gene group B of the upper predetermined number with a large absolute value of the correlation coefficient of "purA" as the gene a and DNA microarray destruction strain data . It is a figure for demonstrating the area | region R calculated | required based on "purA" as the gene a and the gene group B shown in FIG. Specific examples of genes B _R extracted from the region R shown in FIG. 14 is a diagram for explaining the. It is a figure explaining the specific example of the destruction strain in which the expression intensity of "purA" as gene a changed. FIG. 6 is a diagram for explaining a specific example of a change ratio of expression intensity in a disrupted strain of gene “purR” as factor F. For each gene of "purA" and genes B _R ^F as a gene a, takes out character string data included in the window from the matched positions is a diagram for explaining results of multiple alignment. It is a figure explaining the specific example of the information as an experimental result output based on the data of FIG. 13 thru | or FIG. It is a figure explaining the specific example of the gene group B of the upper predetermined number with a large absolute value of the correlation coefficient of "phoD" as the gene a and DNA microarray destruction strain data . It is a figure for demonstrating the area | region R calculated | required based on "phoD" as the gene a and the gene group B shown in FIG. Specific examples of genes B _R extracted from the region R shown in FIG. 21 is a diagram illustrating a. It is a figure explaining the specific example of the destruction strain from which the expression intensity of "phoD" as gene a changed. FIG. 6 is a diagram for explaining a specific example of a change ratio of expression intensity in a disrupted strain of gene “phoP” as factor F. For each gene of "phoD" and genes B _R ^F as a gene a, takes out character string data included in the window from the matched positions is a diagram for explaining results of multiple alignment. It is a figure explaining the specific example of the information as an experimental result output based on the data of FIG. 20 thru | or FIG. It is a figure explaining the structure of a gene. It is a figure explaining the specific model of a Boolean network. It is a figure explaining the specific model of a Bayesian network.

Explanation of symbols

a, b gene Avg _aB [i] similarity average value _B, B _{R, B} ^{R F} genes F factor Max _aB [i] similarity maximum R region _{W a [i], W b} [j] Window wsim _{a, b} (i, j) Window similarity

Claims (9)

When the gene a is controlled by the factor F, in the gene interaction estimation method for estimating the influence on the other gene b controlled by the factor F,
An input step for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation and DNA base sequence data, which is data representing the base sequence, to a computer;
The processor in the computer reads the DNA microarray disrupted strain data input in the input step, and the upper predetermined number of the correlation coefficients between the designated gene a and the DNA microarray disrupted strain data is large. A gene group selection step of selecting a gene group B consisting of another gene b;
Wherein the processor, and data of said genes a, for a pair of each gene b data of said gene group B selected in the genes selection step, based on the DNA base sequence data input in said input step Then, when the control area of the gene a is taken out as a character string, a continuous area of a predetermined length l from each position i in the control area is defined as a window W _a when the position of the i-th character from the beginning is set as the position i [I], and when the control region of each gene b of the gene group B is extracted as a character string, the position of the j-th character from the beginning is defined as the position j, from each position j in the control region. cut out the continuous region of a predetermined length l as a window W _b [j], these windows W _{a [i],} the character string data contained in W b [j] Calculating all the similarity of the string between the predetermined length in the range of a predetermined number shifted back and forth at each position within these windows, the when the maximum value of the calculated degree of similarity is defined as the similarity of its position A window similarity that is the sum of the similarities for each position in the window is calculated, and the control region of the gene a is calculated based on the window similarity data that includes the window similarity of each gene b of the gene group B A window similarity calculating step of calculating a maximum value Max _aB [i] and an average value Avg _aB [i] of the similarity of the gene group B for each position i in FIG .
Maximum value data consisting of the data of each gene b of the gene group B selected in the gene group selection step and the maximum value Max _aB [i] calculated in the window similarity calculation step, and Average value data Avg _aB [i] is read, and among the control region of the gene a, the maximum value Max _aB [i] is statistically specific in the window similarity distribution at the position i. To the change of the position i of the control region of the gene a based on the window similarity data of each gene b of the gene group B calculated in the window similarity calculation step peak position a position which gives the maximum value when viewed the maximum value Max ab window similarity _[i] of each gene b of the genes B An area searching step the determined second area to include the peak position, to explore the region R is the a first region and said second region when a,
Data of each gene b of the gene group B selected in the gene group selection step and the window similarity data of each gene b of the gene group B calculated in the window similarity calculation step The region data consisting of the region R obtained in the region search step is read, and for each region R among the genes b of the gene group B, a peak position exists in the region R and the degree of similarity There the genes extraction step of extracting seeking genes B _R consisting of high gene b than the predetermined value,
The processor uses the change ratio of the expression intensity between the wild strain not subjected to the experimental operation based on the DNA microarray disrupted strain data input in the input step and the disrupted strain subjected to the experimental operation, to each region. for R, among the genes of the genes B _R extracted by said genes extraction step, a second gene to extract only those expression intensity greatly changes in the experimental procedure as genes B _R ^F A method for estimating a gene interaction, comprising: a group extraction step . In the area search step, the processor calculates a difference _{Max aB [i] -Avg aB [} i] and the maximum value _Max aB [i] and the average value _Avg aB [i], the calculated difference _{Max aB} If _{[i] -Avg} aB _[i] is greater than a predetermined value, the said maximum value _Max aB [i] is determined to statistically specifically greater in the distribution of the window similarity at the position i gene interaction estimation method of claim 1, wherein Rukoto determined as the first area. In the second group of genes extraction step, the processor determines a strong influencing factor group gave the gene a and have groups Dzu the DNA microarray fracture strain data input in said input step, each determined factor for F, the factors on the basis of the change ratio of the expression intensity in the disrupted strain F, among the genes of the genes B _R, said factor F influence strongly received a gene group the genes consisting of gene B _R of gene interaction estimation method according to claim 1 to wherein Rukoto extracted as ^F. Data of each gene of the gene group B _R ^F for each region R and each factor F extracted in the second gene group extraction step by the processor, and the DNA base sequence input in the input step Alignment that extracts character string data included in a window from the peak position in each region R and performs multiple alignment with respect to the peak position for each gene of the gene a and the gene group B _R ^F based on the data The method according to claim 3, further comprising a step. When a gene a is controlled by factor F, a computer-executable gene interaction estimation program for estimating an effect on another gene b controlled by the factor F,
The computer,
Input means for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data which is data representing the base sequence;
The DNA microarray-disrupted strain data input by the input means is read, and the gene consisting of the predetermined number of the other genes b having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray-disrupted strain data Gene group selection means for selecting group B,
And data of said genes a, for a pair of each gene b data of said gene group B selected by said genes selecting means, based on the DNA base sequence data inputted by said input means, said gene a When the position of the i-th character from the beginning is taken as position i when the control area is extracted as a character string, a continuous area of a predetermined length l is cut out as window W _a [i] from each position i in the control area. In addition, when the control region of each gene b of the gene group B is taken out as a character string, the position of the jth character from the beginning is defined as a position j, and a predetermined length l from each position j in the control region. cut out region as the window W _b [j], these windows W _{a [i],} the character string data contained in W b [j], these Said window when calculating all the similarity of the string between the predetermined length at a predetermined number shift range back and forth at each position in the window, the maximum value of the calculated degree of similarity is defined as the similarity of its position The window similarity that is the sum of the similarities for each position of the gene group B is calculated, and each window in the control region of the gene a is calculated based on the window similarity data including the window similarities of the genes b of the gene group B. Window similarity calculation means for calculating a maximum value Max _aB [i] and an average value Avg _aB [i] of the similarity of the gene group B to the position i ;
Maximum value data consisting of the data of each gene b of the gene group B selected by the gene group selection means, the maximum value Max _aB [i] calculated by the window similarity calculation means, and the average value Avg _aB Average value data consisting of [i] is read, and the maximum value Max _aB [i] of the control region of the gene a is statistically specifically large in the distribution of the window similarity at the position i. , And based on the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, the gene group B with respect to a change in the position i of the control region of the gene a peak position a position which gives the maximum value when the viewed the window maximum similarity value Max ab _[i] of each gene b of And to determine the second region to include the peak position when the said first region and is and the second area search means for searching for a region R is a region,
Data of each gene b of the gene group B selected by the gene group selection means, the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, and the region search Region data consisting of the region R obtained by the means, and for each region R of the genes b of the gene group B, there is a peak position in the region R, and the similarity is greater than a predetermined value. genes extracting means for extracting even seeking genes B _R consisting of high gene b and,
For each region R, using the ratio of expression intensity change between the wild strain not subjected to the experimental operation based on the DNA microarray disrupted strain data input by the input means and the disrupted strain subjected to the experimental operation, the gene among the genes of the genes B _R extracted by the group extraction unit, functions only what expression intensity greatly changes in the experimental procedure as a second gene group extracting means for extracting a group of genes B _R ^F A gene interaction estimation program characterized by When a gene a is controlled by factor F, in a gene interaction estimation apparatus for estimating an effect on another gene b controlled by the factor F,
An input means for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data which is data representing the base sequence;
The DNA microarray-disrupted strain data input by the input means is read, and the gene consisting of the predetermined number of the other genes b having a large absolute value of the correlation coefficient between the designated gene a and the DNA microarray-disrupted strain data Gene group selection means for selecting group B;
And data of said genes a, for a pair of each gene b data of said gene group B selected by said genes selecting means, based on the DNA base sequence data inputted by said input means, said gene a When the position of the i-th character from the beginning is taken as position i when the control area is extracted as a character string, a continuous area of a predetermined length l is cut out as window W _a [i] from each position i in the control area. In addition, when the control region of each gene b of the gene group B is taken out as a character string, the position of the jth character from the beginning is defined as a position j, and a predetermined length l from each position j in the control region. cut out region as the window W _b [j], these windows W _{a [i],} the character string data contained in W b [j], these Said window when calculating all the similarity of the string between the predetermined length at a predetermined number shift range back and forth at each position in the window, the maximum value of the calculated degree of similarity is defined as the similarity of its position The window similarity that is the sum of the similarities for each position of the gene group B is calculated, and each window in the control region of the gene a is calculated based on the window similarity data including the window similarities of the genes b of the gene group B. Window similarity calculation means for calculating the maximum value Max _aB [i] and the average value Avg _aB [i] of the similarity of the gene group B to the position i ;
Maximum value data consisting of the data of each gene b of the gene group B selected by the gene group selection means, the maximum value Max _aB [i] calculated by the window similarity calculation means, and the average value Avg _aB Average value data consisting of [i] is read, and the maximum value Max _aB [i] of the control region of the gene a is statistically specifically large in the distribution of the window similarity at the position i. , And based on the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, the gene group B with respect to a change in the position i of the control region of the gene a peak position a position which gives the maximum value when the viewed the window maximum similarity value Max ab _[i] of each gene b of An area searching means for searching a region R that determine the second region is the first region is a and the second region include the peak position when the,
Data of each gene b of the gene group B selected by the gene group selection means, the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, and the region search Region data consisting of the region R obtained by the means, and for each region R of the genes b of the gene group B, there is a peak position in the region R, and the similarity is greater than a predetermined value. and genes extracting means for extracting seeking genes B _R consisting of high gene b also,
For each region R, using the ratio of expression intensity change between the wild strain not subjected to the experimental operation based on the DNA microarray disrupted strain data input by the input means and the disrupted strain subjected to the experimental operation, the gene among the genes of the genes B _R extracted by the group extracting means, and a second gene group extracting means for extracting only those that expression intensity greatly changes in the experimental procedure as genes B _R ^F A gene interaction estimation device characterized by comprising: In a binding site estimation method for estimating a binding site to which a common factor that controls gene a and any other gene b binds,
An input step for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation and DNA base sequence data, which is data representing the base sequence, to a computer;
For a pair of the data of the gene a selected by the processor in the computer and the data of each gene b of the gene group B selected based on the DNA microarray disrupted strain data input in the input step, Based on the DNA base sequence data input in the input step, when the control region of the gene a is taken out as a character string, the location of the i-th character from the head is set as the position i. A continuous region of a predetermined length l is cut out from each position i as _a window W _a [i], and when the control region of each gene b of the gene group B is extracted as a character string, the location of the jth character from the beginning is determined. When a position j is set, a continuous area of a predetermined length l is cut out from each position j in the control area as a window W _b [j], For character string data included in these windows W _a [i], W _b [j] , all similarities between character strings of a predetermined length in a range shifted by a predetermined number at each position in these windows are obtained. Calculating a window similarity that is a sum of similarities for each position in the window when the calculated maximum similarity is defined as the similarity of the position, and calculating each gene of the gene group B The maximum value Max _aB [i] and the average value Avg _aB [i] of the similarity of the gene group B for each position i in the control region of the gene a based on the window similarity data consisting of the window similarity of b A window similarity calculation step for calculating
The processor _{sets the} maximum value data including the data of each gene b of the selected gene group B and the maximum value Max _aB [i] calculated in the window similarity calculation step, and the average value Avg _aB [i]. And a first region where the maximum value Max _aB [i] is statistically specifically large in the distribution of the window similarity at the position i among the control regions of the gene a. And obtaining each gene of the gene group B with respect to a change in the position i of the control region of the gene a based on the window similarity data of the gene b of the gene group B calculated in the window similarity calculation step the peak position when the position giving the maximum value and the peak position when viewed the maximum value Max ab window similarity _[i] and b An area searching step of searching the region R such seek second region is the first region is a and the second region to include,
The processor includes an output step of outputting information indicating a binding site to which a factor common to the genes a and b binds based on the region data including the region R obtained in the region searching step.
A method for estimating a binding site, characterized by: In a computer-executable binding site estimation program for estimating a binding site to which a common factor that controls gene a and any other gene b binds,
The computer,
Input means for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data which is data representing the base sequence;
For the pair of the specified gene a data and the data of each gene b of the gene group B selected based on the DNA microarray disruption strain data input by the input means, the data input by the input means Based on the DNA base sequence data, when the control region of the gene a is extracted as a character string, the position of the i-th character from the beginning is defined as the position i, and a predetermined length l from each position i in the control region A continuous region is cut out as _a window W _a [i], and when the control region of each gene b of the gene group B is extracted as a character string, the control is performed when the position of the j-th character from the beginning is set as the position j cut out the continuous region of a predetermined length l as a window _W b [j] from each position j in the area, these windows W _a [i] The character string data contained in W b _[j], the maximum of all the similarity string between the predetermined length calculated, the calculated similarity to the extent that a predetermined number shifted back and forth at each position within the windows A window comprising the window similarity of each gene b in the gene group B, which is calculated as a window similarity that is the sum of the similarities for each position in the window when the value is defined as the similarity of the position. Window similarity calculating means for calculating the maximum value Max _aB [i] and the average value Avg _aB [i] of the gene group B for each position i in the control region of the gene a based on the similarity data ,
The data of each gene b of the selected gene group B, the maximum value data consisting of the maximum value Max _aB [i] calculated by the window similarity calculation means, and the average value consisting of the average value Avg _aB [i] And a first region where the maximum value Max _aB [i] is statistically specifically large in the distribution of the window similarity at the position i among the control regions of the gene a , and Based on the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, the window of each gene b of the gene group B with respect to a change in the position i of the control region of the gene a including the peak position when the position giving the maximum value and the peak position when viewed maximum value Max ab similarity _[i] Obtains a second region as the first region a and and said second area search means for searching for a region R is a region and,
Functioning as output means for outputting information indicating a binding site to which a common factor binds to each of the genes a and b based on the region data consisting of the region R obtained by the region searching means.
A binding site estimation program characterized by In a binding site estimation apparatus for estimating a binding site to which a common factor that controls gene a and any other gene b binds,
An input means for inputting DNA microarray disruption strain data obtained by performing gene disruption as a specific experimental operation, and DNA base sequence data which is data representing the base sequence;
For the pair of the specified gene a data and the data of each gene b of the gene group B selected based on the DNA microarray disruption strain data input by the input means, the data input by the input means Based on the DNA base sequence data, when the control region of the gene a is extracted as a character string, the position of the i-th character from the beginning is defined as the position i, and a predetermined length l from each position i in the control region A continuous region is cut out as _a window W _a [i], and when the control region of each gene b of the gene group B is extracted as a character string, the control is performed when the position of the j-th character from the beginning is set as the position j cut out the continuous region of a predetermined length l as a window _W b [j] from each position j in the area, these windows W _a [i] The character string data contained in W b _[j], the maximum of all the similarity string between the predetermined length calculated, the calculated similarity to the extent that a predetermined number shifted back and forth at each position within the windows A window comprising the window similarity of each gene b in the gene group B, which is calculated as a window similarity that is the sum of the similarities for each position in the window when the value is defined as the similarity of the position. Window similarity calculating means for calculating the maximum value Max _aB [i] and the average value Avg _aB [i] of the gene group B for each position i in the control region of the gene a based on the similarity data When,
The data of each gene b of the selected gene group B, the maximum value data consisting of the maximum value Max _aB [i] calculated by the window similarity calculation means, and the average value consisting of the average value Avg _aB [i] And a first region where the maximum value Max _aB [i] is statistically specifically large in the distribution of the window similarity at the position i among the control regions of the gene a , and Based on the window similarity data of each gene b of the gene group B calculated by the window similarity calculation means, the window of each gene b of the gene group B with respect to a change in the position i of the control region of the gene a including the peak position when the position giving the maximum value and the peak position when viewed maximum value Max ab similarity _[i] An area searching means for searching a region R that determine the second region is the first region is a and the second region as,
Output means for outputting information indicating a binding site to which a common factor is bound to each of the genes a and b based on the area data composed of the area R obtained by the area searching means.
A combined site estimation device characterized by the above.

Images