JP2014035653A - Population size estimation method, population estimation method, gene expression analysis method, data analysis method, sensitivity measurement method, program, recording medium and population size estimation system - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a population size estimation method for precisely estimating the size of population using measurement data obtained by using a measurement method having detection limit.SOLUTION: The total number of data is estimated on measurement data obtained using a measurement method having detection limit is made on a computer. The measured data is sorted in an order from the largest/smallest data. Subsequently, a range of candidate values of the total number of data is set based on the measured data. Then, scatter diagrams of probability distribution of a normalized measured data and an estimated distribution of the measured data are drawn with respect to all candidate values within the range. These scatter diagrams are evaluated on the linearity using a correlation like Pearson's R-value or a function on a linear relationship to estimate an optimum candidate value as the population size.

Description

  The present invention relates to a population size estimation method, a parameter estimation method, a gene expression level analysis method, a data analysis method, a sensitivity measurement method, a program, a recording medium, and a population size estimation apparatus, and more particularly to detection sensitivity. The present invention relates to a population size estimation method, a population estimation method, a gene expression level analysis method, a data analysis method, a sensitivity measurement method, a program, a recording medium, and a population size estimation apparatus when using data with limitations.

The inherited biological information is recorded in the genome in the form of a DNA base sequence. DNA is transcribed into RNA in the cell, and the RNA accumulates in the cell.
These RNAs retain information on specific proteins, but RNAs that are frequently translated maintain high intracellular concentrations (expression levels). That is, the concentration of these RNAs is a stage where genetic information changes to a physical form for the first time. Content related to the type and expression level of RNA contained in a cell is referred to as a transcriptome. The expression level of each RNA is an important clue for knowing the function and state of the cell.

Since there are tens of thousands of types of RNA, it is not easy to quantify after identifying these types. However, the DNA microarray method using a DNA chip appeared in the 1990s, and its comprehensive analysis became possible. This is a method of obtaining relative values having different sensitivities for each gene as measured values from hybridization images and the like.
There are two methods for measuring transcriptome without using a DNA chip: a method using quantitative PCR and a method using a DNA sequencer.
Furthermore, a method for identifying a substance by mass-analyzing an organic compound, protein, metabolite or the like, which is an intracellular metabolite, while separating and quantifying it by a method such as chromatography has been put into practical use. Examples of methods for identifying these substances include LC-MS (Liquid Chromatography Mass Spectrometer).

Here, when analyzing measurement data of a relative value composed of a plurality of samples, the first thing to do is standardization of the measurement data. Samples used for measurement are often arranged in some way in advance, but numerical values obtained in the aligned state are not absolute values. This is because the sensitivity varies from sample to sample on the principle of measurement. Because of this sensitivity difference, data cannot be directly compared between measurements. Therefore, it is necessary to standardize the measurement data separately. This standardization allows each sample to be compared.
In order to standardize measurement data, a “parameter” representing the distribution of the population of data is estimated. Here, all the substances to be measured contained in a certain sample are considered to be a population. For example, even if the data of a “strongly” sample that was lost from a specific intensity range is known, it is possible to determine how strongly the sample was measured if the parameter representing the nature of the population is known. Can be estimated. Therefore, comparison between measurements is possible by using the parameter.
For example, in the case of a normal distribution model, there is a parameter μ at the center of the distribution estimated by the average and a parameter σ of the width of the distribution estimated by the standard deviation. These parameter estimates are important to provide the basis for any subsequent analysis.

With reference to the example of FIG. 9, the relationship between the standardization of measurement data and the population will be described in more detail.
For example, as shown in FIG. 9A, when the original population is data according to a normal distribution with mean μ = 0 and standard deviation σ = 1, the histogram of the data is as shown in FIG. 9A. become.
Conventionally, there are many measurement methods with a detection limit such that only data that can be detected is obtained as data due to a sensitivity limit or the like. For example, those that record data over time, such as electrophoresis, chromatography, and mass spectrometry, always have detection limits. In addition, as with RNA-seq using a next-generation sequencer, data that has not been read cannot be obtained even in the measurement using the Bernoulli trial that data is basically “read / not read”. There is a detection limit. Furthermore, hybridization using fluorescence or an electrical detection method has a detection limit related to the sensitivity limit of the instrument.
FIG. 9B is an example of a histogram of measured values of data having a detection limit. In this way, a part of the distribution smaller than the detection limit is often specifically lost. Conversely, a large signal that exceeds the detection limit may be ignored, or the data may be lost due to failure of the measuring instrument.

Conventionally, when a measurement method having a detection limit is used, the size of the population of measurement data is unknown, and thus distribution characteristics and parameters are often not obtained.
Here, the “population size” is a value indicating the total number of types of objects to be measured. Conventionally, assuming that the data is a random sample from the population, for example, assuming a normal distribution, the number of parameters is estimated, and the size of the generating means has been standardized.
However, for example, when the average is calculated from the data of FIG. 9B, the standard deviation is 0.3 and the standard deviation is about 0.8. This is far from the average 0 and standard deviation 1 of the original population. These differences are caused by the assumption that data is lost from a certain intensity range and is a random sampling.

In addition, when a part of the distribution is specifically lost in this way, the parameter cannot be calculated correctly even by using a so-called robust calculation method. This is because the cause of the disturbance is not a small number of outliers. For example, if the median, which is a commonly used robust μ estimation method, is calculated with the data in FIG. 9B, it is 0.22, and MAD (median absolute deviation), which is also a robust σ estimation method, is 0. 84.
In addition, in practice, the detection limit depends on the sample and measurement. This is because the total amount of signal is not constant for each individual sample being measured, even if the measurement equipment is the same.
That is, the proportion of disappearing areas that cannot be detected below the sensitivity limit, as shown in gray (filled) in the examples of FIGS. 9C, 9D, and 9E, is not constant. Therefore, in FIG. 9C, the average is 0.15 and the standard deviation is 0.88, and in FIG. 9D, the average is 0.30 and the standard deviation is 0.80. Then the average is 0.51 and the standard deviation is 0.70.

This variation in the ratio of the disappearing region becomes remarkable particularly when the population has a lognormal distribution. For example, due to the thermodynamic properties of cells, mRNA is log-normally distributed. Because of this property, many quantitative omics data are expected to have a similar distribution.
In the lognormal distribution, the lowest data appears most frequently. Therefore, a slight difference in detection sensitivity greatly changes the number of data obtained. Naturally, the estimate of the distribution parameter is greatly influenced by the number of data lost.

Here, as a standardization method of conventional DNA microarray data, with reference to Patent Document 1, in a data set obtained from one DNA microarray, a step of rearranging each data in order of size and determining the rank of the data; DNA microarray data comprising: inputting a rank of certain data to an inverse function of a standard normal cumulative distribution function, outputting a predicted value after standardization of the data; and replacing the data with the predicted value (Hereinafter referred to as prior art 1).
According to the processing method of the prior art 1, it is possible to provide a standardization method that avoids the distortion from linearity in the standardization of the conventional parametric method and the necessity of the iterative calculation necessary for parameter estimation.

International Publication No. 2008/056693

Any conventional measurement method with a detection limit has a limit in sensitivity, and the limit differs for each measurement. In addition, the expression level data related to transcriptome and the like basically has a lognormal distribution. For this reason, the difference in the number of measurement data (population size) that is influenced by a slight difference in measurement limit tends to increase.
However, the prior art 1 did not assume an accurate estimation of the size of the population of measurement data when data is lost due to a detection limit.

  This invention is made | formed in view of such a condition, and makes it a subject to eliminate the above-mentioned subject.

