JP2005149037A - Method, apparatus and program for estimating gene expression action - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To obtain a method for estimating interaction of gene expression in which an initial value of a suitable coefficient can be easily obtained and the volume of calculation is reduced. <P>SOLUTION: The method for preparing a gene model and estimating gene expression action on the basis of time-sequence data of a wild strain composed of a gene group which is not operated and destroyed strain obtained by destroying specific genes is provided with individual gene expression coefficient estimation steps (S17-S20) in which a loop for repeating a procedure for estimating a generation coefficient for simultaneously satisfying the time-sequence data of the wild strain and the specific destroyed strain in the gene model by the number of individual genes to be estimated which are included in the wild strain is included, and gene expression coefficient matrix estimation steps (S16-S20) in which a repeating loop for successively performing the individual gene expression coefficient estimation steps by extracting succeeding destroyed strain after finishing the estimation of the specific destroyed strain by the prescribed number of destroyed strains is included. An initial value for estimating a gene expression coefficient matrix is obtained through the two steps. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

  The present invention relates to a method for estimating the amount of change action, in which each gene changes by being influenced by time from other genes among a plurality of genes.

Conventionally, in such a plurality of gene groups, individual genes are influenced by other genes, that is, as a method for estimating gene interaction, a threshold test, a Bayesian network using a statistical method, and a binomial relationship are used. There are various technologies such as a Boolean network and a method using a differential equation.
In particular, an estimation method using a differential equation model that can identify feedback control of an equivalence class is well known.

  The applicant has not obtained prior art documents at the time of filing.

In the conventional method, the initial value is given to the coefficient of action of each gene for estimation, but the estimation is started. There is a problem that it takes time.
In addition, there is a problem that the amount of calculation becomes enormous, or that it is necessary to add a restriction in order to prevent the divergence of numerical values in iterative estimation, and therefore it is necessary to limit the number of genes.

  The present invention has been made in order to solve the above-described problems, and provides a gene expression interaction estimation method in which an initial value of an appropriate coefficient can be easily obtained and a calculation amount is suppressed.

The gene expression action estimation method according to the present invention is a method of creating a gene model and estimating a gene expression action based on time series data of a wild strain consisting of an unmanipulated gene group and a disrupted strain in which a specific gene is destroyed. ,
Expression coefficient for individual genes with a loop that repeats the number of individual genes to be estimated included in the wild strain for the gene model to estimate the expression coefficient that simultaneously satisfies the time series data of the wild strain and the specific disrupted strain An estimation step;
When the estimation with a specific disrupted strain is completed, the next step is to take the next disrupted strain and repeat the expression coefficient estimation step for individual genes. And with
Through the above two steps, an estimated initial value of the gene expression coefficient matrix was obtained.

  Since the initial value of the gene interaction matrix is estimated from the combination of the time course values of the wild strain and the disrupted strain, the time until the convergence of the estimation is shortened and there is an effect that stable and reasonable estimation can be performed.

Embodiment 1 FIG.
In order to estimate gene interaction or the degree of influence from other genes, a semi-optimization method using the least square method is applied to obtain an initial value, and a quasi-optimization using a linear differential equation model is applied to the initial value. A method for reducing the constraints and reducing the estimated time by applying optimization will be described.
In the present embodiment, the following equation (1) is set as a linear differential equation model.

  Here, x is a gene expression level vector, A is a coefficient matrix indicating an interaction between genes, and w is a bias coefficient vector based on uncertain elements. The uncertain element may be a polynomial function of time, but is described here as a bias model. Further, from the above equation (1), in the least square method, the gene expression level vector is represented by the following equation (2) obtained by integrating both sides of the equation (1).

