JP4819382B2 - Expression rate expression generation device, expression rate expression generation method, and expression rate expression generation program - Google Patents

Description

  The present invention relates to an expression rate expression generation device, an expression rate expression generation method, and an expression rate expression generation program in bioinformatics for estimating gene expression interaction.

Elucidation of the expression interaction of each gene between genes is expected, and one of the estimation methods is a method of applying an S-System model. In the method to which the S-System model is applied, in order to reproduce given observation data (time course data), an expression indicating a gene expression rate is expressed by a term indicating a gene generation rate and a term indicating a decomposition rate ( Non-patent document 1).
Kenichi Matsubara and Yoshiyuki Tsuji, “Genome Information Biology”, Nakayama Shoten, July 10, 2001, pp. 165-188

  The gene expression interaction estimation method based on the S-System model requires enormous calculation as both the decomposition system and the generation system operate simultaneously. In addition, in order to estimate the gene expression interaction, time course data for all genes to be estimated is required.

  The present invention has been made to solve the above-described problems, and an object thereof is to estimate expression interaction between genes.

  The expression rate expression generation apparatus of the present invention is configured to determine a gene involved in the degradation of the first gene based on the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene. Under the condition that only the first gene is used, a degradation system coefficient calculation unit that calculates a degradation system coefficient for the first gene and stores it in the storage unit, and expression of each gene in the non-destructive strain for the first gene And a generation system coefficient calculation unit that calculates a generation system coefficient for the first gene based on non-destructive strain time course data that is the amount of data observed, and stores the generation system coefficient in a storage unit.

  In addition, the degradation system coefficient calculation unit, when the data amount of the disrupted strain time course data is insufficient with respect to the number of genes of interest expressing the expression interaction by the expression rate formula, The degradation system coefficient for the first gene is calculated under the condition that only the first gene is involved in the degradation.

  The expression rate expression generating apparatus of the present invention calculates a degradation system coefficient for the first gene based on the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene. Based on non-destructive strain time course data, which is data obtained by observing the expression level of each gene for the non-destructive strain for the first gene, and for the generation of the first gene A generation system coefficient calculation unit that calculates a generation system coefficient for the first gene and stores the generation system coefficient in a storage unit under a condition that limits the number of genes involved to a specific number or less.

  In addition, the generation system coefficient calculation unit is configured to output the first gene when the data amount of the non-destructive strain time course data is insufficient with respect to the number of target genes representing the expression interaction by the expression rate equation. A generation system coefficient for the first gene is calculated under a condition that limits the number of genes involved in the generation of the gene to a specific number or less.

  The expression rate expression generation device of the present invention determines the genes involved in the degradation of the first gene based on the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain of the first gene. Under the condition that only one gene is used, a degradation system coefficient calculation unit that calculates a degradation system coefficient for the first gene and stores it in the storage unit, and an expression level of each gene for a non-destructive strain of the first gene Based on the nondestructive strain time course data, which is the observed data, the generation system coefficient for the first gene is calculated under the condition that the number of genes involved in the generation of the first gene is limited to a specific number or less. And a generation system coefficient calculation unit stored in the storage unit.

  In addition, the degradation system coefficient calculation unit, when the data amount of the disrupted strain time course data is insufficient with respect to the number of genes of interest expressing the expression interaction by the expression rate formula, Under the condition that the gene involved in degradation is only the first gene, a degradation system coefficient for the first gene is calculated, and the generation system coefficient calculation unit is a target for expressing the expression interaction by the expression rate equation. When the data amount of the non-destructive strain time course data is insufficient with respect to the number of genes, the first is under the condition that the number of genes involved in the generation of the first gene is limited to a specific number or less. It is characterized in that a production system coefficient for the gene of is calculated.

  Further, the expression rate expression generation device inputs the destroyed strain time course data and the non-destructive strain time course data, and calculates the spline interpolation calculation for the logarithmic value of each time course value indicated by the destroyed strain time course data. The first derivative is calculated and stored in the storage unit, and the second derivative is calculated and stored for the function obtained by the spline interpolation calculation for the logarithmic value of each time course value indicated by the non-destructive stock time course. The decomposition system coefficient calculation unit acquires the first differential coefficient calculated by the preprocessing unit based on the destroyed stock time course data from the storage unit, and acquires the first A decomposition system coefficient for the first gene is calculated based on the derivative, and the generation system coefficient calculation unit is configured to calculate the second derivative coefficient calculated by the preprocessing unit based on non-destructive strain time course data. Acquired from the storage unit, and calculates a production system coefficients for the first gene based on the obtained second derivative.