The population size estimation method of the present invention sorts the measurement data in a population size estimation method that causes a computer to perform estimation of the total number of measurement data obtained by a measurement method with a detection limit. Prepare a plurality of candidate values of the total number of data within the range calculated from the measurement data, prepare a set of theoretical values related to the population distribution of the measurement data for each candidate value, and sort The measured data and the set of theoretical values are evaluated by a function related to a correlation or linear relationship, and the candidate value that is optimally evaluated by the function is estimated as a population size, To do.
The population size estimation method according to the present invention is characterized in that the function evaluates a linearity or a distance between the value of the measurement data in the same order from the top and the value of the set of theoretical values.
The population size estimation method of the present invention uses the Pearson product moment correlation coefficient as the function, and the optimal candidate value is at the end of the range, or the prepared candidate value has low accuracy Is characterized in that the range of the candidate values is adjusted, and when the sufficient accuracy is obtained, the candidate values that are optimally evaluated are determined.
The parameter estimation method of the present invention uses the set of the theoretical values of the candidate values that are optimally evaluated and the sorted measurement data, and the values of the measurement data in the same order from the top and the theory It is characterized in that a set of values is compared, a first-order approximation intercept and slope are calculated, and a parameter indicating the characteristics of the population distribution is calculated.
The gene expression level analysis method of the present invention uses the gene expression level quantification data as measurement data, and standardizes the gene expression level based on the total number of data estimated by the population size estimation method. And
The gene expression level analysis method of the present invention is characterized in that the quantification data of the gene expression level is a transcript expression level, a cDNA reading count, and derivative data.
The proteome or metabolome data analysis method of the present invention uses the protein amount or substance amount quantification data as measurement data, and the protein amount or substance amount based on the total number of data estimated by the population size estimation method. Is standardized.
The sensitivity measurement method of the DNA microarray of the present invention uses RNAseq and quantification data of the gene expression level of the DNA microarray as measurement data, standardized by the total number of data estimated by the population size estimation method, The sensitivity of each spot of the DNA microarray is obtained by correcting the standardized gene expression level of the RNAseq.
The program of the present invention is a program for causing a computer to execute a population size estimation method for estimating the total number of measurement data obtained by a measurement method having a detection limit, and causing the computer to sort the measurement data. , In a range calculated by the measurement data, to prepare a plurality of candidate values of the total number of data, for each candidate value, to prepare a set of theoretical values related to the population distribution of the measurement data, A population size that evaluates a match between the sorted measurement data and the set of theoretical values by a function related to a correlation or linear relationship, and estimates the candidate value that is optimally evaluated by the function as the size of the population It is characterized in that the estimation method is executed.
The recording medium of the present invention is a computer-readable recording medium that records a program that causes a computer to execute a population size estimation method for estimating the total number of measurement data obtained by a measurement method with a detection limit. Causing the computer to sort the measurement data and prepare a plurality of candidate values of the total number of data within a range calculated by the measurement data; for each candidate value, a distribution of the population of the measurement data A set of theoretical values according to the above is prepared, and the coincidence between the sorted measurement data and the set of theoretical values is evaluated by a function related to a correlation or linear relationship, and the candidate value that is optimally evaluated by the function is determined as a mother A program for executing a population size estimation method for estimating the size of a population is recorded.
A population size estimation device according to the present invention sorts the measurement data in a population size estimation device that performs estimation of the total number of measurement data obtained by a measurement method having a detection limit. Means for preparing a plurality of candidate values for the total number of data in the range calculated from the measurement data, and a theory relating to the distribution of the population of the measurement data for each candidate value Distribution preparing means for preparing a set of values, and evaluating the match between the sorted measurement data and the set of theoretical values by a function related to a correlation or linear relationship, and the candidate value that is optimally evaluated by the function It further comprises candidate value evaluation means for estimating the size of the population.

  According to the present invention, it is possible to accurately estimate the size of a population of measurement data having a detection limit and use it for standardization by setting a range of candidate values and comparing the assumed distribution with a scatter diagram. A possible population size estimation method can be provided.

It is a block diagram which shows the control structure of the analyzer 10 which concerns on embodiment of this invention. It is a flowchart of the population size estimation process which concerns on embodiment of this invention. It is a graph which shows the example of drawing of the measurement data 251 and scatter diagram concerning an embodiment of the invention. It is a graph which shows the result of having applied the population size estimation method to RNAseq in Example 1 which concerns on embodiment of this invention. It is a graph which shows the result of having compared standardization of RNAseq data by the population size estimation method with the conventional standardization in Example 2 which concerns on embodiment of this invention. It is a graph which shows the result of having investigated the coincidence of RNAseq data and DNA microarray data at the time of using the population size estimation method in Example 3 which concerns on embodiment of this invention. It is a graph which shows the result of having performed the sensitivity measurement of a DNA microarray using RNAseq data using the population size estimation method in Example 4 which concerns on embodiment of this invention. In Example 5 which concerns on embodiment of this invention, it is the histogram for computer experiments produced actually using the code of R language. It is a conceptual diagram explaining the conventional method of standardization of the measurement data with a detection limit.

<Embodiment>
[Control Configuration of Analysis Device 10]
First, with reference to FIG. 1, the control configuration of the analysis apparatus 10 (computer, population size estimation apparatus, analysis means) according to the embodiment of the present invention will be described.
The analysis apparatus 10 is a computing apparatus such as a PC / AT compatible machine or a general-purpose machine, and an OS such as Linux (registered trademark) or Windows (registered trademark) is installed therein.
The main components of the analysis apparatus 10 include a control unit 100 that is a control / arithmetic unit such as a CPU (Central Processing Unit) and a GPU (Graphics Processing Unit), a RAM (Random Access Memory), and a ROM (Read Only Memory). Storage device 110 which is a storage device / recording medium such as HDD, Hard Disk Drive (HDD), flash memory, SSD (Solid State Drive), and other external devices such as a pointing device such as a keyboard and a mouse, a touch panel, and a microarray analyzer An input unit 130 including an I / O interface, a liquid crystal display, an organic EL display, a display unit 140 such as a printer for printing, and 1000 Bass A LAN board or wireless LAN board, serial / parallel / USB interface such as standards such as -T and an input-output unit 150.
The analysis apparatus 10 is executed by the control unit 100 mainly using various programs stored in the storage unit 110 and data including a database and the like, so that the population size estimation method according to the embodiment of the present invention is performed. Can be realized using hardware resources.

Storage unit 110 stores a program and data for realizing the population size estimation method according to the embodiment of the present invention. The program and data in the storage unit 110 can be executed / processed by the control unit 100 using hardware resources.
The program and data include the database 200, the measurement data sorting unit 210 (measurement data sorting unit), the candidate value range setting unit 220 (candidate value preparation unit), the distribution creation unit 225 (distribution preparation unit), and the candidate value drawing unit 230 ( Candidate value drawing means) and a candidate value evaluation unit 240 (candidate value evaluation means, parameter calculation means).
The program and data may be separately recorded on a recording medium and configured to be installable in the analysis apparatus 10. It can also be downloaded from the Internet and installed.

  The database 200 is a part for storing a database such as SQL and various data. The database 200 mainly stores measurement data 251, distribution data 252 (theoretical value set), setting data 253, and original data 254.