x ₀ is the gene expression amount initial value (initial value vector), t denotes time.
The least square method is set so that each gene Y changes in the following equation (3) based on the time course data that changes with time for one gene of the gene expression level vector (gene vector) x. Then, in the time change, the matrix B indicating the interaction from other genes and W are estimated.
Y = Y ₀ + B · X + W · t (3)
Here, Y is time course data of one gene in the target gene group, and X is a gene group as an aggregate of individual genes Y, and is a time integral value of the gene expression level x. Y ₀ is the initial value.
FIGS. 1 and 2 show a flow of the estimation method based on the model set in this way. FIG. 3 shows an example of a calculation circuit configuration for developing this operation. In FIG. 3, as a circuit configuration for estimation, an input device 51 that reads actual data, a display device 52 that obtains the degree of influence between genes, which is an analysis result of this device, as a pathway, and displays this pathway. And a gene expression action estimation device having a hardware configuration for performing this processing between the input and output. This gene expression action estimation apparatus can be realized by a computer having a central processing unit. The detailed configuration of the gene expression action estimation apparatus has the following components. In order from the input side, a time course / data input processing unit 41, an interpolation processing unit 42 for interpolating the data obtained thereby, a normalization processing unit 43 for normalizing from the interpolation result, an initial value estimation for estimating the coefficient initial value of the model Unit 44, an interaction coefficient estimation processing unit 45 that estimates an interaction coefficient based on the obtained initial value, a coefficient optimization processing unit 46 that performs further optimization processing, and a time course that stores input time course data. A data storage unit 48, an analysis method storage unit 47 that stores a gene expression model that expresses a mathematical expression of the model, and a pathway drawing processing unit 49 that graphically illustrates the effect between the obtained genes, each with hardware having a corresponding program is there. Or each said part may be implement | achieved only by a program, may be implement | achieved only by hardware, and may be implement | achieved by the combination of a program and hardware.

An estimation method using this computer will be described with reference to FIGS.
First, the time change data of the individual gene Y of interest is read from the experimental value which is time course data. FIG. 4 shows observed values obtained by experiments and the like. In the figure, the leftmost column is an index for explanation, and 1 to 27 are shown in the figure. Further, in this figure, an example of four genes (G01, G02, G03, G04) is shown for easy explanation. Indexes 5 to 7 show time-change data of wild strains, where a gene group composed of natural individual genes that have not been manipulated is a wild strain. The leftmost column represents the time, in this example, three times (1.1, 1.2, 1.3), and the next column (1.11, 1.21, 1.31) represents the gene G01. The time course values at the left time are the following 1.12 to 1.32, the time course value of the sample time of the gene G02, 1.13 to 1.33 are the time course values of G03, 1.14 Thru 1.34 is the time course value of G04. Furthermore, as shown in index 8, indexes 10 to 12 are time course values at three times (1.1, 1.2, 1.3) of the disrupted strain in which the G01 gene is disrupted. The lower columns of 12, 2.13, and 2.14 are time change values of the genes G01, G02, G03, and G04, respectively. In the following, indexes 15 to 17 are the time course values of the disrupted strain in which the gene G02 is disrupted, 20 to 22 are the disrupted strain in which the G03 is disrupted, and 25 to 27 are the disrupted strains in which the G04 is disrupted. The time course / data input processing unit 41 performs this in S11.

In S12, the computer normalizes the time course data with data having a large absolute value for computer processing, and in S13 performed by the interpolation processing unit 42, the intermediate value of the elapsed time in FIG. 4 is splined (smoothed) based on the logarithmic value. ) Interpolate. In S14 performed by the normalization processing unit 43, the interpolated logarithmic value data is returned to the original physical quantity.
In S15, the matrix X and the vector Y are generated by integrating the data so far for each gene at each sampling time. In S15, in addition to this, in the generated matrix X, the elements of the interaction that are not estimated are fixed to zero, and the significant digits of the X matrix and the Y vector are adjusted for calculation accuracy and stabilization.

A loop characteristic of this embodiment will be described. In S16, a large loop surrounded by S21 is formed. This large loop loops over time of the wild strain and one disrupted strain (in FIG. 4, for example, a disrupted strain that has disrupted G01 noted in index 8 and its time course value is shown in indexes 10 to 12). This represents the estimation of B and W for all combinations, such as a combination of a wild strain and the next disrupted strain (destroyed strain that destroyed G02).
The small loop surrounded by S17 and S20 represents that B and W are obtained from the time course values of the disrupted strain and the wild strain that destroyed G01, for example. That is, in S17, the first wild-type gene G01 is picked up, and from the time course values of indexes 5 to 7 and 10 to 12 in FIG. 4, paying attention to G01, the least square method is obtained in S18. Provisionally, the coefficient of interaction regarding the gene G01 in the matrix B and the vector W in the equation (3) is estimated. Thus, in the example of FIG. 4, in the estimation using the time course values of the indexes 5 to 7 and the time course values of 10 to 12, there are four genes, so four B and W are estimated. That is, in this small loop, the estimation using the indexes 10 to 12 in FIG. 4 is performed. In the first loop, the estimation focusing on G01 is performed, and in the second small loop, the estimation focusing on G02 is performed. In the small loop, estimation is performed by paying attention to G03. Finally, in the fourth small loop, estimation is performed focusing on g04. As a result, in S19, using the obtained counts in B and W, the values of the matrices A and w in the equation (1) are temporarily estimated.