In addition, the expression rate expression generation device uses n as the number of target genes that express the expression interaction in the expression rate expression, and the expression level of the i-th gene (1 ≦ i ≦ n) corresponding to the first gene _Is the i-th gene production system coefficient and the coefficient to be multiplied by the term indicating the production rate is α _i , the i-th gene production form coefficient and the j-th gene for the generation of the i-th gene (1 ≦ j ≦ n) ) _Is the decomposition coefficient of the i-th gene and the coefficient to be multiplied by the term indicating the decomposition rate is β _i , the decomposition system coefficient of the i-th gene and the j-th gene When the interaction coefficient is h _ij , the generation system coefficients α _i and g _ij of the i th gene and the degradation system coefficients β _i and h _{ij of the} i th gene are calculated with respect to the expression rate expression expressed by the following formula: It is characterized by doing.

  Furthermore, the expression rate equation generation device inputs a decomposition system coefficient for the first gene calculated by the decomposition system coefficient calculation unit and a generation system coefficient for the first gene calculated by the generation system coefficient calculation unit, The expression rate equation is generated using the degradation system coefficient and the generation system coefficient for the input first gene, the expression level of the first gene is calculated based on the generated expression rate equation, and the generation system The evaluation result of the production system coefficient for the first gene calculated by the coefficient calculation unit is the comparison result between the expression level of the first gene indicated by the non-destructive strain time course data and the calculated expression level of the first gene. A first evaluation unit that outputs based on the first evaluation unit is provided.

  Further, the expression rate expression generation device includes a second evaluation unit that evaluates a gene involved in the generation of the first gene, and the generation system coefficient calculation unit is involved in the generation of the first gene. Gene sets are generated by limiting the number of genes to a specific number or less, a generation coefficient for the first gene is calculated for each gene set, and the first evaluation unit is the generation system coefficient calculation unit An evaluation value is calculated as an evaluation result for each generative coefficient corresponding to each gene set calculated by the first evaluation unit, and the second evaluation unit corresponds to each gene corresponding to each evaluation value calculated by the first evaluation unit. And the evaluation result for the gene involved in the generation of the first gene is output based on the total value for each gene.

  Further, the expression rate expression generation device calculates the expression level of each gene by combining the expression rate formulas representing the expression rate of each gene for a plurality of genes, and based on the calculated expression level of each gene, It has a simultaneous evaluation unit that evaluates the production system coefficient used in the expression rate formula, and the degradation system coefficient calculation unit inputs the time course data of each disrupted strain for each gene and calculates the degradation system coefficient for each gene. The generation system coefficient calculation unit inputs each non-destructive strain time course data for the non-destructive strain of each gene and calculates a generation system coefficient for each gene, and the simultaneous evaluation unit calculates the decomposition system coefficient calculation unit The decomposition system coefficient for each gene calculated by and the generation system coefficient for each gene calculated by the generation system coefficient calculation unit are input, and the input decomposition system coefficient and generation system coefficient for each gene are used. Each expression rate formula is generated, and the expression level of each gene is calculated by combining the generated expression rate formulas. The evaluation results for the generation system coefficient for each gene calculated by the generation system coefficient calculation unit are all strains ( It outputs based on the comparison result of the expression level of each gene which the time course data of the said destroyed strain and the said non-destructive strain show, and the calculated expression level of each gene.

The generation system coefficient calculation unit calculates a plurality of generation system coefficients for the first gene under a condition that limits the number of genes involved in the generation of the first gene to a specific number or less, and the simultaneous evaluation The unit selects one gene from a plurality of genes, registers the selected gene as a simultaneous gene that is a gene corresponding to the simultaneous expression rate expression, and selects one of a plurality of generation system coefficients for the selected gene. A first registration unit that selects a generation system coefficient and registers the generation system coefficient selected as a simultaneous generation system coefficient that is a generation system coefficient used in the simultaneous expression rate equation; and the first registration unit serves as a simultaneous gene One gene is selected from genes other than the registered genes and registered as temporary simultaneous genes, and two or more generating system coefficients for each of the temporary simultaneous genes are temporarily generated. Select as a coefficient, using the decomposition system coefficient for the simultaneous gene, the simultaneous generation system coefficient, the decomposition system coefficient for the temporary simultaneous gene, and each temporary simultaneous generation system coefficient, and the expression rate formula for the simultaneous gene and each of the temporary simultaneous gene An expression rate expression, and the expression rate expression for the generated simultaneous gene is associated with each expression rate expression for the temporary simultaneous gene to express the expression level of the simultaneous gene and the expression level of the temporary simultaneous gene. Calculated for each temporary simultaneous generation system coefficient, and the expression level of the calculated simultaneous genes and the temporary simultaneous genes and the simultaneous genes and the temporary simultaneouss indicated by the time course data of all strains (the disrupted strain and the non-destructed strain) An evaluation value is calculated for each temporary simultaneous generation coefficient based on the result of comparison with the expression level with the gene, and the calculated evaluation value is compared with the specific value. When at least one of the evaluation values is equal to or greater than the specific value, the temporary simultaneous gene is added and registered as a simultaneous gene, and the temporary simultaneous generation coefficient corresponding to the highest evaluation value is used as the simultaneous generation system coefficient. A second registration unit for additionally registering,
The simultaneous evaluation unit executes the second registration unit for all unselected genes as temporary simultaneous genes after executing the first registration unit, and outputs the generation system coefficient registered as the simultaneous generation system coefficient. Features.

  Furthermore, the simultaneous evaluation unit selects one gene registered as a simultaneous gene, selects one generation system coefficient that is a generation system coefficient for the selected gene and is not registered as a simultaneous generation system coefficient, and selects the selected gene The evaluation value for the generative coefficient selected for the gene is compared with the evaluation value for the generative coefficient registered as a simultaneous generative coefficient for the selected gene, and the evaluation value of the selected generative coefficient is higher A generating system coefficient updating unit that updates the simultaneous generating system coefficient registered for the selected gene with the selected generating system coefficient.

  Further, the simultaneous evaluation unit selects one gene that is not registered as a simultaneous gene, a degradation system coefficient for the simultaneous gene, a simultaneous generation system coefficient, a degradation system coefficient for an unregistered gene as a simultaneous gene, and an unrecognized gene as a simultaneous gene. For each generative coefficient for the registered gene, the expression rate expression for the simultaneous gene and the unregistered gene as a simultaneous gene is used to register each gene expression coefficient and the unregistered as the simultaneous gene expression level. As the simultaneous gene and simultaneous gene, the expression level of the calculated simultaneous gene and the unregistered gene as a simultaneous gene and the time course data of all strains (destructed strain and non-destructed strain) are calculated. Based on the result of comparison with the expression level with the unregistered gene, an evaluation value is calculated for each generation coefficient of the unregistered gene as a contiguous gene, Compare the evaluation value and the specific value, and if at least one of the evaluation values is equal to or higher than the specific value, add the unregistered gene as a simultaneous gene and register it as a simultaneous gene. A generation system coefficient additional registration unit that additionally registers a generation system coefficient of a corresponding unregistered gene as a simultaneous generation system coefficient is provided.

  The expression rate expression generation method of the present invention is based on the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene. Under the condition that only the first gene is used, the degradation system coefficient for the first gene is calculated and stored in the storage unit, and the expression level of each gene is observed for the non-destructive strain for the first gene. A generation system coefficient for the first gene is calculated based on certain non-destructive strain time course data and stored in a storage unit.

  The expression rate formula generation method of the present invention calculates the degradation system coefficient for the first gene based on the time course data of the disrupted strain, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene. Based on the non-destructive strain time course data that is stored in the storage unit and the expression level of each gene is observed for the non-destructive strain for the first gene, the number of genes involved in the generation of the first gene is calculated. A generation system coefficient for the first gene is calculated and stored in a storage unit under a condition limited to a specific number or less.

  The expression rate expression generation method of the present invention is based on the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain of the first gene. Under the condition that only one gene is used, the degradation system coefficient for the first gene is calculated and stored in the storage unit, and the non-destructive strain of the first gene is data obtained by observing the expression level of each gene. Based on the disrupted strain time course data, the generation system coefficient for the first gene is calculated and stored in the storage unit under the condition that the number of genes involved in the generation of the first gene is limited to a specific number or less. It is characterized by that.

The expression rate expression generation program of the present invention uses the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene, to determine the genes involved in the degradation of the first gene. Under the condition that only the first gene is used, the degradation system coefficient for the first gene is calculated and stored in the storage unit, and the expression level of each gene is observed for the non-destructive strain for the first gene. The computer is caused to calculate a generation system coefficient for the first gene based on certain non-destructive strain time course data and store it in the storage unit.
.

  The expression rate expression generation program of the present invention calculates a degradation system coefficient for the first gene based on the time course data of the disrupted strain, which is data obtained by observing the expression level of each gene for the disrupted strain for the first gene. Based on the non-destructive strain time course data that is stored in the storage unit and the expression level of each gene is observed for the non-destructive strain for the first gene, the number of genes involved in the generation of the first gene is calculated. Under a condition limited to a specific number or less, the computer is caused to calculate the generation coefficient for the first gene and store it in the storage unit.

  The expression rate expression generation program of the present invention determines the genes involved in the degradation of the first gene on the basis of the disrupted strain time course data, which is data obtained by observing the expression level of each gene for the disrupted strain of the first gene. Under the condition that only one gene is used, the degradation system coefficient for the first gene is calculated and stored in the storage unit, and the non-destructive strain of the first gene is data obtained by observing the expression level of each gene. Based on the disrupted strain time course data, the generation system coefficient for the first gene is calculated and stored in the storage unit under the condition that the number of genes involved in the generation of the first gene is limited to a specific number or less. Let the computer do that.

According to the present invention, in the gene expression interaction estimation apparatus, the degradation system coefficient calculation unit calculates the degradation system coefficient under the condition that the gene involved in the degradation is only its own gene, and the generation system coefficient calculation unit By generating generative coefficients under conditions that limit the number of genes involved in generation to a specific number or less, gene expression rate formulas can be generated even when the amount of time course data is insufficient it can.
Further, the first evaluation unit evaluates the generation system coefficient, so that a more suitable generation system coefficient can be selected based on the evaluation result to generate a gene expression rate expression.
In addition, the simultaneous evaluation unit can evaluate the expression rate expression of each gene in combination, thereby generating a gene expression rate expression by selecting a more suitable generation coefficient based on the evaluation result.
Moreover, the 2nd evaluation part can estimate the gene concerned in the production | generation of each gene based on the evaluation result by collectively evaluating the evaluation result of a production | generation system coefficient for every gene.

Reference Example 1
The gene expression interaction estimation method using the S-System model is expressed as follows: “Expression rate change rate (expression rate) = production rate−decomposition” with a formula for reproducing time course data as a differential equation of gene expression level (expression rate formula) Speed [formula 1] ".

Then, the generation system coefficients α _i and g _ij and the decomposition system coefficients β _i and h _ij in Expression 1 are calculated.
α _i is a generation system coefficient of the i-th gene and is a coefficient to be multiplied by a term indicating the generation speed. G _ij is a generation coefficient of the i-th gene, and is an interaction coefficient of the j-th gene with respect to the generation of the i-th gene. Β _i is a decomposition system coefficient of the i-th gene and is a coefficient by which a term indicating the decomposition rate is multiplied. H _ij is a decomposition system coefficient of the i-th gene, and is an interaction coefficient of the j-th gene with respect to the decomposition of the i-th gene.

Here, the following method (Japanese Patent Application No. 2004-101773) may be taken as a method for calculating the optimum value of each coefficient of Equation 1 without requiring enormous calculation.
In this method, using the time course data of the disrupted strain, using the value converted to the logarithmic space and the derivative obtained from the interpolation data, the optimum value of the decomposition system coefficients β _i and h _ij is obtained by the least square method. calculate. Also, assuming the decomposition system coefficient, the time course data of the non-destructive strain data is used, the derivative and the generation speed profile are calculated from the value converted into the logarithmic space and the interpolation data, and then used to calculate the minimum two The optimum values of the generation system coefficients α _i and g _ij are calculated by multiplication.
The procedure of this method will be described below.

FIG. 1 is a flowchart showing a processing flow of the gene expression interaction estimation method in Reference Example 1.
The procedure will be described with reference to FIG.

  First, the time course data of a wild strain (non-destructive strain) and a disrupted strain are input (S101).

  Next, as preprocessing, calculation of a spline interpolation coefficient based on the logarithmic value of the time course value and calculation of a differential coefficient of logarithmic value data using the spline interpolation coefficient and the time course value are performed (S102). ).

In calculating the spline interpolation coefficient, when the logarithmic value of the time course of the gene i of the k strain is Z _i ^(k) (t _j ) (j = 1, 2,... M), each t _j and t A coefficient (spline interpolation coefficient) that can be expressed by one cubic polynomial between _{j + 1} is obtained. The coefficient to be obtained is a value that passes through each point of Z _i ^(k) (t _j ) and is smoothly connected at the switching point t _j of the polynomial when represented by a graph using a cubic polynomial as a function.

If the coefficient can be obtained as described above and spline interpolation can be performed, the function can be expressed by a time polynomial, and the polynomial can be differentiated. Each coefficient after differential calculation is a derivative.
The derivative is calculated for each of the wild strain and the disrupted strain.

Next, the decomposition system coefficients β _i and h _ij are estimated (S103).

In the estimation of the degradation system coefficients β _i and h _ij , first, for the data of the disrupted strain of gene i,

The following vector b ₁ and matrix A ₁ are obtained.

Then, using the matrix A ₁ and the vector b ₁ , a vector c ₁ such that “b ₁ = A ₁ c ₁ ” is estimated by the least square method, and the decomposition system coefficients β _i and h _ij for the gene i are obtained. calculate. The relationship between the decomposition system coefficients β _i , h _ij and the vector c ₁ is as follows.

The above processing is performed for all genes, and the degradation system coefficients β _i and h _ij of all genes are estimated.

Next, the generation system coefficients α _i and g _ij are estimated (S104).

In the estimation of the production system coefficients α _i and g _ij , first, for the data of the non-destructive strain of gene i,

The following vector b ₂ and matrix A ₂ are obtained.