The measurement data 251 is data created by sorting the values of the original data 254 in ascending order by the measurement data sorting unit 210.
The distribution data 252 is data related to a set of theoretical values of distributions that are assumed to be the same as the original data 254 and the measurement data 251 created by the distribution creating unit 225. The distribution data 252 may store a random value of an assumed distribution for the candidate value.
The setting data 253 includes a setting for creating the measurement data 251 from the original data 254, an assumed distribution type of the measurement data 251, a parameter for setting a candidate value range, a candidate value range parameter, and a candidate value. Data, a set of candidate values, and the like.
The original data 254 is original data acquired from experimental equipment using a measurement method with a detection limit. As the original data 254, for example, quantification data of the gene expression level related to the transcriptome, mass measurement data such as proteins and low molecular organic compounds, and the like can be used. Examples of the quantification data of gene expression level include the number of RNAseq reads, DNA microarray signal intensity, RT-PCR expression level value, real time corresponding to the RNA of the gene to be analyzed and the amount of translated protein. The RT-PCR fluorescence value (fluorescence intensity ratio), the number of tags measured by the SAGE method (Serial Analysis of Gene Expression), derived data thereof, and the like can be used. As mass measurement data having a detection limit, various mass measurement data such as LC-MS data can also be used. Furthermore, data of various biological and non-biological values obtained as a set of numerical values whose measurement values are determined can be used as the original data 254.

The measurement data sorting unit 210 performs initialization such as sorting the original data 254 acquired from the experimental equipment by size, and stores the data as measurement data 251 in the storage unit 110.
In addition, the measurement data sorting unit 210 calculates statistical values such as a maximum value, a minimum value, an average, a median, and a variance from the measurement data 251 in order to determine a range of candidate values.

The candidate value range setting unit 220 is a part that sets a range of theoretical values (hereinafter referred to as “candidate values”) that are assumed as the original complete data number of the measurement data 251.
The candidate value range setting unit 220 sets a new candidate value range in accordance with the evaluation of the candidate value evaluation unit 240. The setting of the range of candidate values will be described later.

The distribution creation unit 225 creates a set of distribution values (theoretical values) corresponding to the set candidate values as the distribution data 252 according to the assumed distribution type setting of the setting data 253. Here, the distribution creation unit 225 can calculate a set of theoretical values of the distribution at a predetermined z score and store it in the distribution data 252.
For example, when the measurement data 251 is assumed to be a normal distribution, the distribution creation unit 225 prepares a standard normal distribution, calculates a set of theoretical values, and stores the set in the distribution data 252. That is, for example, the distribution creating unit 225 prepares a distribution of μ = 0 and σ = 1 for a normal distribution.
Further, for example, when the measurement data 251 is assumed to be a log normal distribution, the distribution creation unit 225 converts the created normal distribution to a log normal distribution, and distributes a set of log normal distribution values as a set of theoretical values. Store in data 252.
The distribution creation unit 225 can also directly create a set of normal distribution values corresponding to the candidate values using random numbers or the like. In this case, the distribution creating unit 225 can sort the created set of distribution values and store them as distribution data 252.
In addition, the distribution creation unit 225 stores the distribution function related to the assumed distribution type and the probability density function itself in the distribution data 252, and calculates the theoretical value of the distribution at a predetermined z score for each candidate value each time. Also good.

The candidate value drawing unit 230 is a part that creates a scatter diagram of the distribution data 252 and the measurement data 251.
For example, the candidate value drawing unit 230 creates and draws a scatter diagram by associating a set of distribution values of the prepared distribution data 252 with a set of values obtained by standardizing the measurement data 251 with each candidate value in descending order. To do.

In this embodiment, the candidate value drawing unit 230 calculates a so-called z-score and performs standardization using the estimated parameter when the data has a normal distribution. That is, the z-score is calculated and standardized by the following equation (1):

z = (measured value−m) / s (1)

In the above equation (1), m is an estimated value of the population mean μ, and s is an estimated value of the population standard deviation.

The candidate value evaluation unit 240 is a part that evaluates the relationship between the distribution data 252 and the measurement data 251.
The candidate value evaluation unit 240 evaluates the match between the sorted measurement data 251 and the distribution data 252 related to the set candidate value by “linearity”, “theoretical ranking value—the distance between the sorted data” or the like. To do.
The candidate value evaluation unit 240 can calculate and evaluate the Pearson product-moment correlation coefficient r value, for example, as a linearity evaluation function.
Further, the candidate value evaluation unit 240 uses the Mahalanobis general distance such as the sum and average of absolute values of the differences between the ranks of the distribution data 252 and the measurement data 251 and the sum and average of the sum of squares as the distance between the sort data. Can do.
The candidate value evaluation unit 240 determines the candidate value with the highest evaluation as an estimated value of the size of the population.
In addition, the candidate value evaluation unit 240 can also calculate a parameter indicating the distribution characteristics of the population using the estimated size of the population.

[Population size estimation process]
Next, a population size estimation process for executing the population size estimation method according to the embodiment of the present invention will be described with reference to FIGS.
Here, in the population size estimation processing according to the present embodiment, the number of types of objects to be measured is estimated as “the size of the population”. As the “population size”, for example, in the case of SDS-PAGE experimental data, which is a kind of electrophoresis, proteins with high hydrophobicity or proteins with high charge are not electrophoresed in the first place and cannot be detected. It is preferred not to enter the group. That is, in this example, all the proteins that can be detected by SDS-PAGE in principle are the population.
That is, when there are many types of objects to be measured and some of them can be measured by a predetermined experimental method, the number of types that can be measured becomes the size of the population.

In the population size estimation process of the present embodiment, the total number of data n is an unknown, and the most suitable candidate value is searched for while substituting various candidate values for the total number of data n. In order to confirm the suitability of the candidate value, the coincidence between the theoretical value and the data is evaluated by a function related to the correlation or linear relationship. In this case, the correlation or linear relationship is evaluated by “linearity”, “theoretical value of rank—distance between sorted data”, or the like. As the distance between the sorted data, Mahalanobis distance, “norm”, or the like can be used.
Here, when the theoretical value and the evaluated candidate value are drawn in a scatter diagram or the like which is a modification of the QQ plot (Quantile-Quantile plot), the evaluation becomes higher as the line approaches, and the evaluation is highest for the total number of data n The value is appropriate. That is, it is an estimation result of a desirable “population size”.

That is, in this embodiment, the type of distribution assumed is known, and the size of the population of measurement data 251 obtained by a measurement method with a detection limit is estimated. Thereby, the measurement data 251 is standardized. For this purpose, the total data number n assumed as the original complete data number of the measurement data 251 is unknown, and the value is analytically solved.
More specifically, the analysis apparatus 10 first sorts the original data 254 using the measurement data sorting unit 210 to create measurement data 251.
Next, the analysis apparatus 10 uses the candidate value range setting unit 220 to determine a range of candidate values for the total number of data n. Then, the analysis apparatus 10 uses the candidate value range setting unit 220 to prepare a plurality of candidate values from within the range.
Then, the analysis apparatus 10 causes the candidate value drawing unit 230 to correspond the distribution data 252 that is an assumed distribution and the measurement data 251 standardized with each candidate value in the range in order of increasing, and transform the QQ plot. Draw a scatter diagram.
Next, the analysis apparatus 10 evaluates the coincidence between the theoretical value related to the total number of data n and the measurement data 251 data using the evaluation function of the candidate value evaluation unit 240. It is also possible to directly evaluate coincidence between data without drawing a scatter diagram.
When the candidate value range setting unit 220 calculates the end of the range as the optimal candidate value or the prepared candidate value has low accuracy, the analysis apparatus 10 performs the calculation while adjusting the range of the candidate value. repeat. Then, the analysis device 10 determines the candidate value when sufficient accuracy is obtained. This determined candidate value becomes the expected value of the original number n of total complete data of the measurement data 251.
Hereinafter, with reference to the flowchart of FIG. 2, each population size estimation process will be described in detail for each step.