Next, the process returns to S16, and the loop of S17 to S19 is repeated for the other disrupted strains in the time course data (in the example of FIG. 4, the other disrupted strains indicated by indexes 13, 18, and 23). In other words, a large loop represents a shift to estimation using other disrupted strains and wild strains. The initial value estimation unit 44 performs processing of a double loop including this large loop.
In this way, the disrupted strains are sequentially taken out, and B and W that are satisfied simultaneously with the disrupted strain and the wild strain are estimated. As a result, the interaction that all the disrupted strains give to all genes can be estimated. At this time, the gene Y is a vector that stores time-series data of the expression level of the gene included in the loop. X stores the integrated value of the expression level of the gene Y vector at the tag time.
When the estimation of the B and W coefficients of the genes affected by all the disrupted strains is completed in S21, the interaction coefficient estimation processing unit 45 uses the remaining A and w coefficients as wild strains for the coefficients that are not affected by the disrupted strains in S23. The time course data is input and estimated by the method of least squares. In this way, the coefficient matrix A and the vector w are provisionally determined as initial values for semi-optimization.

After S24, the coefficient optimization processing unit 46 performs a process for determining the final coefficient. In S24, two or more models, for example, five models are generated as linear differential equation models by slightly varying the coefficients of the elements of A and w for the model of x tentatively determined in S23. To do.
In S25, the maximum number of trials is set in order to limit the calculation time until confirmation, and in S26, the first model is extracted from the five models generated in S24, and the processes of S27 and S28 are performed. In S27, the differential equation of the model is numerically integrated to obtain the time course value of each gene of each individual.
In S28, the time course value obtained by the integration is compared with the corresponding input experimental value in S11 to obtain an estimated relative error.
Returning from S29 to S26, the process of S27 and S28 is repeated for all five models generated in S24, and the estimated relative error of each model is obtained.

If there is a model having a predetermined error or less in S30, the estimation is finished assuming that the model is an A estimation matrix and a w vector.
If none of the models is equal to or smaller than the predetermined error in S30, the model having the smallest error among the five models is selected as an individual in S31, and five models are generated again. In S32, the different coefficients of the matrix A and the vector w are slightly changed for the five models. In S33, if the error does not become the threshold value or less in the predetermined number of times S30 within the maximum number of trials, the minute fluctuation in S32 is reduced.
If a model having a predetermined error or less is likely to be obtained in S30 by the trial loop surrounded by S25 and S34, and a model with a small error is determined, the matrix A and the vector w represented by the model are determined as the final estimation result. Become.

FIG. 5 is a diagram in which the final model obtained in this way and obtained with a small error in S30 is represented by an estimated value and compared with the time course experimental value of FIG.
FIG. 6 is a schematic diagram showing the interaction between genes, that is, the influence relationship to other genes, obtained as a result of estimation. This figure is obtained by the pathway drawing processing unit 49. As a result of estimation according to FIGS. 1 and 2, B is determined, and gene C3 receives a positive interaction from gene B6, receives a negative action from B7, and forms a positive action on C1 in combination with C2. Presumed.
In this way, the gene is modeled with a differential equation, and the initial value of the gene interaction matrix of the differential equation model is estimated from the combination of the time course values of the wild strain and the disrupted strain, so the number of simulations until the estimated convergence is shortened And a stable and reasonable estimate.
It should be noted that a program having the steps of FIGS. 2 and 3 as a procedure is prepared, installed in another computer by a recording medium or online loading, and the program is executed by the central processing unit of the computer to perform the above analysis. You may do it.

It is a figure which shows the matrix and vector estimation processing flow in Embodiment 1 of this invention. FIG. 10 is a diagram showing a continuation of the matrix and vector estimation processing flow in the first embodiment. 3 is a diagram showing a circuit configuration for performing matrix and vector estimation processing in Embodiment 1. FIG. 6 is a diagram illustrating an example of a time course value obtained by an experiment in Embodiment 1. FIG. It is a figure which shows the comparison with the time course value based on the estimated value obtained in Embodiment 1 and the time course value of an experimental value. 3 is a pathway diagram showing an interaction between genes based on an estimation result in Embodiment 1. FIG.