Then, using the matrix A ₂ and the vector b ₂ , a vector c ₂ such that “b ₂ = A ₂ c ₂ ” is estimated by the least square method, and the generation system coefficients α _i and g _ij for the gene i are obtained. calculate. The relationship between the generation system coefficients α _i and g _ij and the vector c ₂ is as follows.

The above processing is performed on all genes, and the generative coefficients α _i and g _ij of all genes are estimated.

  In the least square method, the vector c having the smallest value of “| b−Ac |” is calculated, and the optimum value of each coefficient is estimated.

  With the procedure described with reference to FIG. 1, it is possible to calculate the optimum value of each coefficient of Equation 1 without requiring enormous calculation.

Embodiment 1 FIG.
In the case of the method of Reference Example 1 above, time course data of the disrupted strain is required for all genes to be analyzed, but it is difficult to prepare time course data of the disrupted strain for all genes from the beginning. is there.
Therefore, in the first embodiment, a solution for a case where a disrupted strain does not exist or a case where the time course data of the disrupted strain is insufficient and the least square method cannot be applied, and a gene expression interaction estimation apparatus 100 for executing the method. Will be described. As a result, the solution (each coefficient) can be calculated regardless of the presence or absence of time course data of the destroyed stock.

FIG. 2 is a diagram illustrating an appearance of the gene expression interaction estimation apparatus 100 according to the first embodiment.
2, a gene expression interaction estimation apparatus 100 includes a system unit 910, a CRT (Cathode Ray Tube) display device 901, a keyboard (K / B) 902, a mouse 903, a compact disk device (CDD) 905, a printer device 906, A scanner device 907 is provided, and these are connected by a cable.
Further, the gene expression interaction estimation apparatus 100 is connected to a FAX machine 932 and a telephone 931 via a cable, and is connected to the Internet 940 via a local area network (LAN) 942 and a web server 941.

FIG. 3 is a hardware configuration diagram of the gene expression interaction estimation apparatus 100 according to the first embodiment.
In FIG. 3, the gene expression interaction estimation apparatus 100 includes a CPU (Central Processing Unit) 911 that executes a program. The CPU 911 includes a ROM 913, a RAM 914, a communication board 915, a CRT display device 901, a K / B 902, a mouse 903, an FDD (Flexible Disk Drive) 904, a magnetic disk device 920, a CDD 905, a printer device 906, and a scanner device 907. Connected with.
The RAM 914 is an example of a volatile memory. The ROM 913, the FDD 904, the CDD 905, the magnetic disk device 920, and the optical disk device are examples of nonvolatile memories. These are examples of a storage device or a storage unit.
The communication board 915 is connected to a FAX machine 932, a telephone 931, a LAN 942, and the like.
For example, the communication board 915, the K / B 902, the scanner device 907, the FDD 904, and the like are examples of the information input unit.
Further, for example, the communication board 915, the CRT display device 901, and the like are examples of the output unit.

Here, the communication board 915 is not limited to the LAN 942 but may be directly connected to the Internet 940 or a WAN (Wide Area Network) such as ISDN. When directly connected to the Internet 940 or a WAN such as ISDN, the gene expression interaction estimation apparatus 100 is connected to the Internet 940 or a WAN such as ISDN, and the web server 941 is unnecessary.
The magnetic disk device 920 stores an operating system (OS) 921, a window system 922, a program group 923, and a file group 924. The program group 923 is executed by the CPU 911, the OS 921, and the window system 922.

The program group 923 stores programs that execute functions described as “˜units” in the description of the embodiments described below. The program is read and executed by the CPU 911.
In the file group 924, in the description of the embodiment described below, “determined to”, “result of determining to”, “calculated to”, “calculated result of”, and “to” are processed. The result information described in terms of expressions such as “Shi” and “Result of processing” is stored as “˜File”.
In addition, the arrow portion of the flowchart described in the description of the embodiment described below mainly indicates input / output of data, and for the input / output of the data, the data is the magnetic disk device 920, FD (Flexible Disk cartridge), Recording is performed on an optical disc, a CD (compact disc), an MD (mini disc), a DVD (Digital Versatile Disk), and other recording media. Alternatively, it is transmitted through a signal line or other transmission medium.

  In addition, what is described as “unit” in the description of the embodiment described below may be realized by firmware stored in the ROM 913. Alternatively, it may be implemented by software alone, hardware alone, a combination of software and hardware, or a combination of firmware.

  In addition, a program for implementing the embodiment described below may be stored using a magnetic disk device 920, an FD, an optical disk, a CD, an MD, a DVD, or other recording media.

FIG. 4 is a configuration diagram of the gene expression interaction estimation apparatus 100 according to the first embodiment.
4, the gene expression interaction estimation apparatus 100 according to Embodiment 1 includes a preprocessing unit 110, a decomposition system coefficient estimation unit 120, a generation system coefficient estimation unit 130, a simultaneous evaluation unit 140, a single evaluation unit 150, and a storage unit 190. Is provided.

The preprocessing unit 110 inputs time course data of the wild strain and the disrupted strain, and calculates a differential coefficient Y 1 — _i ^(k) (t _j ) of _logarithmic data of the time course.
The decomposition system coefficient estimation unit 120 performs an estimation process of the decomposition system coefficients β _i and h _ij .
Generating system coefficient estimating section 130 performs a process of estimating the generation system coefficients alpha _{i, g ij,} and calculates the evaluation value for the estimated generated based coefficients alpha _{i, g ij.}
The simultaneous evaluation unit 140 performs an evaluation process when applying the generation system coefficients α _i and g _ij and the decomposition system coefficients β _i and h _ij for each gene to the differential equation [Equation 1] of the S-System model, Find each coefficient. At this time, the differential equation [Equation 1] for each gene is simultaneously performed for the evaluation process.
The single evaluation unit 150 performs an evaluation process of the degree related to the generation of the gene j with respect to the gene i.
The storage unit 190 stores time course data of wild strains and destroyed strains, data of each obtained coefficient, evaluated degree data, processing result data of each unit, processing variable data, and the like.
The preprocessing unit 110, the decomposition system coefficient estimation unit 120, the generation system coefficient estimation unit 130, the simultaneous evaluation unit 140, and the single evaluation unit 150 input various data from the storage unit 190, and output processing result data to the storage unit 190. .

  The simultaneous evaluation unit 140 includes a simultaneous gene registration unit 141, a generation system coefficient update (mountain climbing) unit 142, and a simultaneous gene addition registration unit 143.

The simultaneous gene registration unit 141 performs evaluation in combination with a differential equation [Equation 1] of an S-System model of another gene based on a certain gene, and a combination of genes whose evaluation value is equal to or greater than a reference value (simultaneous gene ).
The generation system coefficient updating (mountain climbing) unit 142 replaces the generation system coefficient to perform an evaluation process, and updates the generation system coefficient based on the evaluation result.
The simultaneous gene addition registration unit 143 performs an evaluation process by adding another gene to the combined genes, and determines a gene to be added based on the evaluation result.

FIG. 5 is a flowchart showing a process flow of the gene expression interaction estimation method according to the first embodiment.
The flow of processing for estimating the expression interaction of each gene in the first embodiment will be described below with reference to FIG.

  First, the preprocessing unit 110 inputs time course data of a wild strain and a disrupted strain (S201). This process is the same as S101.

  The time course data of the wild strain and the disrupted strain may be input from an external device, or may be stored in the storage unit 190 in advance.

Next, the preprocessing unit 110 performs preprocessing based on the input time course data of the wild strain and the disrupted strain. Then, the calculated data of the derivative Y 1 — _i ^(k) (T _j ) is stored in the storage unit 190 (S202). This process is the same as S102.

Next, the decomposition system coefficient estimation unit 120 performs an estimation process of the decomposition system coefficients β _i and h _ij for each gene, and registers (stores) the data of the estimated decomposition system coefficients β _i and h _{ij in} the storage unit 190. (S203). This process corresponds to S103.

FIG. 6 is a flowchart showing a flow of processing for estimating / registering decomposition system coefficients in the first embodiment.
The processing flow of S203 performed by the decomposition system coefficient estimation unit 120 in the first embodiment will be described below with reference to FIG.

One gene i is selected from unprocessed genes, that is, genes for which the degradation system coefficients β _i and h _ij are estimated (S301).

  Next, with reference to the storage unit 190, it is determined whether there is time course data of the disrupted strain of gene i selected in S301 (S302).

  As a result of the determination in S302, if there is time course data of the disrupted strain of gene i, the time course data is obtained from the storage unit 190, and whether or not the obtained time course data is data with a sufficient number of observation points is observed. Determination is made by “score n_sp ≧ number of genes n + 1” (S303).

In the following processing, a vector c satisfying “b = Ac” is obtained for the vector b and the matrix A by the least square method. In the least square method, the constraint condition number (observation of the disrupted strain) is determined from the unknown number (the number of genes + 1). (Score) needs to be large. That is, when the dimension (number of constraint conditions) of the vector b is smaller than the dimension (unknown number) of the vector c, there is no vector c that satisfies “b = Ac”.
Therefore, the determination “number of observation points n_sp ≧ number of genes n + 1” is performed.

As a result of the determination in S303, when the condition “number of observation points n_sp ≧ number of genes n + 1” is satisfied, the vector b ₁ and the matrix A ₁ are calculated (S304), the vector c ₁ is estimated by the least square method (S305), and the vector The decomposition system coefficients β _i and h _ij for the gene i are calculated from the c ₁ and stored in the storage unit 190 (S306).
The processing of S304 to S306 is the same as S103.