(Step S101)
The control unit 100 performs a data reading sort process using the measurement data sort unit 210.
In this process, first, the control unit 100 acquires original data 254 obtained by a measurement method having a detection limit. The original data 254 may be stored in the storage unit 110 from another recording medium. Further, the original data 254 may be acquired from another server on the network, experimental equipment, or the like (not shown) via the input / output unit 150.
The control unit 100 rearranges (sorts) a set of values to be standardized, such as the expression level of the acquired original data 254, in order of decreasing or increasing order according to size.
Then, the control unit 100 stores the sorted set of values in the storage unit 110 as measurement data 251.
Further, the control unit 100 obtains statistical values such as the maximum value, minimum value, average, and variance (standard deviation) of the measurement data 251 and stores them in the storage unit 110.
3A, the control unit 100 draws a histogram with the z score as the horizontal axis and the frequency of values as the vertical axis when the measurement data 251 is standardized as a standard normal distribution, or A frequency distribution table may be created and stored in the storage unit 110. FIG. 3A shows an example in which normal distribution data similar to FIG. 9B is drawn.

Regarding the rearrangement of data, for example, when the “R language” that is a statistical programming language is used, the following codes can be used:

sorted_data <-sort (data, na.last = T, decreasing = T)

Here, the array data is rearranged and assigned to the array sorted_data.

(Step S102)
Next, the control unit 100 uses the candidate value range setting unit 220 to perform candidate value range initial setting processing.
In this process, the control unit 100 uses the measurement data 251 to determine an initial value of a “candidate value” range that is an estimated value of the total number n of measurement data 251.
As the range of candidate values, the control unit 100 first acquires an assumed distribution type from the setting data 253. Then, using a statistical value of the measurement data 251, a possible range from the assumed distribution is set as an initial value and stored in the setting data 253.
As the range of candidate values, for example, when the measurement data 251 is standardized, an integral value of values that have a distribution of −4 to −1 in z score can be calculated. In addition, a predetermined value may be prepared for each type of data, or an initial value in a range of candidate values may be prepared from the same type of data measured previously.

(Step S103)
The control unit 100 uses the candidate value range setting unit 220 to perform candidate value preparation processing.
Here, the control unit 100 prepares a plurality of candidate values from the set range of candidate values and stores them in the setting data 253 of the storage unit 110.
For example, the control unit 100 can equally divide the set range of candidate values by a predetermined value and prepare candidate values for the predetermined value + 1. That is, in the case of the above example, when the predetermined value is 100, the control unit 100 can prepare 101 candidate values by dividing the range of 0 to 650 into 100 equal parts.
If the measurement data 251 is an integer such as the number of leads, the total number of data n is assumed to be an integer, and integer candidate values can be prepared.

For this candidate value preparation process, for example, in the case of preparing 101 candidate values in increments of 20 in the range of 9000 to 11000 as R language codes, the following codes can be used:

n_candidates <-seq (from = 9000, to = 11000, by = 20)

Here, n_candidates is an array of candidate values.
Further, in the linearity evaluation process described later, the linearity of the scatter diagram is determined by calculating, for example, the Pearson r value as the value of the evaluation function. For this reason, a vector rs for storing the value of the evaluation function related to each candidate value is also prepared like the following R language code:

rs <-1: 101 + NA

(Step S104)
The control unit 100 uses the distribution creation unit 225 to perform distribution preparation processing.
In this process, first, the control unit 100 selects one candidate value from the prepared candidate values.
Next, the control unit 100 refers to the setting data 253 and prepares a set of theoretical values of distribution assumed as the measurement data 251.
For the normal distribution, the control unit 100 prepares, for example, a standard normal distribution with N (0, 1), that is, an average of 0 and a variance σ ² = 1, and stores the data as distribution data 252 in the storage unit 110.
Moreover, the control part 100 uses the thing of each standard state, when using a different distribution model. For example, in the case of a lognormal distribution, the control unit 100 prepares a lognormal distribution with γ = 0, μ = 0, and σ = 1, and stores the logarithmic normal distribution in the storage unit 110 as distribution data 252.

(Step S105)
The control unit 100 uses the candidate value drawing unit 230 to perform a scatter diagram drawing process.
Specifically, first, the control unit 100 draws a scatter diagram by associating the measurement data 251 with the prepared distribution data 252 in descending order.
Here, in the conventional QQ plot, the data at the same quantile are associated with each other. On the other hand, the control part 100 draws the figure which matched the number in the same order | rank as a scatter diagram of this embodiment.
For this reason, in the scatter diagram of the present embodiment, in the lowest region, a situation occurs in which the theoretical value or the sorted data does not have a counterpart to be associated. In this case, the control unit 100 draws using only a pair that can be matched.

  If it demonstrates in detail with reference to FIG.3 (b)-(d), when the control part 100 gives a certain probability p, q1 and q2 used as two probability points (quantile) will be each, for example, respectively. A scatter diagram (probability plot) plotted for the X axis (horizontal axis) and the Y axis (vertical axis) is created. Here, the distribution data 252 is drawn as q1. In addition, the control unit 100 standardizes the measurement data 251 with the selected candidate value and renders it as q2. This is a QQ plot (see Gnanadesikan, R .; Wilk, MB (1968), "Probability plotting methods for the analysis of data", Biometrica 55 (1): variants 1 to 17). . That is, in this scatter diagram, the distribution data 252 that is the theoretical value of the distribution prepared for the X-axis coordinates is compared with the measurement data 251 that is the data sorted on the Y-axis coordinates. At this time, the distribution data 252 is correlated in ascending order and the measurement data 251 is correlated in ascending order to examine the relationship.

For example, the following code can be used as an R language code for displaying the scatter diagram:

ideal <--qnorm (ppoints (n_candidates [i]))
s_limit <-min (length (data [! is.na (data)]), n_candidates [i])
plot (sorted_data [1: s_limit], ideal [1: s_limit])

Referring back to FIG. 3, FIG. 3B shows an example of a scatter diagram when the candidate value is 7000 in the data of FIG.
FIG. 3C shows an example of a scatter diagram when the candidate value is 10,000, as in FIG.
FIG. 3D shows an example of a scatter diagram when the candidate value is 14000, similarly to FIGS. 3B and 3C.

Note that the control unit 100 uses the distribution creation unit 225 to directly create a set of distribution values corresponding to the selected candidate value from the (pseudo) random number sequence in the above-described distribution preparation process, and generate this distribution data. It can also be stored as 252. At this time, the control unit 100 can create a normal distribution by using a Box = Muller transform, a method using a central limit theorem, or the like. The control unit 100 can also create a lognormal distribution from the normal distribution.
The control unit 100 also uses the candidate value drawing unit 230 to perform a process of directly drawing a scatter diagram in which the distribution data 252 that is a random number sequence having such a distribution and the measurement data 251 are associated with each other. Is possible.
Moreover, the control part 100 can also perform the following candidate value evaluation process directly, without drawing this scatter diagram.