Explanation of symbols

  S16 Step of picking up specific disrupted strains, S17 Step of selecting individual genes to be estimated, S18 Step of estimating matrix B by coefficient of least squares, vector W, S19 Step of registering elements in matrix A and vector w, S20 specification A step of checking whether there is an unestimated gene in the disrupted strain, a step of checking whether there is an S21 unestimated disrupted strain, S17, S20 A small one that sequentially estimates the coefficient of each individual gene among the specific disrupted strains Loop, S16, S20 large loop for estimating each disrupted strain sequentially, S22 step for confirming whether estimation has been completed for all disrupted strains, S23 step for estimating unestimated coefficients of matrices A and w, S24 a plurality of linears to be estimated Differential equation model creation step, S26 Select one linear differential equation model created in S24 Step S27, obtaining the time series value of each gene by integrating the model, S28 calculating the error between the time integral value and the input gene time series value, and taking the linear differential equation model created in S29 S24 in order, A step of confirming whether the model has been completed, a step of checking whether there is a model having a predetermined error or less in the models created in S30 and S24, and a step S31 for the number of models created in S24 with respect to the model having the smallest error. Step, step S32, changing the values of the elements of the matrix A and the vector w differently for each model, step S33, reducing the variation rate, step S34, checking whether the set maximum number of trials has been exceeded, step 41 Data input processing unit, 42 Interpolation processing unit, 43 Normalization processing unit, 44 Initial value estimation unit, 4 Interaction coefficient estimation processing unit, 46 coefficient optimization processing unit, 47 analysis method storage unit, 48 Time-course data storage unit, 49 pathway drawing processing unit.

Claims (9)

In a method for estimating a gene expression action based on time series data of a wild strain consisting of a gene group that has not been manipulated and a disrupted strain in which a specific gene has been destroyed,
To individual genes having a loop that repeats the number of the estimation target individual genes included in the wild strain for the gene model to estimate the expression coefficient that simultaneously satisfies the time series data of the wild strain and the specific disrupted strain An expression coefficient estimation step of
When the estimation with the specific disrupted strain is completed, the gene expression coefficient matrix having a repeated loop is sequentially performed for a predetermined number of disrupted strains by taking the next disrupted strain and repeating the expression coefficient estimation step for the individual gene. An estimation step,
A gene expression action estimation method characterized in that an estimated initial value of a gene expression coefficient matrix is obtained through the above two steps.   The gene expression action estimation method according to claim 1, wherein the gene model is estimated by approximating the gene expression with a time integration of a gene expression level having a gene expression coefficient matrix.   2. The gene expression action estimation method according to claim 1, wherein, in the expression coefficient estimation step for an individual gene, a least square method is applied to estimate the expression coefficient. In a device that creates a gene model and estimates the gene expression effect based on time series data of a wild strain consisting of a group of genes that has not been manipulated and a disrupted strain that has destroyed a specific gene
With respect to the gene model, a procedure for estimating an expression coefficient that simultaneously satisfies the time series data of the wild strain and a specific disrupted strain is calculated as a loop that repeats the number of individual genes to be estimated included in the wild strain, When the estimation with a specific disrupted strain is completed, an initial coefficient estimator is provided that sequentially takes out a predetermined number of disrupted strains by taking the next disrupted strain and repeating the expression coefficient estimation calculation for the individual genes. Characteristic gene expression action estimation device.   5. The gene expression action estimation apparatus according to claim 4, wherein the coefficient initial value estimation unit estimates the gene model by approximating the gene expression by time integration of the gene expression level having the gene expression coefficient matrix.   The gene expression action estimation device according to claim 4, wherein the coefficient initial value estimation unit applies a least square method to the estimation of the expression coefficient. In a program that makes a computer create a gene model and estimate the gene expression effect based on time series data of a wild strain consisting of a group of genes that have not been manipulated and a disrupted strain that has destroyed a specific gene,
To individual genes having a loop that repeats the number of the estimation target individual genes included in the wild strain for the gene model to estimate the expression coefficient that simultaneously satisfies the time series data of the wild strain and the specific disrupted strain An expression coefficient estimation procedure for
When the estimation with the specific disrupted strain is finished, the next disrupted strain is taken and the expression coefficient estimation procedure for the individual gene is repeated, sequentially for a predetermined number of disrupted strains, and a gene expression coefficient matrix having an iterative loop A gene expression action estimation program comprising a procedure for obtaining an estimated initial value of a gene expression coefficient matrix using an estimation procedure.   8. The gene expression action estimation program according to claim 7, wherein the expression coefficient estimation procedure estimates a gene model by approximating the gene expression level having a gene expression coefficient matrix by time integration.   8. The gene expression action estimation program according to claim 7, wherein the expression coefficient estimation procedure applies a least square method to the expression coefficient estimation.

Images