If the result of determination in S303 does not satisfy the condition “number of observation points n_sp ≧ number of genes n + 1”, it is assumed that “the only gene involved is self (gene i)”, and the following vector corresponding to the assumption: b ₃ and matrix A ₃ are calculated (S307).

Here, when a substance decomposes, a certain percentage of the substance is broken (the gene involved is only you), which is a common natural phenomenon, and the rate of decomposition Assuming that “the only gene involved (the gene i) is involved” is that the effect that is controlled (promoted and suppressed) appears to appear as an adverse effect in the term of the production rate that is generated later. Stood up.
As a result, the condition “number of observation points n_sp ≧ number of genes n (1) +1” is satisfied, and the vector c ₃ can be estimated by the following method with the least square method.

Next, a vector c ₃ that satisfies “b ₃ = A ₃ c ₃ ” is estimated by the least square method (S308).

Then, the decomposition system coefficients β _i and h _ij for the gene i are calculated from the vector c ₃ and stored in the storage unit 190 (S309). The relationship between the decomposition system coefficients β _i , h _ij and the vector c ₃ is as follows.

  At this time, the following equation is established based on the assumption that “the gene involved is only me (gene i)”.

h _ij = 0 (i ≠ j) [Formula 2]

If the time course data of the disrupted strain of gene i does not exist as a result of the determination in S302, the degradation system coefficient β _i is calculated so that the generation rate of gene i is positive, and h _ii = 1, h _ij = 0 (i ≠ j) is registered in the storage unit 190 (S310).

As described above, if the time course data of disruption strain of gene i is not present, it is impossible to determine the vector c ₃ by the least squares method.
Therefore, assuming that “the only gene involved is self (gene i)”, Equation 2 is established as in S309.
Further, h _ii is assumed by the following equation.

h _ii = 1 [Formula 3]

Further, the generation rate of gene i is determined to be positive with respect to the time course data of all non-destructive strains. The necessary condition for the generation rate of gene i to be positive is that the conditional expression [Equation 2] of the decomposition system coefficient h _ij is given in advance for all observation points t of all non-destructive strains k. , Β _i can be expressed by the following equation.

Y 1 — _i ^(k) (t)> β _i [Formula 4]

β _i is calculated so as to satisfy Equation 4.

The reason why β _i is determined so that the generation rate of gene i is positive is that the generation rate and decomposition rate of gene i must be positive in the differential equation [Equation 1] of the S-System model. It is.

  Here, the condition that the generation speed is positive can be expressed by the following equation using the differential equation [Equation 1] of the S-System model.

Here, the following equation is established based on the assumption that “h _ij = 0 (i ≠ j)” [Equation 2].

Further, Y _{1 —} i that is a differential coefficient with respect to the logarithmic value of the expression level of gene i can be expressed by the following equation.

  When expanded, the left side can be expressed as:

  In other words, the condition that the generation rate is positive is

Thus, the result β _i satisfies the condition of “Y _{1 — i} ^(k) (t)> β _i ” [Formula 4].

Through the processes in S302 to S310, the decomposition system coefficients β _i and h _ij for the gene i selected in S301 can be estimated and registered in the storage unit 190.

Then, the unprocessed genes for which the degradation system coefficients β _i and h _ij are not estimated and registered (stored) are processed in S301 to S310, and the degradation system coefficients β _i and h of all the genes of interest. _ij is registered in the storage unit 190 (S311).

Next, the generation system coefficient estimation unit 130 performs an estimation process of the generation system coefficients α _i and g _ij for each gene, and registers (stores) the data of the estimated generation system coefficients α _i and g _{ij in} the storage unit 190. . In addition, an evaluation process is performed on the estimated generation system coefficients α _i and g _ij and the evaluation value is registered in the storage unit 190 (S204). In FIG. 5, the processing of S204 to S207 corresponds to S104.

FIG. 7 is a flowchart showing a flow of processing for estimating / registering a generation system coefficient in the first embodiment.
The flow of the process of S204 performed by the generation system coefficient estimation unit 130 in the first embodiment will be described below with reference to FIG.

One gene i is selected from unprocessed genes, that is, genes for which the generation system coefficients α _i and g _ij are estimated (S401).

  Next, the time course data of the non-destructive strain of the gene i selected in S301 is acquired from the storage unit 190, and whether or not the acquired time course data is data of a sufficient number of observation points is “the number of observation points n_sp ≧ the number of genes n + 1”. (S402).

As a result of the determination in S402, when the condition “number of observation points n_sp ≧ number of genes n + 1” is satisfied, the vector b ₂ and the matrix A ₂ are calculated (S403), the vector c ₂ is estimated by the least square method (S404), and the vector Generating system coefficients α _i and g _ij for gene i are calculated from c ₂ and registered in the storage unit 190 (S405).
The processing of S403 to S405 is the same as S104.

  Further, as a result of the determination in S402, when the condition “number of observation points n_sp ≧ number of genes n + 1” is not satisfied, a gene generator set L indicating a set of participating genes with restriction on the number of genes involved in the generation of gene i Is generated. The data of the generative gene set L thus generated is registered in the storage unit 190 (S406).

If the decomposition system coefficients β _i and h _ij can be calculated in the process of S203, the constraint condition expression “b ₂ = A ₂ c ₂ ” applied to the least square method can be calculated. However, when the number of observation points is small, a sufficient number of constraints cannot be obtained to solve the least squares method. For example, this situation occurs when the number of non-destructive strains is small and the number of genes is large.
Therefore, a limit is imposed on the number of genes involved in the generation of gene i to deal with such a situation.

Hereinafter, the maximum number of genes involved that satisfies the condition “number of observation points n_sp ≧ number of genes n + 1” is n _mx . However, this does not include myself. In this case, the generation gene set L can be expressed as follows.

Here, when the total number of genes is n, the number n _cal of the generation gene sets L satisfying this condition can be expressed by the following equation.

For example, when the total number of genes is 5, the target gene number i is 3, and n _mx is 2, the generation gene set L is only L = {3} when m = 1, and L = {1, when m = 2. 3}, {2, 3}, {3, 4}, {3, 5}, m = 3, L = {1, 2, 3}, {1, 3, 4}, {1, 3 , 5}, {2, 3, 4}, {2, 3, 5}, {3, 4, 5}. Further, the number n _cal of the generation gene set L is 11 cases of 1 + 4 + 6.

In the following, processing is performed for each generative gene set L that is generated. If the number of sets _ncal is very large, it is possible to process the generative gene sets L for an allowable number at random to reduce the processing complexity. Good.

  Next, data indicating the generation gene set L generated in S406 is acquired from the storage unit 190, and a target gene set to be processed is selected from the unprocessed generation gene set L (S407).

Next, the following vector b ₄ and matrix A ₄ are calculated for the target gene group selected in S407 (S408).

Next, a vector c ₄ that satisfies “b ₄ = A ₄ c ₄ ” is estimated by the least square method (S409).

Then, the generation system coefficients α _i and g _ij for the gene i are calculated from the vector c ₄ and registered in the storage unit 190 in association with the target gene set (S410). The relationship between the generation system coefficients α _i , g _ij and the vector c ₄ is as follows.

  At this time, the following equation is established for a gene j not included in the target gene set.

Incidentally, the vector _{b 4} are as defined above vector _{b 2.} Further, the matrix _{A 4} and the vector _{c 4} is that the above matrix _{A 2} and the vector _{c 2} is for all genes 1 to n, that for gene _i 1 through i _m a subject genetically element It is.

Next, the differential equation [Equation 1] of the S-System model is numerically integrated based on the generation system coefficients α _i and g _ij calculated in S409 and the decomposition system coefficients β _i and h _ij calculated in the process of S203. Then, the estimated gene expression level of gene i is calculated for the target gene set (S411).

  For example, if the total number of genes is 5 and the target gene is number 3, the differential equation of the S-System model

In order to calculate the value, the expression levels X ₁ , X ₂ , X ₃ , X ₄ , and X ₅ of each gene are required. Here, spline interpolation data is used for expression levels X ₁ , X ₂ , X ₄ , and X ₅ other than the No. 3 target gene. Then, the expression level _{X 3} of the third is the target gene, generating system coefficients calculated in S410 alpha _{i, g ij,} acquired from the catabolic factor beta _{i, h ij} and spline interpolation data calculated in the processing of the S203 Using X ₁ , X ₂ , X ₄ , and X ₅ , the differential equation [Formula 12] is integrated and determined. The result of this integration is the estimated gene expression level.
Here, the spline interpolation data indicates a value obtained from a function expressed by spline interpolation of the time course data in the preprocessing (S202).

Next, the estimated gene expression level of gene i calculated in S411 is compared with spline interpolation values which are values indicated by the spline interpolation data (non-destructive strain time course data), and the generated shape coefficients α _i , g calculated in S410 are compared. An evaluation value for _ij is calculated. Then, the calculated evaluation value is registered in the storage unit 190 in association with the gene i and the target gene set (S412).

An evaluation function for calculating an evaluation value is represented by f _v ⁽ⁱ⁾ (L _j ⁽ⁱ⁾ ).
The evaluation function can be expressed as follows.

  The evaluation function is defined so that the calculated estimated gene expression level becomes larger as the coincidence with the spline interpolation value becomes higher.

Through the processing of S408 to S412, the generation system coefficients α _i and g _ij and the evaluation value for the gene i can be calculated and registered in the storage unit 190 for the target gene group selected in S407.