(Step S106)
The control unit 100 uses the candidate value evaluation unit 240 to perform candidate value evaluation processing.
In this process, the control unit 100 evaluates the coincidence between the distribution data 252 and the measurement data 251 with a function related to a correlation or linear relationship. As a function related to this correlation or linear relationship, a predetermined evaluation function related to “linearity” or “theoretical value of rank—the distance between sorted data” or the like can be used.
When evaluating by linearity, as this evaluation function, the control part 100 can calculate and use Pearson's r value, for example.
Further, when calculating the Pearson r value, the control unit 100 can also estimate μ and σ as the y-intercept and the slope of a straight line, respectively, in the case of a normal distribution.

Referring again to FIG. 3, in the scatter diagram of FIG. 3 (b), Pearson's r value was 0.9844363. The estimated values of μ and σ of the standardized measurement data 251 were (0.46, 0.71).
In FIG. 3C, Pearson's r value was 0.9998996, and the estimated values of μ and σ of the measurement data 251 used for standardization were (−0.01, 0.99).
In FIG. 3D, Pearson's r value was 0.9981416, and the estimated values of μ and σ of the measurement data 251 used for standardization were (−0.43, 1.14).
As described above, the data of FIG. 3A has a large Pearson r value and high linearity when the candidate value is 10,000 as shown in FIG. It can be seen that the estimated value is good.

  When the control unit 100 evaluates the distribution data 252 and the measurement data 251 based on “theoretical value of the rank—the distance between the sort data”, the total of the absolute values of the differences between the values of the ranks is used as the “distance”. And Mahalanobis general distance such as average, sum of squares, average, etc. can be calculated and used.

(Step S107)
The control unit 100 uses the candidate value evaluation unit 240 to determine whether or not all candidate values prepared within the range have been evaluated.
If Yes, that is, if all candidate values have been evaluated, the control unit 100 advances the process to step S108.
No, that is, if there is a candidate value within the range, the control unit 100 returns the process to step S104, and continues preparation of distribution, drawing of a scatter diagram, and evaluation of linearity. At this time, for example, evaluation may be performed alternately for a large value and a small value within a range of candidate values, and the range may be narrowed to gradually obtain the most probable candidate value, and the calculation may be performed at high speed. Further, Newton's method may be used.

As an R language code for such linearity evaluation, for example, the following codes can be used:

for (i in 1: 101) {
n <-n_candidates [i]
ideal <--qnorm (ppoints (n)) # ideal value
s_limit <-min (length (data [! is.na (data)]), n)
rs [i] <-cor (sorted_data [1: s_limit], ideal [1: s_limit])} # Peasron r value

This code shows an example in which all candidate values are simply evaluated and substituted for the vector rs prepared in the above code.
On this basis, the optimum code is selected from the candidate values using the following code. Here, the control unit 100 selects a candidate that gives the largest evaluation value.

selected <-n_candidates [which (rs == max (rs))]

The variable selected is a candidate value that is most likely to be evaluated in the vector rs.

Referring to FIG. 3 (e), an example of a change in the Pearson r value when a scatter diagram is drawn with 0 to 650 as a range of candidate values for measurement data 251 different from FIG. 3 (a) is shown. ing. In FIG. 3E, the X-axis is each candidate value, and the Y-axis is Pearson's r value.
In this way, the control unit 100 monitors the relationship linearity with an evaluation function such as the Pearson r-value as described above while changing the value of the candidate value, and searches for the most likely candidate value. Can do. This relationship is considered to be simply convex upward or downward.
For this reason, Newton's method can be applied to the task of obtaining the most probable candidate value, and it is possible to save computational resources.

  If this figure is not convex upward, it is considered that the candidate value is inappropriate, and therefore the control unit 100 performs a candidate value range shifting process described later. Alternatively, the control unit 100 performs a candidate value range narrowing process for selecting candidates from a narrower range.

It is also possible to evaluate the above scatter plot using another evaluation function that is not the Pearson r value.
Referring to FIGS. 3 (f) and 3 (g), for example, in the R language code described above, the place to be assigned to the vector rs can be changed as follows:

rs [i] <--mean ((sorted_data [1: s_limit] -ideal [1: s_limit]) ^ 2) # Another evaluation function

In FIG. 3F, the X-axis represents each candidate value, and the Y-axis represents the value of this other evaluation function. As described above, since the value changes depending on the evaluation function, the evaluation function can be selected in accordance with the data characteristics.
FIG. 3G is an example in which a scatter diagram is drawn with the following code for the variable selected to which the candidate value having the maximum value of the other evaluation function is substituted, as in FIG.

ideal <--qnorm (ppoints (selected))
s_limit <-min (length (data [! is.na (data)]), selected)
plot (sorted_data [1: s_limit], ideal [1: s_limit])

(Step S108)
The control unit 100 uses the candidate value evaluation unit 240 to determine whether or not there is an optimal candidate value on the outside (end) of the prepared candidate value.
Here, the control unit 100 determines that there is a more optimal candidate value when the graph of the evaluation function described above has almost no extreme value and monotonously increases / decreases, that is, Yes. That is, the control unit 100 determines Yes when the value of the evaluation function is the highest value at the end of the prepared candidate value. Otherwise, it is determined as No.
In the case of Yes, the control unit 100 proceeds with the process to step S109 and changes the range of candidate values.
In No, the control part 100 advances a process to step S110.

(Step S109)
When the candidate value range needs to be changed, the control unit 100 uses the candidate value range setting unit 220 to perform candidate value range shifting processing.
Here, if the value of the evaluation function is the highest value at any end of the prepared candidate values, the range of candidate values is adjusted and shifted to that side. That is, for example, when the prepared candidate value range is 500 to 1000, and the value of the evaluation function of the candidate value 1000 is the maximum value, processing for resetting the candidate value range to 1001 to 1501 is performed. Do.
The control unit 100 then returns the process to step S103 and proceeds with the calculation again within the range of the shifted candidate values.
When resetting candidate values, the range of the candidate values can be set based on the rate of increase / decrease of the graph of the evaluation function, instead of shifting the range of the candidate values by the same width.

(Step S110)
When there is a most probable candidate value within the range of candidate values, the control unit 100 uses the candidate value evaluation unit 240 to determine whether more accuracy is necessary.
For example, the control unit 100 refers to the prepared candidate value range width and the candidate value range parameter of the setting data 253, and if the search needs to be performed with a narrower width, Yes. Is determined. In other cases, for example, candidate values are prepared as integers, and if there is an optimum candidate value in the evaluated consecutive integers, it is determined as No.
In Yes, control part 100 advances a process to Step S111.
In No, the control part 100 advances a process to step S112.

(Step S111)
The control unit 100 uses the candidate value range setting unit 220 to perform candidate value range narrowing processing.
The control unit 100 narrows and sets the range of candidate values again. Here, the control unit 100 can set values before and after the candidate value having the highest value of the linearity evaluation function as a range of new candidate values.
The control unit 100 then returns the process to step S103 and proceeds with the calculation again within the narrowed range of candidate values.

(Step S112)
When the optimal candidate value is obtained with sufficient accuracy, the control unit 100 uses the candidate value evaluation unit 240 to perform candidate value determination processing.
Here, the control unit 100 determines the candidate value when sufficient accuracy is obtained.
This is estimated as the expected value of the original complete data number n of the measurement data 251, that is, the size of the population.

(Step S113)
Next, the control unit 100 uses the candidate value evaluation unit 240 to estimate a parameter based on the determined candidate value, and standardizes and stores the measurement data 251 using this parameter.
At this time, the control unit 100 can perform calculation using, for example, linear approximation while considering the number of lost data.
At this time, the control unit 100 compares the sort data and the theoretical value with values in the same order from the top, and performs a first-order approximation to calculate the intercept and the slope. The control unit 100 estimates the parameter using the calculated intercept and slope. The controller 100 can use a least square method, another robust method, or the like in this linear approximation.

For example, the following code can be used as the R language code for setting parameters for parameter estimation:

u_limit <-selected-length (data [! is.na (data)])
params <-coef (lm (sorted_data [u_limit: s_limit] ~ ideal [u_limit: s_limit]))

On top of this, the estimated values for μ and σ are:

(mu <-params [1]) # μ
(sigma <-params [2]) # σ

It can be calculated with a code such as That is, in the case of a normal distribution, μ is an intercept and σ is a slope as in this example.
The control unit 100 standardizes and stores the measurement data 251 using these calculated parameters.

(Step S114)
The control unit 100 uses the candidate value evaluation unit 240 to perform analysis processing.
Next, if necessary, the control unit 100 compares the actual measurement data with other standardized measurement data 251. For example, the control unit 100 can compare RNAseq, which is the same expression level data, with data of the DNA microarray. As a result, the sensitivity of each spot of the DNA microarray can be obtained by correcting with the normalized gene expression level of RNAseq.
The control unit 100 can also perform data analysis such as principal component analysis.
The population size estimation process is thus completed.

With the configuration described above, the following effects can be obtained.
As described above, it is difficult to accurately estimate the parameter with the data of the measurement method that loses data due to the detection limit.
On the other hand, the population size estimation method according to the embodiment of the present invention uses the size of the population based on the relationship between the distribution model of the measurement data 251 of the measurement method having a detection limit and the assumed distribution. A total number n of data can be estimated, and a population parameter such as μ or σ can be obtained. This increases the accuracy and accuracy of standardized values. In particular, in the measurement data 251 having a lognormal distribution, the total number of data n is important in accurately obtaining the parameter. Accordingly, the present invention can accurately standardize the data of the measuring method having a detection limit.
That is, by estimating and standardizing the population size by the population size estimation method of the present embodiment, RNAseq, metabolome analysis, transcriptome analysis such as microarray using qPCR, data analysis of next generation sequencer, etc. Can be performed more accurately than in the past.
Although the detection limit may change in measurement, the population size estimation method according to the present embodiment can compare standardized values even in such a case. Therefore, for example, it is possible to increase the measurement items from which a significant difference is obtained.
In addition, with the population size estimation method according to the embodiment of the present invention, more accurate analysis can be performed by standardizing transcriptome data relating to diseases and the like and using it for principal component analysis or the like. .

Conventionally, there was no method for measuring the sensitivity of each spot of a DNA microarray. Therefore, it was difficult to directly compare the data of DNA microarrays with different designs, even if they measured the same gene. This is because the base sequence used for the probe differs depending on which sequence of the gene is measured, and the sensitivity of the probe differs.
The population size estimation method of the present embodiment can estimate the sensitivity of each spot by standardizing such DNA microarray data in correspondence with another measurement method. Thus, by obtaining the sensitivity of each spot, it becomes possible to correct data of DNA microarrays of other organs and patients. Thereby, the microarray data can be used as a kind of absolute value. For this reason, it becomes possible to directly compare data of chips having different designs.

Also, DNA microarrays cannot detect RNA that is not on the chip design. For this reason, there has been a problem that it is not possible to measure new / unknown objects such as RNA that does not encode proteins, which has been particularly important in recent years.
On the other hand, the population size estimation method of the present embodiment can be applied to analysis of a transcriptome or the like that can measure a new / unknown target such as RNAseq. Further, by comparing with well-known DNA microarray data, the expression level of the new / unknown target can be reliably measured. In addition, RNAseq and the results can be directly compared regardless of the sensitivity difference in the measurement results of the DNA microarray.

Note that the population size estimation method of the present embodiment can be applied to other than biological data.
For example, the present invention can be applied to a case where substances contained in food, soil, water, and air are identified and quantified in combination with a separation method of substances such as chromatography and a detection method such as mass spectrometry. In that case, information on the distribution form of these substances is required. For example, since a food is derived from living organisms, a lognormal distribution can be used.

  EXAMPLES Next, although an Example demonstrates this invention further based on drawing, the following specific examples do not limit this invention.

[Application to standardization of RNAseq data]
RNAseq is a transcriptome measurement method that reverse-transcribes RNA into DNA and determines the sequence of the DNA fragment. In RNAseq, since a complementary sequence of RNA present in a larger amount is observed at a higher frequency, the concentration of each RNA is estimated from this frequency. RNAseq can be measured even with unknown RNAs with different splicing and the like. In addition, it is a digital measurement method in which data is observed / not observed, and since there is no intermediary such as an image, a more accurate measurement result is expected.
However, RNAseq is multivariate data with large measurement items, and analysis is not easy. It also has binomial noise, derived from Bernoulli trials that count the number of observations. Initially, it was even said that standardization was unnecessary, but of course, standardization is necessary for data comparison in order to match the total number of leads.
Therefore, in this example, the RNAseq data was standardized by estimating the population size by the population size estimation method of this embodiment.

(Use data)
The data “John C. Marioni et. Al., RNA-seq: An assessment of technical reproductivity and comparison with gene expression arrays, Gen. Res. Used as 254.
This is data obtained by measuring kidney and liver samples from a single human specimen with RNAseq using a next-generation sequencer manufactured by Illumina and a DNA microarray manufactured by Affymetrix.

(result)
For the liver (R1L4Liver) and kidney (R1L1Kidney) samples, specify 12000 to 23000 as the initial range, use the lognormal distribution as the distribution of the original data 254, and use the population size estimation method of this embodiment to set each candidate value. A scatter diagram was drawn and Pearson's r value was measured. 10 was used for the base of the logarithm.
As a result, the optimal candidate value for the kidney sample was 12620. That is, the total number of data n is estimated to be 12620. The population of this kidney sample had an average μ = 1.77 and a standard deviation σ = 0.512.
In the liver sample, the optimal candidate value was 13060. That is, the total number of data n is estimated to be 13060. The parameters of this liver sample were mean μ = 1.515 and standard deviation σ = 0.610.
The relationship between the measurement data 251 and the distribution model of the distribution data 252 will be described below with reference to FIG.

FIG. 4A is a scatter diagram, which is a modification of the above-described QQ plot, drawn for a liver sample. In FIG. 4A, the X-axis is a z-score of a prepared normal distribution, and the Y-axis is a z-score of standardized measurement data 251 created from a liver sample. Each point indicated by a circle indicates a calculated probability point.
FIG. 4B is a histogram when the liver sample is standardized. The X axis shows the z-score of the prepared normal distribution, and the Y axis bar shows the frequency of each sample. The dotted line indicates the frequency expected from the optimal candidate value distribution.
As shown in the results of FIGS. 4A and 4B, the data up to about 0.5 below the z score may coincide with the model of the parameter obtained from the calculated optimal candidate value. Recognize. This corresponds to about 40% of the total. Data with a z-score of 0.5 or less is below the measurement limit and is considered to be affected by noise from the sample and noise from the measuring instrument / experimental method according to the binomial distribution.