Then, the processes of S407 to S412 are performed on the unprocessed generation gene set L for which the generation system coefficients α _i and g _ij and the evaluation values are not calculated, and the generation of the gene i is performed for all the generation gene sets L The system coefficients α _i , g _ij and the evaluation value are registered in the storage unit 190 (S413).

Next, the data of the generation system gene set L (generation system gene set, generation system coefficients α _i , g _ij , evaluation value) stored in the storage unit 190 is sorted by evaluation value (S414).

The generation system coefficients α _i , g _ij (and evaluation values) for the gene i selected in S 401 can be calculated and registered in the storage unit 190 by the processing in S 402 to S 414.

Then, the unprocessed genes for which the generation system coefficients α _i , g _ij (and evaluation values) are not calculated are processed in S401 to S414, and the generation system coefficients α _i , g _{ij of} all the genes of interest are obtained. (And the evaluation value) are registered in the storage unit 190 (S415).

By the process of S204, the generation system coefficients α _i and g _ij were calculated for each gene.
At this time, if the number of observation points of the time course data of the non-destructive strain does not satisfy the condition “number of observation points n_sp ≧ number of genes n + 1”, the generation system coefficient α _i , for each set of participating genes (generation system gene set L), g _ij was calculated and evaluated by an evaluation function. Further, the generation system coefficients α _i and g _ij of each generation system gene set L were sorted by the evaluation values calculated by the evaluation function. That is, a plurality of generation system coefficients α _i and g _ij were calculated for each gene (the number of generation system gene sets L).

Hereinafter, by the processes of S205 to S207, the generation system coefficients α _i and g _ij for each gene are coupled and evaluated by integrating the differential equation [Equation 1] of the S-System model, and the optimal generation system for each gene Genes α _i and g _ij are acquired.

First, the first simultaneous evaluation process is performed, and the generation system coefficients α _i and g _ij for each gene are registered (simultaneous evaluation 1) (S205).

FIG. 8 is a flowchart showing the flow of processing for simultaneous evaluation 1 in the first embodiment.
The flow of the process of S205 performed by the simultaneous gene registration unit 141 in the first embodiment will be described below based on FIG.

First, one gene i _{itg (1)} is selected at random from the gene of _interest . The selected gene i _{itg (1)} is registered in the simultaneous gene _Mitg indicating the set of genes for simultaneous evaluation and stored in the storage unit 190 (S501).

The simultaneous gene M _itg at this time is _represented by “M _itg = {i _{itg (1)} }”.

Next, with _respect to the gene i _{itg (1)} selected in S501, the gene _{having the} highest evaluation value among the generation system genes α _i and g _ij corresponding to each generation system gene set L registered in the storage unit 190 in S204 Select. The selected generator coefficient α _i , g _ij (simultaneous generator coefficient) is registered in the storage unit 190 in association with the gene i _{itg (1)} (S502).

  Next, “0” is set to the variable “highest score” for setting the maximum evaluation value by simultaneous evaluation and stored in the storage unit 190 (S503).

Next, a gene i _{itg (j)} (provisional simultaneous gene) is selected at random from genes other than _those registered in the simultaneous gene _Mitg , and is registered in addition to the simultaneous gene _Mitg (S504).

Next, with _respect to the gene i _{itg (j)} selected in S504, the variable “pass flag” indicating that the generation system coefficients α _i and g _ij whose evaluation values by the simultaneous evaluation are equal to or higher than the reference value exists is set to “0 (not set). Pass) ”is set.
Also, the variable “W generating system coefficient” for setting the generating system coefficients α _i and g _ij corresponding to the “highest score” is initialized.
Each variable is stored in the storage unit 190 (S505).

Next, simultaneous evaluation of each gene, each evaluation value is equal to or higher than a reference value, and the generation system coefficient α _i , for the gene i _{itg (j)} , which has the highest evaluation value _all indicating the evaluation of _{all the} simultaneous genes, g _ij is extracted (S506).

FIG. 9 is a flowchart showing a flow of extraction processing of generation system coefficients α _i and g _ij for gene i _{itg (j)} in the first embodiment.
The flow of the process of S506 performed by the simultaneous gene registration unit 141 in the first embodiment will be described below based on FIG.

First, the gene i _{itg (j)} selected in S504 is selected in descending order of the evaluation value for the generation system genes α _i and g _ij corresponding to each generation system gene set L registered in the storage unit 190 in S204. The selected generator coefficients _αi and g _ij (provisional simultaneous generator coefficients) are stored (temporarily registered) in the storage unit 190 in association with the gene i _{itg (j)} (S601).

Next, for each gene registered in the simultaneous gene M _itg , the differential equation [Equation 1] of the S-System model is integrated and integrated using the generation system coefficients α _i and g _ij corresponding to each gene, The estimated gene expression level is calculated (S602).

The calculation process of the estimated gene expression level of each gene is the same as S411.
For example, it is _assumed that the number of genes is 5, the number of genes registered in the simultaneous gene _Mitg is No. 3, and the number of genes newly added in the processing of S504 is No. 5. In this case, spline interpolation data is used for expression levels X ₁ , X ₂ , and X ₄ other than the 3rd and 5th simultaneous genes. The expression amounts X ₃ and X ₅ of the No. 3 and No. 5 simultaneous genes are the generation system coefficients α _i and g _ij registered in S502 and S601, respectively, and the decomposition system coefficients β _i and h calculated in the processing of S203. _The differential equation for gene No. 3 and the differential equation for gene No. 5 are simultaneously integrated using X ₁ , X ₂ , X ₄ obtained from _ij and spline interpolation data, and the estimated gene of gene No. 3 and No. 5 Each expression level is calculated.
In addition, it is _assumed that the genes already registered in the simultaneous gene _Mitg are the 3rd and 5th genes, and the gene newly added in the process of S504 is the 2nd gene. In this case, spline interpolation data is used for expression levels X ₁ and X ₄ other than the simultaneous genes No. 3, No. 5, and No. 2. The expression amounts X ₃ , X ₅ and X ₂ of the simultaneous genes No. 3, No. 5 and No. ₂ are respectively registered generation system coefficients α _i , g _ij , and decomposition system coefficient β _i calculated in the process of S203. , H _ij and X ₁ and X ₄ acquired from the spline interpolation data, the differential equations are integrated and integrated to calculate the estimated gene expression levels of gene No. 3, gene No. 5, and gene No. 2, respectively.

Next, for each gene registered in the simultaneous gene _Mitg , the estimated gene expression level calculated in S602 is evaluated, and each evaluation value is calculated (S603).

  The process of evaluating the estimated gene expression level of each gene is the same as the process (S412).

Next, for each gene registered in the simultaneous gene _Mitg , the evaluation value calculated in S603 is compared with the reference value, and it is determined whether all the evaluation values are equal to or higher than the reference value (pass) (S604).

If the result of the determination in S604 is a pass, all genes registered in the simultaneous gene _Mitg are evaluated, and an evaluation value _all is calculated (S605).

Next, the evaluation value _all calculated in S605 is compared with the highest score stored in the storage unit 190, and it is determined whether the evaluation value _all is larger than the highest score (S606).

As a result of the determination in S606, when the evaluation value _all is larger than the highest score, the highest score is updated with the evaluation value _all . Further, “1 (pass)” is set in the pass flag, and the generation system coefficients α _i and g _ij for the gene i _{itg (j)} provisionally registered in S601 are set in the W generation system coefficient (S607).

Then, the above-described processing of _{S601 to S607} is performed on the generation system gene set L of the _unprocessed gene i _{itg (j)} that has not been simultaneously evaluated for the generation system coefficients α _i and g _ij , and each of the genes i _{itg (j)} is processed. The generation system coefficients α _i and g _ij for the generation system gene set L are simultaneously evaluated (S608).

At this time, the number of generation system gene sets L to be evaluated simultaneously may be limited to the specific number n _test, and only the generation system gene sets L having higher evaluation values may be subjected to simultaneous evaluation. Thereby, the amount of calculation can be suppressed.

By the process of S506 described above based on FIG. 9, the W generation system coefficients α _i and g _ij that passed the simultaneous evaluation for the gene i _{itg (j)} registered in S505 and had the highest evaluation value. Is set. The highest score is set to the highest evaluation value. In the pass flag, presence / absence of passed generation system coefficients α _i and g _ij is set.

Next, the presence / absence of the generation system coefficients α _i and g _ij that have passed with the pass flag is determined (S507).

As a result of the determination in S507, if there is a passed generation system coefficient α _i , g _ij (pass flag = 1), the generation system coefficient α _i , g _ij (simultaneous generation system coefficient) set in the W generation system coefficient is The gene is registered in the storage unit 190 in _association with the gene i _{itg (j)} registered in the simultaneous gene M _itg (S508).

As a result of the determination in S507, when there is no passed generation system coefficient α _i , g _ij (pass flag = 0), the gene i _{itg (j)} is deleted from the simultaneous gene M _itg . Further, the generation system coefficients α _i and g _ij for the temporarily registered gene i _{itg (j)} are deleted (S509).

Next, the gene i _{itg (j)} deleted from the simultaneous gene M _itg is registered in the non-simultaneous gene M _fail and registered in the storage unit 190 (S510).

Then, the processing of S504 to S510 is performed on the remaining genes that have not been evaluated simultaneously, and (maximum) simultaneous genes M _itg and the corresponding generation system coefficients α _i and g _ij are registered in the storage unit 190 (S511). .

Thereafter, by the processes of S206 to S207, the generation system coefficients α _i and g _ij for each gene are updated (S206), and the additional registration of the genes registered in the non-simultaneous gene M _fail to the _coupled gene M _itg (S207) is repeated. Then, the optimum generation system genes α _i and g _ij for each gene are obtained.

In the second simultaneous evaluation process, the generation system coefficients α _i and g _ij for each gene are updated (mountain climbing process) (simultaneous evaluation 2) (S206).

FIG. 10 is a flowchart showing a process flow of simultaneous evaluation 2 in the first embodiment.
A processing flow of S206 performed by the generation system coefficient updating (mountain climbing) unit 142 in the first embodiment will be described below with reference to FIG.

First, “0 (no hill climbing)” is set in the variable “hill climbing flag” indicating that the generation system coefficients α _i and g _ij for any of the genes registered in the simultaneous gene M _itg are updated, and stored in the storage unit 190. (S701).

Next, an unprocessed gene, that is, a gene i _{itg (j)} for updating the generation system coefficients α _i and g _ij is randomly selected from the genes registered in the simultaneous gene M _itg (S702).

Next, “0 (no update)” is set to a variable “update flag” indicating whether or not the generation system coefficients α _i and g _ij to be updated with respect to the gene i _{itg (j)} selected in S702.
Also, the variable “W generating system coefficient” for setting the generating system coefficients α _i and g _ij to be updated is initialized.
Also, it sets the registered product based coefficients alpha _i, g _ij registered product based coefficients alpha _i, g _ij gene variable "R generated based coefficient" for saving the i gene i.
Each variable is stored in the storage unit 190 (S703).

Next, each gene is evaluated in a simultaneous manner, each evaluation value is equal to or higher than a reference value, and the generation system coefficient α corresponding to the gene i _{itg (j)} having the highest evaluation value _all indicating the evaluation of _{all the} simultaneous genes. _i and g _ij are extracted (S704).

The process of S704 is the same as the process of S506.
The “update flag” corresponds to the “pass flag” in S506.
By the processing of S704, the generation system coefficients α _i and g _ij that pass the simultaneous evaluation for the gene i _{itg (j)} selected in S702 and have the highest evaluation value _all are set as the W generation system coefficient. The highest evaluation value _all is set as the highest score. In addition, in the update flag, presence / absence of generation system coefficients α _i and g _{ij to be} updated is set.

Next, it is determined whether or not the generation system coefficients α _i and g _ij to be updated by the update flag are present (S705).

As a result of the determination in S705, when there are generation system coefficients α _i and g _ij to be updated (update flag = 1), the gene i _{itg (j) is} generated with the generation system coefficients α _i and g _ij set in the W generation system coefficient. The generation system coefficients α _i and g _ij registered in (1) are updated. Further, “1 (with mountain climbing)” is set in the mountain climbing flag (S706).

If the generation system coefficients α _i and g _{ij to} be updated are not found as a result of the determination in S705 (update flag = 0), the generation system coefficients α _i and g _ij saved in the R generation system coefficient in S703 are used as the simultaneous gene M _itg . Corresponding to the registered gene i _{itg (j)} , it is registered in the storage unit 190 and returned to the original generation system coefficients α _i , g _ij (S707).

Then, the processing of S702 to S707 is performed on the remaining genes that have not been processed, and the updated generation system coefficients α _i and g _ij are registered in the storage unit 190 (S708).

The simultaneous registration in the simultaneous evaluation 2 (mountain climbing) (S206) described above and the simultaneous evaluation in the simultaneous evaluation 1 (S205) are different because the genes registered in the simultaneous gene _Mitg are different. The evaluation results are also different, and there are cases where updating can be performed with more suitable generation system coefficients α _i and g _ij .
Similarly, in some cases, each gene registered in the non-simultaneous gene M _fail in the simultaneous evaluation 1 (S205) can be additionally registered in the simultaneous gene M _itg .

Therefore, by the third simultaneous evaluation process, an additional registration process to the simultaneous gene M _itg is performed for each gene registered in the non-simultaneous gene M _fail (simultaneous evaluation 3) (S207).

FIG. 11 is a flowchart showing a flow of processing of simultaneous evaluation 3 in the first embodiment.
The process flow of S207 performed by the simultaneous gene addition registration unit 143 in the first embodiment will be described below with reference to FIG.

First, with reference to the storage unit 190, it is determined whether or not there is a gene registered in the non-simultaneous gene M _fail (S801).

As a result of the determination in S801, if there is a gene registered in the non-simultaneous gene M _fail , “0 (no additional registration)” is set in the “simultaneous increase flag” indicating the presence / absence of the gene additionally registered in the simultaneous gene M _itg , It memorize | stores in the memory | storage part 190 (S802).

Next, the gene unprocessed from the gene registered in the non-simultaneous gene _{M fail,} that is, one randomly selects a gene _{i ITG (j)} for additionally registering the simultaneous gene _{M itg} (S803).

Hereinafter, the processing of S804 to S811 is the same as the processing of S505 to S511, except for S808.
In S808, “1 (additional registration exists)” is set in the simultaneous increase flag.

Next, referring to the simultaneous increase flag stored in the storage unit 190, it is determined whether there is a gene additionally registered from the non-simultaneous gene M _fail to the simultaneous gene M _itg (S812).

As a result of the determination in S812, when there is a gene additionally registered from the non-simultaneous gene M _fail to the simultaneous gene M _itg (simultaneous increase flag = 1), the simultaneous evaluation 2 (mountain climbing) process is performed again (S206).

If the genes registered in the simultaneous gene M _itg are different, the results of the simultaneous evaluation are also different. Therefore, there are cases where the generative coefficients α _i and g _ij can be updated more appropriately.
Thereafter, simultaneous evaluation 2 (S206) and simultaneous evaluation 3 (S207) are repeated until the generation system coefficients α _i and g _ij are not improved.

As a result of the determination in S812, when there is no gene additionally registered from the non-simultaneous gene M _fail to the simultaneous gene M _itg (simultaneous increase flag = 0), or, as a result of the determination in S801, a gene registered in the non-simultaneous gene M _fail If not, the hill-climbing flag stored in the storage unit 190 is referred to and it is determined whether the generation system coefficients α _i and g _ij updated in the simultaneous evaluation 2 (S206) are present (whether hill-climbed) (S813).

As a result of the determination in S813, when the updated generation system coefficients α _i and g _ij are present (hill climbing), simultaneous evaluation 2 (hill climbing) processing is performed again (S206).

If the generation system coefficients α _i and g _ij for the genes registered in the simultaneous gene Mitg are different, the result of the simultaneous evaluation is also different, so that it may be possible to update the generation system coefficients α _i and g _ij to be more suitable.
Thereafter, simultaneous evaluation 2 (S206) and simultaneous evaluation 3 (S207) are repeated until the generation system coefficients α _i and g _ij are not improved.

Even when the time course data of the disrupted strain or the non-destroyed strain is insufficient by the processing of S201 to S207 described above, each coefficient (generation system) more suitable for the differential equation [Formula 1] of the S-System model. Coefficient, decomposition system coefficient). In addition, it is possible to reduce the amount of calculation in each coefficient acquisition process.
In addition, since genes are selected in a random order in each coefficient acquisition process, different results may be obtained when each coefficient acquisition process is performed again. Thus, by performing each coefficient acquisition process a plurality of times, a more suitable coefficient can be acquired based on each result.
Then, a differential equation [Equation 1] of the S-System model can be generated from each acquired coefficient, and the gene expression interaction can be estimated.

Next, the generation system related gene estimation process (single evaluation) performed in S208 will be described.
FIG. 12 is a flowchart showing the flow of the single evaluation process in the first embodiment.
In the single evaluation process, for each gene, a gene (gene controlling gene) involved in gene generation (gene control) is estimated.
The single evaluation process (S208) performed by the single evaluation unit 150 in the first embodiment will be described below with reference to FIG.

Here, the generation system coefficients α _i and g _ij can be calculated for a set of genes involved in generation (generation system gene set L) for each gene by the processing of S201 to S204 described above, and The evaluation value can be calculated by the evaluation function.

Therefore, this evaluation value is added for each gene involved, and a score f _sgl ⁽ⁱ⁾ (i _x ) is calculated for each gene (S901).

The score f _sgl ⁽ⁱ⁾ (i _x ) can be expressed as:

The score f _sgl ⁽ⁱ⁾ (i _x ) is the sum of the evaluation values of the generation system coefficients for the gene i when the gene i _x is included in the generation system gene set L, and the score f _sgl ⁽ⁱ⁾ ( It is considered that the larger i _x ) is, the higher the possibility that the gene i _x controls the gene i. However, since the gene i is included in any gene series L, the score f _sgl ⁽ⁱ⁾ (i) of the gene i is not evaluated.

Then, the score _f ^sgl for gene _{i x} other than gene _{i (i) (i x)} [i x: 1, ···, i-1, i + 1, ···, n] to be likened to statistics Z score ^{_{Z sgl (i) (i x}} ) [i x: 1, ···, i-1, i + 1, ···, n] is calculated (S902).

When the Z-score value is greater than the reference value, and the gene i _x are controlled gene i.

For example, if the number of genes is 5 and the gene i is number 3, the Z score Z _sgl ⁽ⁱ⁾ (i _x ) can be expressed by the following equation.

Next, the gene involved in the generation system is estimated by statistical processing based on the Z score Z _sgl ⁽ⁱ⁾ (i _x ) calculated in S902 (S903).

  By the single evaluation process (S208) described above, it is possible to estimate a gene for a gene for each gene with a small amount of calculation.