FIG. 4C is a diagram in which the same liver samples as in FIGS. 4A and 4B are directly plotted in order of size, not a scatter diagram. Also in this figure, the X-axis is a prepared normal distribution z-score, and the Y-axis is a standardized measurement data 251 z-score. In this figure, data for the optimum candidate values calculated in accordance with the distribution data 252 are created, arranged in order of size, and drawn and compared with the measurement data 251 in pairs.
Similar to FIGS. 4 (a) and 4 (b), the result is that the data of 0.5 or less is rather affected by noise, and reproducibility of the data cannot be expected.
For this reason, in this sample, it is also possible to compare the expression levels without using data with a z score of 0.5 or less.

4 (d) and 4 (e) are similar to FIGS. 4 (a) and 4 (b) described above, with respect to the kidney sample, the measurement data 251 and the normal distribution standardized by the population size estimation method of the present embodiment. The comparison by a scatter diagram and a histogram are shown.
As a result, when the z-score is about −1, a deviation from the model of the parameter obtained from the calculated optimal candidate value is observed. Data below this suggests that it is affected by noise from the sample or by binomial noise.

[Comparison with conventional standardization method of RNAseq data]
Next, referring to FIG. 5, a conventional standardization method (hereinafter referred to as “constant summation”) that is often used, using the same RNAseq data as the above-mentioned Marioni et al. The expression level was compared with the population size estimation method of this embodiment.
The conventional “total sum constant” standardization method is a standardization method in which the whole is divided by the total number of samples (leads) and the intensities of the measured values are combined.

(Gene selection)
The liver and kidney RNAseq data measured 7 times were t-tested for each gene, and there was a difference in the expression level (number of reads) at a significance level of 1% on both sides in either the liver or kidney. A gene that was positive was selected.
Focusing on these genes, a combination of expression levels paired in the liver and kidney was obtained. That is, the logarithmic ratio (log ratio) of the expression levels of the liver and kidney was compared for genes that became 1% significant by statistical test.

(result)
FIG. 5 (a) shows the difference between the kidney and the liver from two sets of data using the expression levels of each gene standardized by the conventional standardization method with the sum total divided by the total number, and compared each. It is a plot. The X-axis and Y-axis indicate the difference in gene expression between the liver and kidney using the log ratio, respectively. Each point is a plot of the difference in the expression level of the kidney and liver of each gene as one data combination on the X coordinate and the other combination on the Y coordinate. “R” in the figure indicates the Pearson correlation coefficient r value.
As shown in FIG. 5 (a), it is understood that the reproducibility of the expression change is not so high when the summation is constant. In particular, there are some genes that vary in the y = −x direction near the origin (X coordinate, Y coordinate) = (0, 0), that is, exhibit an inverse correlation. This indicates a tendency that the test tends to generate false positives, particularly in a gene group with a small difference in expression level. This phenomenon occurs when there are many noise components in the measurement.

FIG. 5B shows the difference in the expression level between the kidney and the liver in the measurement data 251 in which the expression level is standardized by the population size estimation method of the present embodiment in the same two sets of measurements as in FIG. It is the compared plot. The X-axis and Y-axis indicate the difference in gene expression between the liver and kidney using the log ratio, respectively. Each point is a plot of the expression level of each gene in the kidney and liver in the same manner as in FIG. 5A, and “r” in the figure is the Pearson correlation coefficient r value.
In FIG. 5B, each plot is separated into two on the straight line y = x with the origin as a boundary. That is, almost no genes that are inversely correlated are detected, indicating that the test is unlikely to produce false positives. This is a phenomenon that occurs when the noise component of the measurement result is smaller than the signal. The reproducibility of measurement is high, and the r value is increased from the conventional 0.897 to about 0.949.
It can also be seen that the number of genes having a large value at y = x is clearly larger in the plot of the population size estimation method of the present embodiment. This is thought to be because the number of genes for which a significant difference was detected by the test increased because the noise component decreased. That is, as a result, not only the correlation coefficient is increased, but also the value of the data can be increased.

[The agreement between RNAseq data and DNA microarray data]
Next, regarding the same RNAionq and DNA microarray data of Marioni et al. As in Examples 1 and 2, the RNAseq measured for the expression level of the same gene and the DNA microarray data were compared. Specifically, it was examined whether or not the population size estimation method of this embodiment affects.
In Example 3, the same data as in Example 2 above, either in the expression level of RNAseq in either liver or kidney, a t-test was performed at a significance level of 1% on both sides, and there was a difference in expression level (number of reads). The gene expression level was used as the original data 254. In addition, a DNA microarray having a difference in expression level at a 1% significance level was selected and used as original data 254.

  FIG. 6A is a control in which original data of the original Marioni et al. Is plotted. The Y-axis is the logarithmic ratio of the expression level in RNAseq, and the X-axis is the expression level data of DNA microarray manufactured by Affymetrix. This is a value calculated by standardizing DNA microarray expression data using a conventional RMA (Robust Multi-Array Average) method. Genes that have 250 or more reads in RNAseq are indicated by dark colored dots. Less data is shown in light-colored (lightly filled) dots. Genes with 250 reads or more are limited to the top 6.4% of data. Pearson's r value was about 0.87.

FIG. 6B shows the same gene plotted as in FIG. 6A after standardization by the population size estimation method of the present embodiment. In particular, there is less vertical variation at a relatively low number of leads of less than 250 leads, indicated by the light color (light fill). Specifically, Pearson's r value is improved to about 0.914.
That is, even in a gene with a small expression level, the expression level of RNAseq and DNA microarray can be analyzed with high reproducibility and with less error. It is also possible to compare the expression levels of RNAseq and DNA microarray.

[Measurement of DNA microarray sensitivity]
Next, using the same data as in Example 3 above, the sensitivity of each spot of the DNA microarray was measured by comparing RNAseq with the DNA microarray.
In the analysis of conventional DNA microarray expression level data, the amount of change in gene expression is usually calculated by calculating the “magnification” between the measurement level of the control tissue and the measurement level when measuring the expression level of the tissue to be examined. Seeking. For this reason, the difference in sensitivity caused by the Tm value (melting temperature) or yield of the DNA sequence formed in each spot of the DNA microarray can be offset. However, there has been a desire to measure the sensitivity of each spot in order to obtain a more general quantitative value.
In contrast, in this example, the sensitivity of each spot of the DNA microarray was measured using RNAseq data of the same gene.

(Estimation of sensitivity difference)
The liver data of Example 3 described above is standardized for RNAseq and DNA microarrays using the population size estimation method of this embodiment, and the sensitivity of spots corresponding to each gene in each DNA microarray is estimated.
The sensitivity difference of each spot of each DNA microarray is corrected using the sensitivity of each spot. Specifically, for each spot, the difference in sensitivity obtained by liver RNAseq and DNA microarray is estimated, and the value of this sensitivity difference is multiplied by the kidney data and standardized.

FIGS. 7A and 7B are values obtained by comparing the values obtained by standardizing the measurement values of the microarray on the Y axis and the values obtained by standardizing RNA-seq on the X axis by the population size estimation method of the present embodiment. It is. Both the X-axis and the Y-axis indicate the logarithmic value z-score, not the sensitivity magnification. FIG. 7A shows liver data, and FIG. 7B shows kidney data.
As a result, although there is a positive correlation, it varies greatly. That is, in Example 3, since the sensitivity was canceled by observing the magnification, the results matched. However, as described above, since microarrays have different sensitivities for each spot, it can be seen that the results do not always match when compared with absolute sensitivity.

FIG. 7 (c) estimates the sensitivity as “standardized data of RNA-seq-standardized data of microarray” based on the result of FIG. 7 (a), and this sensitivity is the result of the microarray of FIG. 7 (b). Is corrected by adding to That is, FIG. 7C is a plot corrected using the sensitivity difference of each spot (gene) estimated from the liver data based on the result of FIG. 7A. Similar to FIGS. 7A and 7B, both the X-axis and the Y-axis indicate absolute values (log scale).
As a result, the correlation with RNA-seq is clearly improved as compared with FIG. Thus, by measuring the sensitivity of the DNA microarray, it is possible to accurately measure the expression level for each gene while suppressing the dispersion of results. Further, if the number of data for sensitivity estimation is increased, the accuracy can be further increased.

[Example of implementation in R language]
Next, an example in which the estimation of the parameter according to the above-described embodiment is actually implemented in the R language and a computer experiment is performed will be described below.
First, 10,000 random numbers with normal distribution were prepared and stored in the array rand with the following code.
In addition, array data prepared by deleting data with a z score smaller than −0.5 from the array rand was prepared to simulate the measurement limit:

rand <-rnorm (10000) #
data <-rand; data [data <-0.5] <-NA

A histogram was created for the prepared random number rand and array data:

hr <-hist (rand, breaks = -50: 50/10) #Create histogram
hd <-hist (data, breaks = -50: 50/10)

On top of this, the histogram was overwritten with the following code:

plot (hr, col = "gray50", xlab = "", ylab = "", main = "")
par (new = T)
plot (hd, col = "white", xlab = "", ylab = "", main = "")

FIG. 8 is an example of a drawn histogram. The shaded portion in FIG. 8 corresponds to data that has not been detected, as in FIG.
Next, the parameter of the resulting array rand was estimated with the following code:

# Average
> mean (rand, na.rm = T)
[1] -0.0009032517

# Median (Robust μ estimate)
> median (rand, na.rm = T)
[1] 0.0079797

# Standard deviation (Robust μ estimate)
> sd (rand, na.rm = T)
[1] 0.9975465

# Standard deviation (Robust MAD estimation)
> mad (rand, na.rm = T)
[1] 0.9996847

The values of the parameters obtained from these arrays rand reflected the well-created normal distribution settings.
Next, the parameter of the obtained array data was estimated with the following code, which is a conventional method:

# Average
> mean (data, na.rm = T)
[1] 0.5033331

# Median (Robust μ estimate)
> median (data, na.rm = T)
[1] 0.3849361

# Standard deviation (Robust μ estimate)
> sd (data, na.rm = T)
[1] 0.6977276

# Standard deviation (Robust MAD estimation)
> mad (data, na.rm = T)
[1] 0.7157046

As described above, when the array data is calculated by the conventional method, it is far from the original normal distribution setting even by the robust calculation method.
On the other hand, the R language code of the above-described embodiment estimated the parameter using linear approximation while taking into account the number of lost data:

(mu <-params [1]) # μ
-0.0003264361

(sigma <-params [2]) # σ
0.9854621

  As a result, we succeeded in estimating the parameters close to the set values.

  Note that the configuration and operation of the above-described embodiment are examples, and it is needless to say that the configuration and operation can be appropriately changed and executed without departing from the gist of the present invention.

  The population size estimation method of the present invention can be widely used for transcriptome analysis using next-generation sequencers, DNA microarrays, and the like, and can be used industrially.

DESCRIPTION OF SYMBOLS 10 Analysis apparatus 100 Control part 110 Storage part 130 Input part 140 Display part 150 Input / output part 200 Database 210 Measurement data sort part 220 Candidate value range setting part 225 Distribution creation part 230 Candidate value drawing part 240 Candidate value evaluation part 251 Measurement data 252 Distribution data 253 Setting data 254 Original data

Claims (11)

In a population size estimation method that causes a computer to perform estimation of the total number of measurement data obtained by a measurement method with a detection limit,
Sort the measurement data,
In the range calculated from the measurement data, preparing a plurality of candidate values for the total number of data,
For each candidate value, prepare a set of theoretical values related to the population distribution of the measurement data,
The matching between the sorted measurement data and the set of theoretical values is evaluated by a function related to a correlation or linear relationship, and the candidate value that is optimally evaluated by the function is estimated as a population size. And population size estimation method. The population size estimation method according to claim 1, wherein the function evaluates linearity or distance between the value of the measurement data in the same order from the top and the value of the set of theoretical values. Using the Pearson product moment correlation coefficient as the function,
If the optimal candidate value is at the end of the range, or if the prepared candidate value is less accurate, adjust the range of the candidate value,
The population size estimation method according to claim 1, wherein the candidate value with the best evaluation is determined when sufficient accuracy is obtained. Using the set of theoretical values of the candidate values with the best evaluation according to any one of claims 1 to 3, and the sorted measurement data,
Comparing the value of the measurement data in the same order from the top and the value of the set of theoretical values, calculating a first-order approximation intercept and slope, and calculating a parameter indicating the characteristics of the population distribution A characteristic parameter estimation method. As measurement data, quantification data of gene expression level is used,
The gene expression level analysis method, wherein the gene expression level is standardized by the total number of data estimated by the population size estimation method according to any one of claims 1 to 3. The gene expression level analysis method according to claim 5, wherein the quantification data of the gene expression level is a transcript expression level, a cDNA read count, and derivative data. As measurement data, use quantification data of protein amount or substance amount,
The method for analyzing data of a proteome or metabolome, wherein the protein amount or the substance amount is standardized based on the total number of data estimated by the population size estimation method according to any one of claims 1 to 3. As measurement data, using quantification data of RNAseq and DNA microarray gene expression level,
Standardized by the total number of data estimated by the population size estimation method according to any one of claims 1 to 3,
A sensitivity measurement method for a DNA microarray, wherein the sensitivity of each spot of the DNA microarray is determined by correcting the standardized gene expression level of the RNAseq. In a program for causing a computer to execute a population size estimation method for estimating the total number of measurement data obtained by a measurement method with a detection limit,
In the computer,
Sort the measurement data;
In the range calculated from the measurement data, a plurality of candidate values for the size of the total data number are prepared,
For each candidate value, a set of theoretical values related to the population distribution of the measurement data is prepared,
A population size that evaluates a match between the sorted measurement data and the set of theoretical values by a function related to a correlation or linear relationship, and estimates the candidate value that is optimally evaluated by the function as the size of the population A program characterized by causing the estimation method to be executed. In a computer-readable recording medium recording a program for causing a computer to execute a population size estimation method for estimating the total number of measurement data obtained by a measurement method having a detection limit,
In the computer,
Sort the measurement data;
In the range calculated from the measurement data, a plurality of candidate values for the size of the total data number are prepared,
For each candidate value, a set of theoretical values related to the population distribution of the measurement data is prepared,
A population size that evaluates a match between the sorted measurement data and the set of theoretical values by a function related to a correlation or linear relationship, and estimates the candidate value that is optimally evaluated by the function as the size of the population A recording medium on which is recorded a program for executing the estimation method. In a population size estimation device that performs estimation of the total number of measurement data obtained by a measurement method with a detection limit,
A measurement data sorting means for sorting the measurement data;
Candidate value preparation means for preparing a plurality of candidate values of the size of the total number of data within a range calculated from the measurement data;
Distribution preparation means for preparing a set of theoretical values related to the population distribution of the measurement data for each candidate value;
Candidate value evaluation for evaluating coincidence between the sorted measurement data and the set of theoretical values by a function related to a correlation or linear relationship, and estimating the candidate value that is optimally evaluated by the function as the size of a population A population size estimation device.

Images