The flowchart which shows the flow of a process of the gene expression interaction estimation method in Reference Example 1. FIG. 3 shows an appearance of a gene expression interaction estimation apparatus 100 according to the first embodiment. FIG. 3 is a hardware configuration diagram of the gene expression interaction estimation device 100 according to the first embodiment. 1 is a configuration diagram of a gene expression interaction estimation apparatus 100 according to Embodiment 1. FIG. 5 is a flowchart showing a processing flow of a gene expression interaction estimation method according to Embodiment 1. 5 is a flowchart showing a flow of processing for estimation / registration of decomposition system coefficients in the first embodiment. 6 is a flowchart showing a flow of processing for estimating / registering a generation system coefficient in the first embodiment. 5 is a flowchart showing a flow of processing for simultaneous evaluation 1 in the first embodiment. 9 is a flowchart showing a flow of extraction processing for generating system coefficients α _i and g _ij for gene i _{itg (j)} in the first embodiment. 4 is a flowchart showing a flow of simultaneous evaluation 2 processing in the first embodiment. 5 is a flowchart showing a flow of processing for simultaneous evaluation 3 in the first embodiment. 5 is a flowchart showing the flow of a single evaluation process in the first embodiment.

Explanation of symbols

  DESCRIPTION OF SYMBOLS 100 Gene expression interaction estimation apparatus, 110 Pre-processing part, 120 Decomposition system coefficient estimation part, 130 Generation system coefficient estimation part, 140 Simultaneous evaluation part, 141 Simultaneous gene registration part, 142 Generation system coefficient update (mountain climbing) part, 143 Simultaneous Gene addition registration unit, 150 single evaluation unit, 190 storage unit, 901 CRT display device, 902 K / B, 903 mouse, 904 FDD, 905 CDD, 906 printer device, 907 scanner device, 910 system unit, 911 CPU, 912 bus 913 ROM, 914 RAM, 915 communication board, 920 magnetic disk device, 921 OS, 922 window system, 923 program group, 924 file group, 931 telephone, 932 FAX machine, 940 Internet, 941 web server, 94 2 LAN.

Claims (15)

This is an expression rate expression of an S-System model for obtaining an expression rate of a gene i (1 ≦ i ≦ n) included in n genes (n is an integer of 1 or more), and includes a decomposition system coefficient of the gene i. -In an expression rate expression generation device for calculating a decomposition system coefficient of gene i included in an expression rate expression of a System model,
A preprocessing unit for inputting time course data of the disrupted strain indicating the expression level of the gene i observed at each time from an external device;
The number of expression levels of gene i indicated in the time course data of the disrupted strain input by the pre-processing unit and n + 1 are compared using a CPU (Central Processing Unit), and the gene indicated in the time course data of the disrupted strain It is determined whether or not the number of expression levels of i is n + 1 or more using a CPU, and when the number of expression levels of gene i indicated in the time course data of the disrupted strain is not n + 1 or more, the decomposition of gene i Decomposition system coefficient that uses a CPU to calculate the degradation system coefficient of gene i by the least squares method based on the expression level of gene i shown in the time course data of the disrupted strain under the condition that only gene i is involved Estimator and
The expression rate type | formula production | generation apparatus characterized by the above-mentioned. The expression rate equation of the S-System model includes a generation coefficient of gene i,
The pre-processing unit inputs time course data of a non-destructive strain indicating the expression level of the gene i observed at each time,
The expression rate expression generator further includes:
The number of expression levels of gene i indicated in the time course data of the non-destructive strain input by the pre-processing unit and n + 1 are compared using a CPU, and the expression of gene i indicated in the time course data of the non-destructive strain is expressed. It is determined whether or not the number of quantities is n + 1 or more by using a CPU, and when the number of expression levels of gene i indicated in the time course data of the non-destructive strain is not n + 1 or more, the time course of the non-destructive strain In the time course data of the non-destructive strain, the number of genes i expressed in the data is m + 1 (m is an integer of 1 or more) and m is the number of genes involved in the generation of gene i. A generation system coefficient estimation unit is provided that calculates a generation system coefficient of gene i using a CPU based on the expression level of gene i indicated by the least square method.
The expression rate type | formula production | generation apparatus of Claim 1 characterized by the above-mentioned. The expression rate equation of the S-System model is expressed by a decomposition system coefficient β of gene i. _i Interaction coefficient h of gene j (1 ≦ j ≦ n) with respect to gene i _ij And an expression represented by the expression (1) including:
The preprocessing unit calculates a differential coefficient of logarithmic value corresponding to the expression level of the gene i indicated in the time course data of the input disrupted strain,
When the number of expression levels of the gene i indicated in the time course data of the disrupted strain is not n + 1 or more, the degradation system coefficient estimation unit determines the expression level of the gene i indicated in the time course data of the disrupted strain and the preprocessing Vector b expressed by equation (2) using the logarithmic value derivative calculated by the unit ₃ And the matrix A represented by the equation (3) ₃ And "b ₃ = A ₃ c ₃ C satisfying " ₃ Is estimated by the least square method, and the estimated vector c ₃ And the decomposition coefficient β of the gene i using the equations (4) and (5) _i Coefficient h of the degradation system of gene j with respect to gene i _ij And calculate
The expression rate type | formula production | generation apparatus of Claim 2 characterized by the above-mentioned.

When the number of expression levels of the gene i indicated in the time course data of the non-destructive strain is not n + 1 or more, the generation system coefficient estimation unit determines the number of expression levels of the gene i indicated in the time course data of the non-destructive strain A plurality of generative gene sets L are generated based on m with m + 1 or more and formula (6), and a vector b represented by formula (7) for each generative gene set L ₄ And a matrix A represented by equation (8) ₄ And "b ₄ = A ₄ c ₄ C satisfying " ₄ Is estimated by the least square method, and the estimated vector c ₄ And the expression coefficient α of the gene i included in the expression rate expression (1) of the S-System model based on Expression (9) and Expression (10) _i Interaction coefficient g of gene j with respect to gene i _ij And calculate
The expression rate type | formula production | generation apparatus of Claim 3 characterized by the above-mentioned.

The generating system coefficient estimation unit calculates the generating system coefficient α calculated for each generating system gene set L. _i And generation system interaction coefficient g _ij And the decomposition system coefficient β calculated by the decomposition system coefficient estimation unit _i Coefficient h _ij Based on the above, the expression rate expression (1) of the S-System model is numerically integrated to calculate the estimated gene expression level of the gene i for each generation gene set L using the CPU, and the estimated gene expression of the calculated gene i is calculated. The production system coefficient α for each production system gene set L based on the expression level and the expression level indicated in the non-destructive strain time course data and the equation (11) _i And generation system interaction coefficient g _ij The evaluation value for and is calculated using the CPU
The expression rate expression generation device according to claim 4.

The expression rate expression generator further includes:
The production system coefficient α calculated for each production system gene set L by the production system coefficient estimation unit _i And generation system interaction coefficient g _ij The gene i included in the gene family L based on the evaluation value for and the expression (12) _x Every gene i _x Is calculated using the CPU, and the gene i _x Gene i calculated for each _x And the gene i based on the score of (13) _x Every gene i _x The Z score of is calculated using a CPU and gene i _x Gene i calculated for each _x A single evaluation unit that uses a CPU to estimate a gene involved in the generation system involved in the generation of gene i based on the Z score of
The expression rate expression generation device according to claim 5.

The generating system coefficient estimation unit generates a generating system gene set L for each of the remaining genes excluding the gene i, and the generating system coefficient α of the generated generating system gene set L _i And generation system interaction coefficient g _ij And
The generating system coefficient estimation unit generates the generating system coefficient α for each generating system gene set L of each of the remaining genes excluding the gene i. _i And generation system interaction coefficient g _ij And the evaluation value for
The expression rate expression generator further includes:
The generative coefficient α of each generative gene set L of gene i calculated by the generative coefficient estimation unit _i And generation system interaction coefficient g _ij The generation coefficient α having the highest evaluation value calculated by the generation coefficient estimation unit _i And generation system interaction coefficient g _ij Are registered in the storage device as simultaneous generation system coefficients of the gene i used in the expression rate equation (1) of the S-System model, and genes other than the gene i from the n genes are randomly selected as temporary simultaneous genes one by one Select the selected temporary simultaneous gene as the simultaneous gene M _itg And the generation system coefficient α of each generation system gene set L of the temporary simultaneous genes calculated by the generation system coefficient estimation unit _i And generation system interaction coefficient g _ij The generation coefficient α of the temporary simultaneous gene in descending order of the evaluation value calculated by the generation coefficient estimation unit _i And generation system interaction coefficient g _ij Is selected as the simultaneous generation coefficient of the temporary simultaneous gene, the expression rate expression (1) of the gene i is generated using the simultaneous generation system coefficient of the gene i, and the temporary simultaneous generation is performed using the simultaneous generation system coefficient of the temporary simultaneous gene. Gene expression rate expression (1) is generated, and expression rate expression (1) for gene i and expression rate expression (1) for temporary simultaneous genes are integrated and integrated, and the estimated gene expression level of gene i and temporary simultaneous expression An estimated gene expression level of the gene is calculated, an evaluation value for the simultaneous generation system coefficient of the gene i is calculated based on the estimated gene expression level of the gene i and the formula (11), and the estimated gene expression level of the temporary simultaneous gene Based on the equation (11), an evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene is calculated, and an evaluation value for the simultaneous generation system coefficient of the gene i and an evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene are set to a predetermined criterion. Compared to the value, simultaneous generation of gene i When the evaluation value for the system coefficient and the evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene are greater than or equal to the reference value, the evaluation value for the simultaneous generation system coefficient of the gene i and the evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene Total value is evaluated value _all And the evaluation value for the simultaneous generation system coefficient of the gene i and the evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene are equal to or greater than the reference value. If there are no simultaneous generative coefficients for the gene set L of the temporary gene, the temporary gene is assigned to the gene M _itg The evaluation value for the simultaneous generation system coefficient of the gene i and the evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene are equal to or greater than the reference value. If there is a simultaneous generation system coefficient of the generation system gene set L, the evaluation value _all A simultaneous evaluation unit that registers the highest simultaneous generation coefficient in the storage device as the simultaneous generation coefficient of the temporary simultaneous gene
The expression rate expression generation device according to claim 5 or 6, characterized in that The simultaneous evaluation unit randomly selects a new temporary simultaneous gene, and selects the selected new temporary simultaneous gene as the simultaneous gene M. _itg In addition, a new temporary simultaneous gene coefficient is selected in descending order of the evaluation value calculated by the generating system coefficient estimation unit, and the expression rate of the new temporary simultaneous gene is selected using the selected simultaneous generation system coefficient. Generating the expression (1), integrating the expression rate expression (1) of the gene i, the expression speed expression (1) of the temporary simultaneous gene, and the expression speed expression (1) of the new temporary simultaneous gene, Calculate the estimated gene expression level of i, the estimated gene expression level of the temporary simultaneous gene, and the estimated gene expression level of the new temporary simultaneous gene, the evaluation value for the simultaneous generation system coefficient of gene i and the simultaneous generation system coefficient of the temporary simultaneous gene And the evaluation value for the simultaneous generation system coefficient of the new temporary simultaneous gene is calculated, and when each calculated evaluation value is equal to or higher than the reference value, the total value of the calculated evaluation values is calculated as the evaluation value. _all If there is no simultaneous generation coefficient of a new temporary simultaneous gene whose calculated evaluation value is equal to or greater than the reference value, the temporary simultaneous gene is determined as the simultaneous gene M. _itg If there is a simultaneous generation coefficient of a new temporary simultaneous gene whose calculated evaluation value exceeds the reference value, the evaluation value _all Register the highest simultaneous generative coefficient to the storage device as a new temporary gene gene
The expression rate type | formula production | generation apparatus of Claim 7 characterized by the above-mentioned. The simultaneous evaluation unit is a simultaneous gene M _itg Gene is selected one by one from the genes registered as, and the generation coefficient α of each generation gene set L of the selected gene _i And generation system interaction coefficient g _ij From the above, the simultaneous generation system coefficients of the genes selected in descending order of the evaluation value calculated by the generation system coefficient estimation unit are selected, and the expression rate expression (1) of the selected gene is generated using the selected simultaneous generation system coefficients Expression rate formula (1) of gene i, expression rate formula (1) of the selected gene, and simultaneous genes M other than the selected gene _itg The expression rate expression (1) is integrated and integrated, and the estimated gene expression level of gene i, the estimated gene expression level of the selected gene, and the simultaneous gene M other than the selected gene _itg The estimated gene expression level of the gene i is calculated, the evaluation value for the simultaneous generation system coefficient of the gene i, the evaluation value for the simultaneous generation system coefficient of the selected gene, and the simultaneous gene M other than the selected gene _itg The evaluation value for the simultaneous generation system coefficient is calculated, and when each calculated evaluation value is greater than or equal to the reference value, the total value of the calculated evaluation values is the evaluation value. _all And calculated evaluation value _all Evaluation value of registered simultaneous system coefficient of gene selected by _all If the value is larger, update the registered simultaneous generator coefficient of the selected gene with the selected simultaneous generator coefficient.
The expression rate expression generation device according to claim 8. The simultaneous evaluation unit includes the gene i and the simultaneous gene M among n genes. _itg Non-cognate gene M excluding _fail The temporary simultaneous genes are selected one by one at random, and the selected temporary simultaneous genes are selected as the simultaneous genes M. _itg And the generation coefficient α of each generation gene set L of the provisional simultaneous gene _i And generation system interaction coefficient g _ij From the above, the simultaneous generation system coefficients of the temporary simultaneous gene are selected in descending order of the evaluation value calculated by the generation system coefficient estimation unit, and the expression rate expression (1) of the temporary simultaneous gene is generated using the selected simultaneous generation system coefficient And the expression rate expression (1) of gene i and the coupled gene M _itg The expression rate equation (1) of the gene and the expression rate equation (1) of the temporary simultaneous gene are integrated and integrated, and the estimated gene expression level of the gene i and the simultaneous gene M _itg The estimated gene expression level of the gene and the estimated gene expression level of the temporary simultaneous gene are calculated, and the evaluation value for the simultaneous generation system coefficient of the gene i and the simultaneous gene M _itg Calculate the evaluation value for the simultaneous generation system coefficient of and the evaluation value for the simultaneous generation system coefficient of the temporary simultaneous gene. If each evaluation value is equal to or greater than the reference value, the total value of the calculated evaluation values is the evaluation value. _all If there is no simultaneous generation coefficient of the temporary simultaneous gene for which each calculated evaluation value is equal to or greater than the reference value, the temporary simultaneous gene is determined as the simultaneous gene M. _itg If there is a simultaneous generation system coefficient of a temporary simultaneous gene that is deleted from _all Register the highest simultaneous generative coefficient to the storage device as the simultaneous generative coefficient of the temporary gene
The expression rate type | formula production | generation apparatus of Claim 9 characterized by the above-mentioned. This is an expression rate expression of an S-System model for obtaining an expression rate of a gene i (1 ≦ i ≦ n) included in n genes (n is an integer of 1 or more), and includes a decomposition system coefficient of the gene i. In the expression rate expression generation method of the expression rate expression generation device that calculates the decomposition system coefficient of the gene i included in the expression rate expression of the system model,
The preprocessing unit inputs time course data of the disrupted strain indicating the expression level of the gene i observed at each time from an external device,
The degradation system coefficient estimation unit compares the number of expression levels of gene i indicated in the time course data of the disrupted strain input by the preprocessing unit with n + 1 using a CPU (Central Processing Unit), and It is determined whether or not the number of gene i expressed in the time course data is n + 1 or more using the CPU, and the number of gene i expressed in the time course data of the disrupted strain is not n + 1 or more. In this case, under the condition that only the gene i is involved in the degradation of the gene i, the degradation coefficient of the gene i is calculated by the least square method based on the expression level of the gene i indicated in the time course data of the disrupted strain. Calculate using
An expression rate expression generation method characterized by the above. The expression rate equation of the S-System model includes a generation coefficient of gene i,
The pre-processing unit inputs time course data of a non-destructive strain indicating the expression level of the gene i observed at each time,
The expression rate expression generator further includes:
The generation system coefficient estimation unit compares the number of expression levels of the gene i indicated in the time course data of the non-destructive strain input by the preprocessing unit with n + 1 using a CPU, and the time course data of the non-destructive strain Whether the number of expression levels of gene i shown in (b) is n + 1 or more using a CPU, and if the number of expression levels of gene i shown in the time course data of the non-destructive strain is not n + 1 or more, On the condition that m is the number of genes involved in the generation of gene i, where the number of expression levels of gene i indicated in the time course data of the non-destructive strain is m + 1 (m is an integer of 1 or more) Based on the expression level of the gene i shown in the time course data of the disrupted strain, the gene i generation system coefficient is calculated by the CPU using the least square method.
The expression rate expression generation method according to claim 11. The expression rate equation of the S-System model is expressed by a decomposition system coefficient β of gene i. _i Interaction coefficient h of gene j (1 ≦ j ≦ n) with respect to gene i _ij And an expression represented by the expression (1) including:
The preprocessing unit calculates a differential coefficient of logarithmic value corresponding to the expression level of the gene i indicated in the time course data of the input disrupted strain,
When the number of expression levels of the gene i indicated in the time course data of the disrupted strain is not n + 1 or more, the degradation system coefficient estimation unit determines the expression level of the gene i indicated in the time course data of the disrupted strain and the preprocessing Vector b expressed by equation (2) using the logarithmic value derivative calculated by the unit ₃ And the matrix A represented by the equation (3) ₃ And "b ₃ = A ₃ c ₃ C satisfying " ₃ Is estimated by the least square method, and the estimated vector c ₃ And the decomposition coefficient β of the gene i using the equations (4) and (5) _i Coefficient h of the degradation system of gene j with respect to gene i _ij And calculate
The expression rate expression generation method according to claim 12.

When the number of expression levels of the gene i indicated in the time course data of the non-destructive strain is not n + 1 or more, the generation system coefficient estimation unit determines the number of expression levels of the gene i indicated in the time course data of the non-destructive strain A plurality of generative gene sets L are generated based on m with m + 1 or more and formula (6), and a vector b represented by formula (7) for each generative gene set L ₄ And a matrix A represented by equation (8) ₄ And "b ₄ = A ₄ c ₄ C satisfying " ₄ Is estimated by the least square method, and the estimated vector c ₄ And the expression coefficient α of the gene i included in the expression rate expression (1) of the S-System model based on Expression (9) and Expression (10) _i Interaction coefficient g of gene j with respect to gene i _ij And calculate
The expression rate expression generation method according to claim 13.

An expression rate expression generation program for causing a computer to execute the expression rate expression generation method according to claim 11.

Images