WO2014080447A1 - データ解析装置、データ解析方法 - Google Patents

データ解析装置、データ解析方法 Download PDF

Info

Publication number
WO2014080447A1
WO2014080447A1 PCT/JP2012/080003 JP2012080003W WO2014080447A1 WO 2014080447 A1 WO2014080447 A1 WO 2014080447A1 JP 2012080003 W JP2012080003 W JP 2012080003W WO 2014080447 A1 WO2014080447 A1 WO 2014080447A1
Authority
WO
WIPO (PCT)
Prior art keywords
data
cluster
clustering
sample data
experimental error
Prior art date
Application number
PCT/JP2012/080003
Other languages
English (en)
French (fr)
Inventor
白井　正敬
神原　秀記
妃代美 谷口
Original Assignee
株式会社日立製作所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 株式会社日立製作所 filed Critical 株式会社日立製作所
Priority to US14/443,118 priority Critical patent/US20150302042A1/en
Priority to JP2014548349A priority patent/JP6029683B2/ja
Priority to PCT/JP2012/080003 priority patent/WO2014080447A1/ja
Publication of WO2014080447A1 publication Critical patent/WO2014080447A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/23—Updating
    • G06F16/2365—Ensuring data consistency and integrity
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/28—Databases characterised by their database models, e.g. relational or object models
    • G06F16/284—Relational databases
    • G06F16/285—Clustering or classification
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
    • G16B40/30—Unsupervised data analysis
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics

Definitions

  • the present invention relates to a technique for clustering analysis of sample data.
  • Non-Patent Document 1 discloses a method that can measure a specific gene expression level with sufficient accuracy by a method using a quantitative PCR apparatus for gene expression analysis in a single cell.
  • Non-Patent Document 2 discloses a method for quantifying the expression of almost all genes using a large-scale DNA sequencer (next-generation sequencer) for gene expression analysis in a single cell.
  • Non-Patent Document 2 also discloses a data analysis method for specifying a cell type. In the future, genome sequences, intracellular proteins, and various biomolecules in cells are expected to be identified at the single cell level.
  • an algorithm for classifying cells using only data and an analysis / diagnosis apparatus using the algorithm do not always have sufficient characteristics to classify cells and use them in medical diagnosis.
  • necessary characteristics include, for example, cell grouping (hereinafter referred to as clustering), even if the optimal number of groups (hereinafter referred to as cluster number) is not known in advance. It is conceivable to perform appropriate clustering. In particular, it is difficult for a conventional analysis / diagnosis apparatus to determine whether an exceptional cluster having a small number of data is an independent cluster or a part of another cluster having a large number of data.
  • the present invention has been made in view of the above-described problems, and an object of the present invention is to provide a data analysis apparatus that can appropriately perform clustering even when exceptional data resulting from experimental error or the like occurs.
  • the data analysis apparatus determines a cluster range parameter for relaxing the cluster boundary in advance according to the range of the experimental error described by the experimental error data.
  • the clustering process for exception data that does not belong to any cluster, if any cluster includes an area further away from the exception data by a distance determined by the cluster range parameter, the exception data belongs to that cluster. If it is determined that it is not included in any cluster that is separated by the distance, it is determined that it constitutes an independent cluster.
  • the data analysis apparatus can appropriately classify cells using the analysis result of a single cell.
  • the number of types of classified cells can be determined with high accuracy.
  • FIG. 3 shows the results of calculating logarithmic likelihood representing the likelihood of the number of clusters by sweeping the CR value to 1, 2, and 4 for the simulated sample data described in FIG.
  • FIG. 6 is a functional block diagram of a data analysis apparatus 906 according to Embodiment 2.
  • FIG. It is a figure which shows the gene expression level of Bglap1, PParg, and Col2a1 of each cluster (indicated by numbers 1 to 5).
  • FIG. 13 is a diagram showing expression levels for the three genes described in FIG. 12 with and without clustering, with and without differentiation induction.
  • 10 is a functional block diagram of a data analysis apparatus 906 according to Embodiment 3.
  • FIG. 10 is a functional block diagram of a data analysis device 906 according to a fourth embodiment.
  • Cell data obtained by single cell analysis can be analyzed using, for example, main factor analysis. Data obtained by principal factor analysis is often used to visually determine groups. In order to specifically examine the problems of the conventional clustering method, simulated sample data will be used below. Specifically, assuming that a single cell analysis was performed on four genes and 180 cells, cell data was generated using random numbers on a computer.
  • FIG. 1 is a diagram showing a clustering result of simulated sample data.
  • FIG. 1A shows the configuration of each cluster in the simulated data.
  • ⁇ is a random number and follows a distribution of mean value 0 standard deviation ⁇ .
  • the random number ⁇ is assumed to follow an isotropic Gaussian distribution in a four-dimensional space, and the one-sided standard deviation is defined as ⁇ .
  • Fig. 1 (b) shows data obtained by conducting a main factor analysis on the simulated data shown in Fig. 1 (a) and visualizing the top two main factors.
  • the axis corresponding to the first principal factor was PC1
  • the axis corresponding to the second principal factor was PC2.
  • the numbers attached to the respective data sets are the cluster numbers described in FIG. As shown in FIG. 1B, the six clusters can be easily visually confirmed.
  • the principal factor analysis method is not an algorithm for performing clustering, it does not clearly give a clustering result by means other than visual confirmation.
  • the reliability of the clustering result is not given quantitatively.
  • the reliability of the clustering result can be used when comparing cell data obtained from a plurality of subjects. For example, when the clustering result of cell data collected from a healthy subject is reliable, the cell data of other subjects can be verified based on the clustering result. However, if the reliability of the clustering result is unknown, it cannot be determined whether it is appropriate to use this as a reference. Therefore, it is very useful to obtain the reliability of clustering results in cell analysis.
  • FIG. 2 is a diagram for explaining the hierarchical clustering method.
  • the hierarchical clustering method is often used as a method for classifying data in the same manner as the main factor analysis method.
  • the result of clustering the same data as in FIG. 1 by the hierarchical clustering method is shown. Since the data of each cell is expressed as a vector in the four-dimensional Euclidean space, the distance between each cell was evaluated using the Euclidean distance.
  • the optimal number of clusters cannot be obtained.
  • the number of clusters obtained by the hierarchical clustering method does not necessarily match the number of clusters obtained by visual confirmation in the main factor analysis method. Therefore, in the hierarchical clustering method, it is difficult to obtain the optimum number of clusters and the clustering reliability at that time.
  • ⁇ To evaluate the reliability of clustering, it is most natural to fit sample data to a probability distribution.
  • a method of fitting to a probability distribution a method called a mixed Gaussian model is best known.
  • each cell data is assumed to have a Gaussian distribution, and each cell data is considered to belong to any cluster.
  • the log likelihood for the Gaussian probability density function is calculated, and the number of clusters, cluster position (average value), and distribution (standard deviation) are determined by the maximum likelihood method.
  • FIG. 3 is a diagram showing a result of applying the Akaike information criterion to the same simulated sample data as in FIG.
  • various information criterion may be applied to the mixed Gaussian model.
  • the Akaike information criterion is well known as a basic information criterion.
  • clustering algorithms include support vector machine method and k-means method. However, even if these methods are applied to data whose number of clusters is not known in advance, it is difficult to obtain effective results. Even if the optimum number of clusters is obtained by these methods, it is still difficult to evaluate the reliability of the clustering by numerical values.
  • sample data obtained in single cell analysis has a problem that experimental errors cannot be ignored. Since sample data including a large amount of experimental error deviates from the clustering result of sample data with a small error, it is difficult to determine which cluster belongs to, or should be regarded as an independent cluster in the first place. Therefore, it is considered to be important as a cell analysis / diagnosis device to execute clustering with significant resolution on sample data including experimental errors.
  • FIG. 4 is a diagram showing a schematic flow of processing for clustering sample data by the data analysis apparatus according to the first embodiment. Hereinafter, each step shown in FIG. 4 will be described. Details of each step will be described again with reference to FIGS.
  • the data analysis apparatus acquires sample data in the format described in FIG. 1 (S401).
  • the data analysis apparatus acquires data relating to the experimental error when the sample data is acquired (S402). These data may be input to the data analyzer by the operator.
  • the data relating to the experimental error is data that numerically indicates how much the sample data varies due to the experimental error. For example, the standard deviation vector value of each gene data can be used as experimental error data.
  • the data analysis apparatus relaxes the cluster boundary in the subsequent clustering process in order to include the exception data that does not belong to any cluster due to the experimental error in any cluster. Specifically, when any cluster exists at a position away from the exception data in a predetermined range in the clustering space, the exception data is regarded as belonging to the cluster.
  • the predetermined range is called a CR (Clustering Resolution) value in the present invention. Since exceptional data is considered to be caused by experimental error, the CR value is set to a value that is greater than or equal to the experimental error described by the experimental error data. For example, when the experimental error data is represented by the standard deviation ⁇ of the error, the CR value can be a value of about ⁇ to 4 ⁇ .
  • Step S404 The data analyzer performs the following step S405 for each CR value while sweeping the CR value.
  • the CR value sweep range may be set to ⁇ to 4 ⁇ as described above, for example, when the experimental error data is represented by the standard deviation ⁇ of the error.
  • the data analysis apparatus evaluates the optimum number of clusters for each clustering result while temporarily performing clustering by temporarily setting the number k of clusters from 2 to the maximum expected value. Specifically, the likelihood of the current number of clusters is evaluated using the log likelihood of the probability that the sample data belongs to each cluster and the log likelihood of the probability that the sample data does not belong.
  • Step S406 When the likelihood of the number of clusters obtained in step S405 takes an extreme value, the data analysis apparatus considers that the number of clusters is optimum, and adopts the number of clusters as the final number of clusters.
  • Step S407 The data analysis apparatus outputs the clustering result based on the number of clusters determined in step S406 and the reliability of the clustering result. As the reliability of the clustering result, the likelihood value of the number of clusters obtained in step S405 can be used.
  • FIG. 5 is a process flow diagram illustrating the details of step S405. Hereinafter, each step of FIG. 5 will be described.
  • the data analysis apparatus provisionally clusters the sample data for a given temporary cluster number k.
  • the clustering method in this step may be any method such as a hierarchical clustering method or a k-means method.
  • Step S502 The data analysis apparatus replicates k data sets obtained by clustering sample data. Each replicated data set is used to evaluate the likelihood of each clustering result in the following steps.
  • the i-th duplicate data set only the i-th cluster is left, all other data are assumed to belong outside the i-th cluster, and clusters other than the i-th cluster are deleted. In other words, it is assumed that all data not belonging to the i-th cluster belongs to a single cluster.
  • Step S504 The data analysis device determines whether or not the i-th cluster is configured by exception data. An example of exception data will be described later with reference to FIG. If it is not exception data, steps S505 to S507 are executed, and if it is exception data, steps S508 to S509 are executed.
  • Step S504 Judgment Criteria No. 1 of Whether Exception Data or Not
  • the i-th cluster holds a sufficient number of sample data, it is determined that the cluster is not due to exception data.
  • the data structure determined using the correlation matrix calculated from the sample data belonging to the cluster is recognized as having high reliability.
  • a threshold th for determining whether the number of sample data is sufficient is determined in advance. For example, a determination criterion such as considering that the number of sample data in the cluster is 2 or less is regarded as exceptional data can be considered.
  • the threshold th can be selected at random within a certain range, for example. Alternatively, an appropriate probability distribution can be assumed and the probability distribution can be selected at random. In this case, the number of clusters located at the center of the probability distribution is most easily selected.
  • the parameter of the probability distribution can be arbitrarily determined, or a desired probability distribution can be determined by optimization calculation.
  • the data analysis apparatus determines the probability distribution of the sample data using the distribution of the sample data belonging to the i-th cluster.
  • the data analysis apparatus evaluates the validity of the judgment whether each sample data belongs to the i-th cluster or not and whether or not the judgment is correct using the probability distribution.
  • the specific method is as follows.
  • Steps S505 to S507 The data analysis apparatus calculates the standard deviation of the cluster center (average of sample data belonging to the cluster) of the i-th cluster and the sample data in the cluster, and normalizes the sample data (S505).
  • the data analysis apparatus calculates an inverse matrix of the correlation matrix of the sample data in the i-th cluster (S506).
  • the data analyzer calculates the Mahalanobis distance between each sample data and the i-th cluster center (S507). The reason for calculating the Mahalanobis distance from the cluster center will be described with reference to FIGS.
  • Steps S508 to S509 The data analysis apparatus calculates the cluster center of the i-th cluster and normalizes the sample data using the cluster center and the CR value (S508).
  • the data analyzer calculates the Euclidean distance between each sample data and the i-th cluster center (S509). The reason for calculating the Euclidean distance from the cluster center will be described with reference to FIGS.
  • Step S510 Step 1
  • the data analysis apparatus increases the probability that the sample data does not belong to the i-th cluster as the distance from the cluster center increases. Is used to calculate the probability that the sample data does not belong to the i-th cluster.
  • the probability distribution function is such that the probability that the sample data belongs to the i-th cluster decreases as the distance from the cluster center increases. The probability that the sample data belongs to the i-th cluster is calculated.
  • a probability value is calculated according to a cumulative probability distribution function of an x square distribution with n degrees of freedom
  • a probability value is calculated according to a function obtained by subtracting the cumulative probability distribution function from 1.
  • An example of this function is illustrated in FIGS. 6 to 7 described later.
  • Step S510 Step 2
  • the data analysis device calculates the likelihood that each sample data belongs to the i-th cluster or the likelihood that it does not belong by calculating the log likelihood of the probability value calculated according to the above function.
  • the likelihood of the current number of clusters k is calculated by calculating the sum of the log likelihoods for all sample data and all clusters and dividing by the number of clusters k. After this step, the same processing is performed from step S503 on the (i + 1) th cluster.
  • Step S510 Supplement
  • Step S511 Optimization Loop and Cluster Number Sweep Loop
  • the data analysis apparatus repeatedly performs clustering using an optimization method such as the Monte Carlo method so that the log likelihood calculated in step S510 is reduced, thereby obtaining an optimal clustering result and an optimal number of clusters. decide.
  • an optimization method such as the Monte Carlo method
  • the same processing is performed from step S503 while randomly replacing sample data belonging to each cluster.
  • the termination condition for the optimization loop may be, for example, when the likelihood of the current number of clusters k calculated in step S510 reaches a predetermined threshold. After completing the optimization loop, the number k of clusters is incremented and the process returns to step S501 to perform the same processing.
  • Step S511 CR value sweep loop
  • FIG. 6 is a diagram showing a processing image of steps S505 to S507.
  • FIG. 1B a clustering space in which axes corresponding to two principal factors are set is illustrated.
  • FIG. 6 since a relatively large amount of sample data is collected at the center of the cluster, it is assumed that the data included in the cluster is not exceptional data affected by experimental errors.
  • the cluster shape is not necessarily circular in the clustering space
  • linear conversion is performed from the data distribution shown on the left of FIG. 6A to the data distribution shown on the right, and each sample data from the cluster center in the cluster shape after conversion is performed.
  • Find the distance to. corresponds to obtaining the Mahalanobis distance from the center of each data cluster. In other words, this corresponds to the assumption that the cluster forms a unit space in the Mahalanobis Taguchi method.
  • Step S505 In ⁇ S507, for a sample data provisionally determined not to belong to the i-th cluster (e.g. the upper left of the data group in FIG. 6 (a)), with x 2 cumulative probability distribution shown in FIG. 6 (b) Then, the probability that the sample data does not belong to the i-th cluster is calculated. as belonging to the i-th cluster for temporary decision sample data (a data group within a circle in FIG. 6 (a)), using 1-x 2 cumulative probability distribution shown in FIG. 6 (b), the sample data The probability of belonging to the i-th cluster is calculated.
  • the probability value increases as it is closer to the cluster center, and for sample data that is provisionally determined not to belong to the i-th cluster, the further away from the cluster center. Probability value increases.
  • FIG. 7 is a diagram showing a processing image of steps S508 to S509.
  • a clustering space in which axes corresponding to two principal factors are set is illustrated.
  • the number of data belonging to the central cluster is small (two in exception 1 and one in exception 2), it is tentatively determined as exception data. This corresponds to the case where the number of data belonging to the i-th cluster is small and the data structure in the cluster cannot be determined sufficiently.
  • the CR value is the size of the cluster.
  • a predetermined number (exception 1)
  • the sample data that does not belong to the exception cluster has a larger probability value.
  • the likelihood that the exception data is an independent cluster is highly evaluated.
  • the CR value in the present invention has a function of returning sample data deviated from the cluster due to experimental error to the original cluster.
  • the CR value is preferably set to a value larger than the experimental error (because the error cannot be corrected if it is smaller than the experimental error), and is set within a limit range obtained from other medical and biological knowledge.
  • the data analysis apparatus is tentatively determined not to belong to any cluster using the cluster range parameter (CR value) that relaxes the cluster boundary determined based on the experimental error. It is determined whether the exception data belongs to another cluster. As a result, clustering can be performed with high accuracy even when exceptional data resulting from experimental error occurs.
  • CR value cluster range parameter
  • the data analysis apparatus evaluates the likelihood of the clustering result while sweeping the CR value with reference to the Mahalanobis distance or the Euclidean distance from the cluster center, and the number of clusters in which the value takes an extreme value Is considered optimal. Thereby, even if the optimum number of clusters is not known in advance, a good clustering result can be obtained.
  • the data analysis apparatus outputs the reliability of the clustering result together with the clustering result.
  • the diagnostic accuracy by the number of cells belonging to a specific type and the gene expression marker belonging to that type can be improved. That is, conventionally, in a method other than single cell analysis, medical evaluation was performed on a plurality of types of cells.
  • a biomarker for a specific population is set. The evaluation is expected to improve the accuracy of determination of the type and condition of a disease using a biomarker.
  • the clusters determined by the data analysis apparatus according to the first embodiment correspond to cell types derived from the sample data, respectively.
  • a data analysis apparatus that can determine the optimum number of clusters is considered to have the ability to clearly indicate the boundary between one cell type and another cell type. Therefore, the number of cell types can be output from the sample data obtained by cell analysis and diagnosis.
  • a population for outputting a statistic such as an average value of a biomarker typified by a specific gene expression level can be determined. Furthermore, since the reliability is given to the determination, it is easy to determine that data having a certain reliability or less is not adopted.
  • the data analysis apparatus can be applied to an analysis / diagnosis apparatus in the bio / medical field for analyzing individual cell characteristics one by one.
  • the present invention can be applied to an analysis / diagnosis device for analyzing blood cells, an analysis / diagnosis device for cells in urine, or an analysis / diagnosis device for tissue sections. The same applies to the following embodiments.
  • a cDNA library in a single cell is prepared on a magnetic bead based on the method described in Non-Patent Document 1, and a single cell A method of quantifying mRNA as small as 0.5 pg in cells with a qPCR apparatus was used.
  • FIG. 9 is a functional block diagram of the data analysis apparatus 906 according to the second embodiment.
  • the data analysis device 906 includes a calculation unit 904, a data input / output unit 905, a sample data input unit 908, and an experimental error data input unit 909.
  • the measurement system for measuring each data is also shown in FIG.
  • FIG. 9 shows the configuration of the measurement system in the case where the operation of acquiring experimental error data and the single cell analysis to be clustered are performed in parallel.
  • the functional block 901 is a functional unit that acquires sample data to be clustered.
  • the functional block 902 is a functional unit that acquires experimental errors.
  • a functional block 903 is a functional unit that performs processing common to both. In order to accurately evaluate the experimental error, it is desirable that the data for evaluating the experimental error and the sample data to be clustered have many parts in common. When all or a part of the data related to the experimental error is acquired in advance, the data may be stored in the database 907 and read out by the calculation unit 904 when clustering is performed.
  • the functional block 901 collects cells one by one, dispenses each cell into a reaction container, and beads containing poly T probe and a lysis buffer for extracting mRNA into the reaction container into which the cells have been dispensed. Introduce.
  • the functional block 902 introduces a large number of cells into the reaction vessel, and then performs mRNA extraction in the same manner as the functional block 901. Thereafter, the mRNA solution is diluted, and an equivalent amount of a single cell is collected. Dispense this solution into the solution.
  • the functional block 903 removes the lysis buffer, dispenses a reaction solution containing reverse transcriptase, removes the reaction solution after the reverse transcription reaction, adds mRNA-degrading enzyme, and after washing, introduces qPCR reagent, PCR Quantification is performed by performing fluorescence measurements while applying a temperature cycle for.
  • the gene expression level in each sample is calculated based on the cycle number Ct value in which the fluorescence intensity exceeds a certain threshold.
  • the Ct value is read as the number of molecules.
  • Experimental error includes errors in all processes other than single cell sampling until quantification is performed. Since it is desirable that the quantitative value follows a Gaussian distribution, a value obtained by taking the logarithm of the number of molecules is output to the calculation unit 904 as sample data. The same applies to the experimental error data.
  • the sample data input unit 908 acquires sample data via the function blocks 901 and 903.
  • the experimental error data input unit 909 acquires experimental error data via the function blocks 902 and 903.
  • the calculation unit 904 performs clustering using these data according to the processing flow described in the first embodiment.
  • the data input / output unit 905 controls the input / output of these data.
  • steps S508 to S509 are performed when the number of elements in the temporary cluster is smaller than 3 (2 or less).
  • the expression level (sample data value) of a specific gene is known to fluctuate biologically or medically
  • the fluctuation range may be used as the CR value.
  • Such a correspondence relationship between a gene and a CR value and a correspondence relationship between a category of sample data and a CR value are stored in the database 907, and the database 907 is selected from the gene and sample data categories input to the data analysis device 906. If there is a corresponding item that matches the one stored in, the calculation unit 906 reads it from the database 907.
  • the determination threshold value in step S504 may be stored in advance in the database 907 instead of the experimental error data, and an appropriate value may be used based on the sample data or gene information input to the data analysis device 904.
  • the determination threshold value in step S504 is a fixed value corresponding to the number of sample data in the cluster. However, a plurality of determination threshold values are provided and clustering is performed using each threshold value. A threshold having the lowest likelihood (most likely) may be selected. Furthermore, the threshold value may be selected randomly, or may be selected randomly on the assumption that the threshold value follows a certain probability distribution. A plurality of types of probability distribution functions at this time may be stored in the database 907 in advance, and an appropriate probability distribution function may be selected according to the contents of the sample data.
  • ⁇ Embodiment 2 Sampling result>
  • 92 mesenchymal stem cells (C3H10T1 / 2) of mice were used. That is, such cells were collected in order to know the differentiation induction state of the bone of the mouse from which the cells were obtained.
  • other cells for example, immune cells in blood, tissue sections of cancer, etc.
  • the gene expression data used as sample data was quantitative values for five types of genes: Bglap1, Col1a1, Pparg, Col2a1, and Eef1g. The kind and number of genes are only one example, and all possible genes may be measured.
  • mRNA extracted from a large number of cells of about 5 ⁇ 10 5 or more was sampled for 96 pieces of 0.5 pg which is equivalent to a single cell, and 5 ⁇ 10 5 cells
  • the qPCR apparatus was used for quantification.
  • the variation in the quantification value is the total amount of experimental errors composed of handling errors when preparing the cDNA library and quantification errors during qPCR quantification. This error varies from gene to gene. This error was evaluated as a five-dimensional vector.
  • CR value (vector value) is assumed to be ⁇ .
  • the experimental error matches the minimum CR value.
  • the expression of a specific gene changes with time even in the same cell state, and there are cases where the fluctuation range of data is known biologically in addition to the variation in data due to experimental errors. If it is not desired to use these variations in clustering, the value of ⁇ may be partially changed based on the sample data. Specifically, for example, when there is a fluctuation of about 10 times with respect to a specific gene, the numerical value of the fluctuation range may be used as the CR value instead of the experimental error only for this gene. However, this is limited to cases where the experimental error is sufficiently smaller than biological or medical fluctuations.
  • FIG. 10 shows the result of calculating the log likelihood representing the likelihood of the number of clusters in the second embodiment.
  • FIG. 11 is a diagram comparing the result of hierarchical clustering with the result of clustering by the data analysis apparatus 906 according to the second embodiment. Fifteen clusters are represented by square boxes below the tournament diagram corresponding to the hierarchical clustering result. The 15 clusters were found to be composed of 5 significant clusters and 10 exception clusters. Five significant clusters are consistent with the result of hierarchical clustering. It can be seen that the boundary between the significant cluster and the exceptional cluster has been clarified by the second embodiment.
  • FIG. 12 is a diagram showing gene expression levels of Bglap1, PPar, and Col2a1 in each cluster (indicated by numbers 1 to 5).
  • the vertical axis represents the number of molecules of each gene in a single cell.
  • the horizontal axis indicates the cluster number.
  • the black bar is the differentiation inducer, and the gray bar is the gene expression level for the corresponding cells without the differentiation inducer.
  • the distinction between presence and absence of differentiation inducers is an example of cell information that is known in advance.
  • FIG. 13 is a diagram showing the expression levels of the three genes described in FIG. 12 with and without clustering, with and without differentiation induction. As shown in FIG. 13, there is a change corresponding to the prior information for Blap1, but there is no change in the expression level corresponding to the prior information for Pparg and Col2a1, and no information is obtained for these two genes. From the above, it can be seen that more detailed information on gene expression in the living body can be obtained by clustering.
  • the data analysis apparatus 906 is an apparatus for estimating the state of a living body by estimating how many cells belonging to what kind (cluster) exist in the living body. is there.
  • the second embodiment is effective when the type of cluster and the number of cells belonging to the cluster change when the health state of the living body to be measured changes.
  • FIG. 14 is a functional block diagram of the data analysis apparatus 906 according to the third embodiment of the present invention.
  • a configuration example using a large-scale DNA sequencer will be described as a gene expression quantification method. Thereby, the number of genes to be measured can be greatly increased.
  • the functional block 901 extracts the clustering target samples by dispensing them one by one into a reaction vessel containing poly T probe beads on a plate, crushing cells in the reaction vessel, and trapping mRNA on the bead surface. .
  • the function block 902 dispenses known amounts of mRNA for each of a plurality of genes into individual reaction vessels in order to measure experimental errors. By calculating the standard variation deviation from the sequencing data corresponding to the reaction vessel into which a known amount of RNA has been introduced, the experimental error relating to sample processing and sequencing can be quantified.
  • Function block 903 performs reverse transcription reaction and mRNA degradation.
  • a primer for batch amplification is introduced into the end of the cDNA, and then batch amplification is performed using this primer. Further, fragmentation is performed, and amplification is performed using an amplification primer having a cell recognition tag that is different in sequence for each reaction tank (primer sequence is common to all containers) to construct a sequencing library. Since a tag having a different sequence is inserted for each container, that is, for each cell, at the end of the sequencing library, the sample can be mixed in the following processing. In order to sequence the mixed samples with a large-scale sequencer, sequencing is performed using individual amplification such as emulsion PCR or bridge PCR.
  • any device such as a device using fluorescence measurement, a device using FET, or a device using nanopores may be used.
  • mapping the sequencing data obtained from the fragmented sample to a known gene sequence it is determined which sequence of which region of which gene can be measured. Thereafter, the data is aggregated for each gene, and the expression level data for each gene is calculated.
  • an algorithm generally used by the expert may be used. As a result, expression level data for each cell / gene is obtained.
  • the clustering procedure is the same as in Embodiments 1 and 2, but since the number of measured genes is as large as about tens of thousands, cells can be classified in more detail.
  • the data input to the data analysis device 906 is the number of mutations obtained by dividing the genome into a plurality of regions and counting each divided region.
  • the measurement purpose in this case may be, for example, to elucidate the mechanism of cancer occurrence or metastasis by measuring cells from a tissue section of cancer, or to diagnose for selecting a molecular target drug.
  • Supplementary data obtained by analyzing genome data The whole genome or a part thereof is sequence-analyzed for each single cell (assuming that data in mRNA sequence analysis is used), for example, a region is divided every 50 k bases, and mutations are measured.
  • the measurement object includes, for example, single base substitution, deletion, insertion, gene copy number abnormality, and the like.
  • the input data is the number of each mutation. Experimental error can be assessed by artificially creating a sample without mutation and evaluating it.
  • Embodiment 4 of the present invention instead of characterizing a cell by quantifying the expression level (number of molecules) of a gene in the cell, a cell sample obtained by immunostaining or the like is imaged with a fluorescence microscope, and the fluorescence intensity is measured.
  • a fluorescence microscope instead of characterizing a cell by quantifying the expression level (number of molecules) of a gene in the cell, a cell sample obtained by immunostaining or the like is imaged with a fluorescence microscope, and the fluorescence intensity is measured.
  • An example of a configuration in which cells are classified by quantifying the amount of protein in the cells by the correspondence data of the number of molecules or single molecule fluorescence counting will be described.
  • the type of gene corresponds to the type of protein.
  • the amount of protein (number of molecules) for each cell is input to the data analyzer 906 as sample data.
  • sample preparation such as immunostaining
  • an error at the time of fluorescence measurement is evaluated and input to the data analyzer 906 as an experimental error. Thereby, clustering of cells can be executed.
  • FIG. 15 is a functional block diagram of the data analysis apparatus 906 according to the fourth embodiment.
  • the functional block 901 collects a tissue section as a cell sample, and performs immunostaining after deparaffinization. At this time, a known amount of a fluorescent dye having a fluorescence wavelength different from that of immunostaining is introduced into the cell, and at the same time, nuclear staining and cell membrane staining are performed. By carrying out the measurement of the fluorescence of 7 colors, images corresponding to the target protein (3 types), the experimental error measurement dye, the nucleus, and the cell membrane are obtained. When increasing the type of target protein, it is necessary to increase the measurement color.
  • the present invention is not limited to the above-described embodiment, and includes various modifications.
  • the above embodiment has been described in detail for easy understanding of the present invention, and is not necessarily limited to the one having all the configurations described.
  • a part of the configuration of one embodiment can be replaced with the configuration of another embodiment.
  • the configuration of another embodiment can be added to the configuration of a certain embodiment. Further, with respect to a part of the configuration of each embodiment, another configuration can be added, deleted, or replaced.
  • the above components, functions, processing units, processing means, etc. may be realized in hardware by designing some or all of them, for example, with an integrated circuit.
  • Each of the above-described configurations, functions, and the like may be realized by software by interpreting and executing a program that realizes each function by the processor.
  • Information such as programs, tables, and files for realizing each function can be stored in a recording device such as a memory, a hard disk, an SSD (Solid State Drive), or a recording medium such as an IC card, an SD card, or a DVD.
  • 906 data analysis device
  • 904 calculation unit
  • 905 data input / output unit
  • 908 sample data input unit
  • 909 experimental error data input unit.

Abstract

　本発明は、実験誤差などに起因する例外データが生じる場合においても、クラスタリングを適切に実施することができるデータ解析装置を提供することを目的とする。　本発明に係るデータ解析装置は、実験誤差データが記述している実験誤差の範囲に応じて、クラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ定めておく。クラスタリングの過程において、いずれのクラスタにも属さない例外データについては、例外データからさらにクラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、例外データがそのクラスタに属するものと判定し、前記距離だけ離れてなおいずれのクラスタにも含まれない場合は、独立したクラスタを構成するものと判定する。

Description

データ解析装置、データ解析方法
　本発明は、サンプルデータをクラスタリング解析する技術に関するものである。
　再生医療、遺伝子を用いた診断、あるいは生命現象の基本的な理解のため、組織の平均としての遺伝子の発現量ばかりでなく、組織を構成する１つ１つの細胞の中味を定量的に分析することが重要視され始めている。このような個々の細胞の特性を１つずつを解析することを単一細胞解析と呼んでいる。
　単一細胞中の生体分子の量が極めて微量であるため、単一細胞解析は、これまで細胞膜上の蛋白質など一部の生体分子を対象とする解析にしか用いられてこなかった。しかし近年では、技術の発展により単一細胞中の微量な生体分子を定量評価できるようになりつつある。
　下記非特許文献１には、単一細胞中の遺伝子発現解析に関して、定量ＰＣＲ装置を用いた方法によって、十分な精度で特定の遺伝子発現量を計測できる方法が開示されている。下記非特許文献２には、同様に単一細胞中に遺伝子発現解析に関して、大規模ＤＮＡシーケンサ（次世代シーケンサ）を用いて、ほぼすべての遺伝子の発現を定量化する方法が開示されている。非特許文献２には、細胞の種類を特定するためのデータ解析方法も開示されている。今後、ゲノム配列、細胞内の蛋白質、さらに細胞内の様々な生体分子が単一細胞レベルで同定されるようになることが予想される。
Ｎａｔｕｒｅ　Ｍｅｔｈｏｄ，Ｖｏｌ．６，Ｎｏ．７（２００９），ｐｐ．５０３ Ｇｅｎｏｍｅ　Ｒｅｓｅａｒｃｈ，Ｖｏｌ．２１，Ｎｏ．７（２０１１），ｐｐ．１０８８
　上記のような単一細胞解析技術の進歩によって、これまで均一であると仮定され解析されてきた細胞組織は、単一細胞解析で得られるデータを用いて、これまで知られているよりも詳細なグループ、下位組織を構成していることが明らかとなるであろう。その結果、非常に多数の細胞から構成されている人体などの個体中の複雑な生命現象は、細胞のデータによって分類された細胞のグループで構成され、このグループ間を様々な生化学的なシグナルがやりとりされるネットワークとして生命が把握され、生命科学分野、特に医療あるいは創薬分野において大きなインパクトを与えると考えられる。
　たとえば、これまで均一であると考えられてきたがん組織をグループ分けし、各グループに対して遺伝子変異を解析することによって、より適切な分子診断薬を選択することができるようになる可能性がある。また、血中の免疫細胞中の遺伝子発現量を解析することにより、様々な病気を診断することができる可能性が示唆されているが、免疫細胞の詳細な分類によって、より精度の高い診断ができるようになると考えられる。
　しかし、データのみを用いて細胞を分類するアルゴリズムおよびこれを応用した解析／診断装置は、細胞を分類して医療向けの診断において用いるために必要な特性が必ずしも十分ではない。ここでいう必要な特性の例としては、例えば細胞のグループ化（以下ではクラスタリングと呼ぶ）において、最適なグループ数（以下ではクラスタ数と呼ぶ）が事前に分かっていない場合であっても、適切なクラスタリングを実施することなどが考えられる。特に、データ数の少ない例外的なクラスタを独立のクラスタとするのか、他のデータ数の多いクラスタの一部とするのかを判定することは、従来の解析／診断装置においては困難である。
　本発明は、上記のような課題に鑑みてなされたものであり、実験誤差などに起因する例外データが生じる場合においても、クラスタリングを適切に実施することができるデータ解析装置を提供することを目的とする。
　本発明に係るデータ解析装置は、実験誤差データが記述している実験誤差の範囲に応じて、クラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ定めておく。クラスタリングの過程において、いずれのクラスタにも属さない例外データについては、例外データからさらにクラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、例外データがそのクラスタに属するものと判定し、前記距離だけ離れてなおいずれのクラスタにも含まれない場合は、独立したクラスタを構成するものと判定する。
　本発明に係るデータ解析装置によれば、単一細胞の解析結果を用いて細胞を適切に分類することができる。また、分類された細胞の種類数を精度よく決定することができる。
模擬サンプルデータのクラスタリング結果を示す図である。 階層的クラスタリング法について説明する図である。 赤池情報量基準を図１と同じ模擬サンプルデータに対して適用した結果を示す図である。 実施形態１に係るデータ解析装置がサンプルデータをクラスタリングする処理の概略フローを示す図である。 ステップＳ４０５の詳細を説明する処理フロー図である。 ステップＳ５０５～Ｓ５０７の処理イメージを示す図である。 ステップＳ５０８～Ｓ５０９の処理イメージを示す図である。 図１で説明した模擬サンプルデータについてＣＲ値を１、２、４に掃引し、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。 実施形態２に係るデータ解析装置９０６の機能ブロック図である。 実施形態２において、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。 階層的クラスタリングの結果と本実施形態２に係るデータ解析装置９０６によるクラスタリング結果を比較する図である。 各クラスタ（１～５の番号で表示）のＢｇｌａｐ１、ＰＰａｒｇ, Ｃｏｌ２ａ１の遺伝子発現量を示す図である。 図１２で説明した３つの遺伝子について、クラスタリングを実施せず、分化誘導ありの場合となしの場合について発現量を示した図である。 実施形態３に係るデータ解析装置９０６の機能ブロック図である。 実施形態４に係るデータ解析装置９０６の機能ブロック図である。
＜従来のデータ解析手法＞
　以下では本発明の理解を促進するため、まず初めに従来のデータ解析手法およびその課題について、具体例を挙げて説明する。その後に、本発明に係るデータ解析装置の具体的な構成について説明する。
　単一細胞解析で得られた細胞データは、たとえば主因子分析法を用いて解析することができる。主因子分析法によって得られたデータは、目視によってグループを判定するために用いられることが多い。従来のクラスタリング方法の課題を具体的に見るため、以下では模擬サンプルデータを用いることにする。具体的には、４つの遺伝子、１８０個の細胞について単一細胞解析を実施したことを想定し、コンピュータ上で乱数を用いて細胞データを生成した。
　図１は、模擬サンプルデータのクラスタリング結果を示す図である。図１（ａ）は模擬データ内の各クラスタの構成を示す。εは乱数であり、平均値０標準偏差σの分布に従うものとした。クラス２とクラスタ３の間の距離を制御するためのパラメータαを設定し、α＝６×σに設定した。乱数εは４次元空間において等方的なガウス分布に従うと仮定しており、その片側標準偏差を上記σとした。
　図１（ｂ）は、図１（ａ）に示す模擬データについて主因子分析を実施し、上位２個の主因子について可視化したデータである。第１主因子に対応する軸をＰＣ１とし、第２主因子に対応する軸をＰＣ２とした。各データ集合に添えている数字は、図１（ａ）に記載している各クラスタ番号である。図１（ｂ）に示すように、６つのクラスタを容易に目視確認することができる。しかし、主因子分析法はクラスタリングを実施するためのアルゴリズムではないため、目視確認以外の手段によって明確にクラスタリング結果を与えるものではない。また、クラスタリング結果の信頼度も定量的に与えられるものではない。
　クラスタリング結果の信頼度は、複数の被験者から得た細胞データを比較する際に用いることができる。例えば、健常な被験者から採取した細胞データのクラスタリング結果が信頼できる場合、そのクラスタリング結果を基準として他の被験者の細胞データを検証することができる。しかしクラスタリング結果の信頼度が不明であれば、これを基準として用いることが適切であるのか否か判断することができない。したがって、細胞解析においてクラスタリング結果の信頼度を得ることは、非常に有用である。
　図２は、階層的クラスタリング法について説明する図である。階層的クラスタリング法は、主因子分析法と同様にデータを分類する手法としてよく用いられる。ここでは、図１と同じデータを階層的クラスタリング法によってクラスタリングした結果を示す。各細胞のデータは４次元ユークリッド空間内のベクトルとして表現されているため、各細胞間の距離はユークリッド距離を用いて評価した。
　階層的クラスタリング法においては、データ間の距離が最も近い２つのデータをペアにして、その距離に応じた高さのトーナメント図に現れる矩形を設定する。以下、ペアの平均値の位置に２つのデータを代表するデータがあるものと仮定して次のデータに進み、すべてのデータがペアとして結びつけられるまで同様の処理を繰り返す。トーナメント図における矩形の高さが高いほど、クラスタ間の距離が大きいことを示す。ある高さにおいてトーナメント図を水平方向に切断した時に得られる長い鉛直線との交点の数がクラスタ数に相当すると考えられる。
　しかし、階層的クラスタリング法においては、最適なクラスタ数は求めることができない。図２に示す例においては、５個のクラスタあるいは６個のクラスタいずれも妥当なクラスタ数であるように思われる。すなわち、階層的クラスタリング法によって得られるクラスタ数は、主因子分析法において目視確認により得られるクラスタ数とは必ずしも一致しない。したがって、階層的クラスタリング法は、最適なクラスタ数およびその時のクラスタリングの信頼度を求めることは困難である。
　クラスタリングの信頼度を評価するためには、サンプルデータを確率分布へフィッティングさせることが最も自然である。確率分布へフィッティングさせる手法としては、混合ガウシアンモデルと呼ばれる方法が最もよく知られている。混合ガウシアンモデルにおいては、各細胞データはガウス分布していると仮定し、各細胞データはいずれかのクラスタに所属していると見なす。次に、ガウシアン確率密度関数に対する対数尤度を計算し、最尤法によってクラスタ数、クラスタ位置（平均値）、分布（標準偏差）を決定する。
　一般に、単純に対数尤度を用いてクラスタ数などを決定する場合においては、クラスタ数とデータ数が一致するまでクラスタ総数を増やせば良いフィッティングが得られる。したがって、単一細胞解析のように最適なクラスタ総数を決定したい場合においては不向きである。
　図３は、赤池情報量基準を図１と同じ模擬サンプルデータに対して適用した結果を示す図である。上記のような混合ガウシアンモデルにおける問題を回避するため、各種情報量規準を混合ガウシアンモデルに対して適用することが考えられる。赤池情報量基準は、基本的な情報量基準としてよく知られている。
　赤池情報量基準を最適化において用いる場合は、多くの情報量を与えて最適化を実施するほど多くのペナルティ値を課す。これを混合ガウシアンモデルに対して適用する場合、クラスタ数を増やすとその分だけペナルティを課すので、必ずしもクラスタ数を増やすほど良いフィッティングが得られるとは限らないことになる。したがって、評価値が極値を有するクラスタ数が最適であろうと推測することができる。
　ＭａｔＬａｂ２００８ａが搭載しているＥＭアルゴリズムを用いて、赤池情報量基準を混合ガウシアンモデルに対して適用した最適化計算を実施したところ、図３に示すように明確な極値を得ることはできなかった。したがって、混合ガウシアンモデル、および赤池情報量基準をこれに適用した修正手法のいずれも、最適なクラスタ数を決定することは困難であることが分かった。
　その他のクラスタリングアルゴリズムとしては、サポートベクトルマシン法、ｋ－ｍｅａｎｓ法などがある。しかしこれらの手法は、事前にクラスタ数が分からないデータに対して適用しても、有効な結果を得ることは困難である。仮に、これら手法によって最適なクラスタ数が得られたとしても、そのクラスタリングの信頼度を数値によって評価することはやはり困難である。
　事前に情報を与える必要がないデータ解析手法として、多くのデータマイニング手法が知られている。たとえば、非特許文献２において採用されている自己組織化マップ（Ｓｅｌｆ　Ｏｒｇａｎｉｚｉｎｇ　Ｍａｐｓ）がある。しかし、自己組織化マップを用いてクラスタリングを実施したとしても、クラスタリング結果の信頼度を得ることができない。
　上記課題の他にも、単一細胞解析において得られるサンプルデータは実験誤差が無視できないという課題がある。実験誤差を多く含むサンプルデータは、誤差の少ないサンプルデータのクラスタリング結果から外れてしまうため、いずれのクラスタに属するのか、あるいはそもそも独立のクラスタとみなすべきであるのかを判断することが難しい。したがって、実験誤差を含むサンプルデータについて、有意な分解能でクラスタリングを実行することも、細胞解析／診断装置として重要であると考えられる。
＜実施の形態１＞
　以上、従来のデータ解析手法およびその課題について説明した。以下では本発明の実施形態１に係るデータ解析装置について説明する。単一細胞ごとの生体分子に関する解析、特に遺伝子発現解析において得られるデータは、各細胞の生体分子の定量値を要素値とする行列を用いて表現することができる。個々の細胞に対する各遺伝子の発現量のサンプルデータは、遺伝子数をｎ、細胞数をｍとすると、ｍ行ｎ列の行列で表示できる。以下ではこの形式によって記述されたサンプルデータを前提とする。
　図４は、本実施形態１に係るデータ解析装置がサンプルデータをクラスタリングする処理の概略フローを示す図である。以下、図４に示す各ステップを説明する。各ステップの詳細については後述の図５～図６を用いて改めて説明する。
（図４：ステップＳ４０１～Ｓ４０２）
　データ解析装置は、図１で説明した形式のサンプルデータを取得する（Ｓ４０１）。データ解析装置は、サンプルデータを取得したときの実験誤差に関するデータを取得する（Ｓ４０２）。これらのデータは、オペレータがデータ解析装置へ入力してもよい。実験誤差に関するデータとは、サンプルデータが実験誤差に起因してどの程度ばらついているかを数値的に示すデータである。たとえば、各遺伝子データの標準偏差ベクトル値を実験誤差データとすることができる。
（図４：ステップＳ４０３）
　データ解析装置は、実験誤差に起因していずれのクラスタにも属さなくなっている例外データをいずれかのクラスタに包含させるため、以後のクラスタリング処理においてクラスタ境界を緩和する。具体的には、クラスタリング空間において例外データから所定範囲離れた位置にいずれかのクラスタが存在する場合は、その例外データは当該クラスタに属するものとみなす。上記所定範囲のことを、本発明においてはＣＲ（Ｃｌｕｓｔｅｒｉｎｇ　Ｒｅｓｏｌｕｔｉｏｎ）値と呼ぶことにする。例外データは実験誤差によって生じると考えられるため、ＣＲ値は実験誤差データが記述している実験誤差以上の値として設定する。例えば実験誤差データを誤差の標準偏差σによって表す場合、ＣＲ値はσ～４σ程度の値とすることができる。
（図４：ステップＳ４０４）
　データ解析装置は、ＣＲ値を掃引しながら各ＣＲ値について以下のステップＳ４０５を実施する。ＣＲ値の掃引範囲としては、実験誤差データを誤差の標準偏差σによって表す場合、例えば上述のようにσ～４σなどとすればよい。
（図４：ステップＳ４０５）
　データ解析装置は、クラスタ数ｋを２～予想最大値まで仮設定して実際にクラスタリングを実施しながら、各クラスタリング結果についてクラスタ数の最適度を評価する。具体的には、サンプルデータが各クラスタに属する確率の対数尤度および属さない確率の対数尤度を用いて、現在のクラスタ数の尤もらしさを評価する。
（図４：ステップＳ４０６）
　データ解析装置は、ステップＳ４０５において求めたクラスタ数の尤もらしさが極値を取るとき、そのクラスタ数が最適であるものとみなし、そのクラスタ数を最終的なクラスタ数として採用する。
（図４：ステップＳ４０７）
　データ解析装置は、ステップＳ４０６において決定したクラスタ数に基づくクラスタリング結果とともに、そのクラスタリング結果の信頼度を出力する。クラスタリング結果の信頼度としては、ステップＳ４０５において求めたクラスタ数の尤もらしさの値を用いることができる。
　図５は、ステップＳ４０５の詳細を説明する処理フロー図である。以下、図５の各ステップについて説明する。
（図５：ステップＳ５０１）
　データ解析装置は、与えられた仮のクラスタ数ｋについて、サンプルデータを暫定的にクラスタリングする。本ステップにおけるクラスタリング手法は、例えば階層的クラスタリング法やｋ－ｍｅａｎｓ法など任意の手法でよい。
（図５：ステップＳ５０２）
　データ解析装置は、サンプルデータをクラスタリングしたデータセットをｋ個分複製する。複製された各データセットは、それぞれ以下のステップにおいて各クラスタリング結果の尤もらしさを評価するために用いる。データ解析装置は、以下のステップにおいて用いるカウンタｉを初期化する（ｉ＝１）。
（図５：ステップＳ５０３）
　データ解析装置は、ｉ番目（ｉ＝１～ｋ）のクラスタについて、ｉ番目の複製データセットを用いて以下のステップを実施する。ｉ番目の複製データセットに関してはｉ番目のクラスタのみを残し、その他データは全てｉ番目のクラスタ外に属するものと仮定し、ｉ番目のクラスタ以外のクラスタについては削除する。すなわち、ｉ番目のクラスタに属さないデータについては、全て単一のクラスタに属するものと仮定する。
（図５：ステップＳ５０４）
　データ解析装置は、ｉ番目のクラスタが例外データによって構成されているものであるか否かを判定する。例外データの例については後述の図６で説明する。例外データでなければステップＳ５０５～Ｓ５０７を実施し、例外データであればステップＳ５０８～Ｓ５０９を実施する。
（図５：ステップＳ５０４：例外データであるか否かの判定基準その１）
　ｉ番目のクラスタがその内部に十分なサンプルデータ数を保持している場合は、そのクラスタは例外データによるものではないと判定する。この場合は、クラスタに属するサンプルデータから計算される相関行列などを用いて決定したデータ構造は、信頼度が高いと認められるからである。サンプルデータ数が十分であるか否かの閾値ｔｈは、あらかじめ定めておく。例えばクラスタ内のサンプルデータ数が２個以下であれば例外データとみなす、などの判定基準が考えられる。
（図５：ステップＳ５０４：例外データであるか否かの判定基準その２）
　閾値ｔｈは、例えばある範囲内でランダムに選択することもできる。あるいは、適当な確率分布を仮定して、その確率分布を前提としてランダムに選択することもできる。この場合は、確率分布の中央に位置するクラスタ数が最も選択され易くなる。確率分布のパラメータは任意に定めることもできるし、最適化計算によって望ましい確率分布を決定することもできる。
（図５：ステップＳ５０５～Ｓ５０７の総括）
　データ解析装置は、ｉ番目のクラスタ内に属するサンプルデータの分布を用いてそのサンプルデータの確率分布を決定する。データ解析装置は、各サンプルデータがｉ番目のクラスタに属する、あるいは属さないと判定した判断が正しいか否かについて、上記確率分布を用いてその妥当性を評価する。具体的な手法は以下の通りである。
（図５：ステップＳ５０５～Ｓ５０７）
　データ解析装置は、ｉ番目のクラスタのクラスタ中心（クラスタ内に属するサンプルデータの平均）とクラスタ内のサンプルデータの標準偏差を算出し、サンプルデータを規格化する（Ｓ５０５）。データ解析装置は、ｉ番目のクラスタ内のサンプルデータの相関行列の逆行列を計算する（Ｓ５０６）。データ解析装置は、各サンプルデータとｉ番目のクラスタ中心との間のマハラノビス距離を計算する（Ｓ５０７）。クラスタ中心からのマハラノビス距離を計算する理由については後述の図６～図７で説明する。
（図５：ステップＳ５０８～Ｓ５０９）
　データ解析装置は、ｉ番目のクラスタのクラスタ中心を算出し、クラスタ中心とＣＲ値を用いてサンプルデータを規格化する（Ｓ５０８）。データ解析装置は、各サンプルデータとｉ番目のクラスタ中心との間のユークリッド距離を計算する（Ｓ５０９）。クラスタ中心からのユークリッド距離を計算する理由については後述の図６～図７で説明する。
（図５：ステップＳ５１０：ステップ１）
　データ解析装置は、ステップＳ５０１においてｉ番目のクラスタに属さないと判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがｉ番目のクラスタに属さない確率が高くなるような確率分布関数を用いて、当該サンプルデータがｉ番目のクラスタに属さない確率を算出する。同様にステップＳ５０１においてｉ番目のクラスタに属すると判定したサンプルデータについて、クラスタ中心からの距離が離れるほど当該サンプルデータがｉ番目のクラスタに属する確率が低くなるような確率分布関数を用いて、当該サンプルデータがｉ番目のクラスタに属する確率を算出する。例えば、前者については自由度ｎのｘ２乗分布の累積確率分布関数にしたがって確率値を計算し、後者については１から上記累積確率分布関数を引いた関数にしたがって確率値を計算する。この関数の例については後述の図６～図７で例示する。
（図５：ステップＳ５１０：ステップ２）
　データ解析装置は、上記関数にしたがって計算した確率値の対数尤度を計算することにより、各サンプルデータがｉ番目のクラスタに属することの尤もらしさ、または属さないことの尤もらしさを計算する。すべてのサンプルデータおよびすべてのクラスタについてこの対数尤度の和を求め、クラスタ数ｋで除算することにより、現在のクラスタ数ｋの尤もらしさを計算する。本ステップの次は、ｉ＋１番目のクラスタについてステップＳ５０３から同様の処理を実施する。
（図５：ステップＳ５１０：補足）
　対数尤度の評価に用いているのはクラスタリングが妥当である確率の評価値であるため、最適パラメータが得られた時の対数尤度の値から、クラスタリングの信頼度を確率値で出力することが可能である。
（図５：ステップＳ５１１：最適化ループおよびクラスタ数掃引ループ）
　データ解析装置は、例えばモンテカルロ法などのような最適化手法を用いて、ステップＳ５１０で計算した対数尤度が小さくなるようにクラスタリングを繰り返し実施することにより、最適なクラスタリング結果と最適なクラスタ数を決定する。最適化手法として例えばモンテカルロ法を用いる場合は、各クラスタに属するサンプルデータをランダムに入れ替えながらステップＳ５０３から同様の処理を実施する。最適化ループの終了条件は、例えばステップＳ５１０で計算した現在のクラスタ数ｋの尤もらしさが所定閾値に達した時点などとすればよい。最適化ループを完了した後はクラスタ数ｋをインクリメントしてステップＳ５０１に戻り、同様の処理を実施する。
（図５：ステップＳ５１１：ＣＲ値掃引ループ）
　データ解析装置は、現在のＣＲ値について最適化ループおよびクラスタ数掃引ループを終了した後は、ＣＲ値をインクリメントしてステップＳ５０１に戻り、同様の処理を実施する。ＣＲ値をインクリメントする幅は、想定しているＣＲ値の最小値と最大値の差に応じて適宜定める。
　図６は、ステップＳ５０５～Ｓ５０７の処理イメージを示す図である。ここでは図１（ｂ）と同様に、２つの主因子に対応する軸を設定したクラスタリング空間を例示した。図６に示すクラスタリング例においては、クラスタ中心に比較的多くのサンプルデータが集まっているため、当該クラスタ内に含まれるデータは実験誤差の影響を受けた例外データではないと仮定する。
　クラスタの形状は必ずしもクラスタリング空間において円形ではないため、図６（ａ）の左に示すデータ分布から右に示すデータ分布へと線形変換を実行し、変換後のクラスタ形状においてクラスタ中心から各サンプルデータまでの距離を求める。これは、各データのクラスタの中心からのマハラノビス距離を求めることに相当する。換言すると、クラスタがマハラノビス田口法における単位空間を形成していると仮定することに対応する。
　ステップＳ５０５～Ｓ５０７においては、ｉ番目のクラスタに属さないと仮判定したサンプルデータ（例えば図６（ａ）の左上のデータ群）については、図６（ｂ）に示すｘ２累積確率分布を用いて、そのサンプルデータがｉ番目のクラスタに属さない確率を計算する。ｉ番目のクラスタに属すると仮判定したサンプルデータ（図６（ａ）の円内のデータ群）については、図６（ｂ）に示す１－ｘ２累積確率分布を用いて、そのサンプルデータがｉ番目のクラスタに属する確率を計算する。これにより、ｉ番目のクラスタに属すると仮判定したサンプルデータについては、クラスタ中心に近いほど確率値が高くなり、ｉ番目のクラスタに属さないと仮判定したサンプルデータについては、クラスタ中心から離れるほど確率値が高くなる。
　図７は、ステップＳ５０８～Ｓ５０９の処理イメージを示す図である。ここでは図１（ｂ）と同様に、２つの主因子に対応する軸を設定したクラスタリング空間を例示した。図６に示すクラスタリング例においては、中央のクラスタに属するデータ数が少ないため（例外１においては２個、例外２においては１個）、例外データであると仮判定される。これはすなわち、ｉ番目のクラスタに属するデータ数が少ないため、そのクラスタ内のデータ構造を十分に決定することができない場合に相当する。
　ステップＳ５０８～Ｓ５０９においては、ＣＲ値を当該クラスタのサイズであると仮定する。ｉ番目のクラスタの周辺に存在する、同クラスタに属さないと判定されたデータ個数が所定個数未満である場合（例外１）、例外データは独立したクラスタを構成している可能性が高い。この場合、例外クラスタに属さないサンプルデータはその確率値が大きくなる。結果として、例外データが独立したクラスタであることの尤もらしさが高く評価されることになる。
　一方、ｉ番目のクラスタの周辺に、同クラスタに属さないと判定されたデータ個数が所定個数以上存在する場合（例外２）、例外データは本来他のクラスタに属していたが実験誤差によってクラスタから外れた可能性が高い。この場合、近傍のクラスタが例外クラスタに属する確率値が相応に高くなる。結果として、例外データが独立したクラスタであることの尤もらしさが低く評価されることになる。
　図７に示すように、本発明におけるＣＲ値は、実験誤差によってクラスタから外れたサンプルデータを元のクラスタに戻す機能を備える。そのためＣＲ値は、実験誤差よりも大きな値とし（実験誤差より小さいと誤差を修正できないため）、他の医学的、生物学的な知見から得られる制限範囲内で設定することが望ましい。
　図８は、図１で説明した模擬サンプルデータについてＣＲ値を１、２、４（それぞれ４×σ、８×σ、１６×σに対応）に掃引し、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。クラスタ数ｋ＝６において対数尤度が最小値を示し、その最小値はＣＲ値を変えても安定であることが観察される。すなわち、クラスタ数の対数尤度が極値となるクラスタ数が最適であると判定することができる。
＜実施の形態１：まとめ＞
　以上のように、本実施形態１に係るデータ解析装置は、実験誤差に基づき定めた、クラスタ境界を緩和するクラスタ範囲パラメータ（ＣＲ値）を用いて、いずれのクラスタにも属さないと仮判定された例外データが他のクラスタに属するか否かを判定する。これにより、実験誤差に起因する例外データが生じる場合であっても、精度よくクラスタリングを実施することができる。
　また、本実施形態１に係るデータ解析装置は、クラスタ中心からのマハラノビス距離またはユークリッド距離を基準として、ＣＲ値を掃引しながらクラスタリング結果の尤もらしさを評価し、その値が極値を取るクラスタ数を最適とみなす。これにより、最適なクラスタ数があらかじめ分かっていない場合であっても、良好なクラスタリング結果を得ることができる。
　また、本実施形態１に係るデータ解析装置は、クラスタリング結果の信頼度をクラスタリング結果と併せて出力する。これにより、特定種類に属する細胞数やその種類に属する遺伝子発現マーカによる診断精度を向上させることができる。すなわち、従来は単一細胞解析以外の方法においては複数種類の細胞に対して医学的評価を実施していたところ、本実施形態１に係るデータ解析装置のクラスタリング結果に基づき特定集団に対するバイオマーカを評価することによって、バイオマーカによる疾患の種類や状態の判定の精度が向上するものと期待される。
　本実施形態１に係るデータ解析装置によって決定したクラスタは、それぞれサンプルデータから導かれた細胞の種類に対応する。最適なクラスタ数を決定できるデータ解析装置は、ひとつの細胞種と別の細胞種との間の境界を明示できる能力を有すると考えられる。それゆえ、細胞の解析や診断で得られたサンプルデータから細胞種毎の個数を出力することができる。また、個々の細胞種において、特定の遺伝子発現量などに代表されるバイオマーカの平均値などの統計量を出力するための母集団を決定することもできる。さらに、その決定には信頼度が付与されているため、ある信頼度以下のデータを採用しないなどの判定も容易である。
　本実施形態１に係るデータ解析装置は、具体的には、個々の細胞の特性を１細胞ずつ解析するバイオ・医療分野の解析・診断装置において適用することができる。特に、血球を解析対象とする解析／診断装置、または尿中の細胞を対象とする解析／診断装置、または組織切片を対象とする解析／診断装置において適用することができる。以下の実施形態についても同様である。
＜実施の形態２：装置構成＞
　本発明の実施形態２では、実施形態１で説明したデータ解析装置の具体的な適用例として、単一細胞中の遺伝子発現解析によるデータによって細胞を分類するデータ解析装置について説明する。
　本実施形態２においては、単一細胞中の遺伝子発現解析を実施するため、非特許文献１に記載されている方法に基づき、磁性ビーズ上に単一細胞中のｃＤＮＡライブラリを作製し、単一細胞中の０．５ｐｇという微量なｍＲＮＡをｑＰＣＲ装置によって定量する方法を用いた。
　図９は、本実施形態２に係るデータ解析装置９０６の機能ブロック図である。データ解析装置９０６は、演算部９０４、データ入出力部９０５、サンプルデータ入力部９０８、実験誤差データ入力部９０９を備える。各データを測定する測定系についても図９内に併記した。
　図９は、実験誤差データを取得する動作と、クラスタリング対象となる単一細胞解析とを平行して実施する場合における測定系の構成を示した。機能ブロック９０１は、クラスタリング対象となるサンプルデータを取得する機能部である。機能ブロック９０２は、実験誤差を取得する機能部である。機能ブロック９０３は、両者に共通する処理を実施する機能部である。正確な実験誤差を評価するためには、実験誤差を評価するためのデータとクラスタリング対象とするサンプルデータは共通である部分が多い方が望ましい。実験誤差に関する全部または一部のデータを事前に取得している場合は、そのデータをデータベース９０７に保存しておき、クラスタリングを実施するとき演算部９０４がこれを読み出すようにしてもよい。
　機能ブロック９０１は、細胞を１つずつ採取し、個々の細胞ごとに反応容器に分注し、この細胞を分注した反応容器に、ｍＲＮＡを抽出するためのポリＴプローブ付きビーズとリシスバッファを導入する。
　機能ブロック９０２は、多数の細胞を反応容器に導入した後、機能ブロック９０１と同様にｍＲＮＡ抽出を実施し、その後ｍＲＮＡ溶液を希釈し、単一細胞相当量を採取して、ポリＴプローブ付ビーズ溶液にこの溶液を分注する。
　機能ブロック９０３は、リシスバッファを除去し、逆転写酵素を含む反応溶液を分注し、逆転写反応後に反応溶液を除去し、ｍＲＮＡの分解酵素を加えて洗浄後、ｑＰＣＲ試薬を導入し、ＰＣＲのための温度サイクルを印加しながら蛍光測定を実施することによって定量を実施する。
　個々のサンプル中の遺伝子発現量は、蛍光強度がある閾値を超えたサイクル数Ｃｔ値に基づき算出される。分子数が既知の複数のＤＮＡに対してｑＰＣＲ定量を実行することによって、Ｃｔ値は分子数に読み替えられる。実験誤差は、定量を実施するまでにおける単一細胞のサンプリング以外のすべてのプロセス内の誤差を含む。定量値はガウス分布に従っているほうが望ましいので、分子数の対数をとった値をサンプルデータとして演算部９０４に出力する。実験誤差データについても同様である。
　サンプルデータ入力部９０８は、機能ブロック９０１、９０３を介してサンプルデータを取得する。実験誤差データ入力部９０９は、機能ブロック９０２、９０３を介して実験誤差データを取得する。演算部９０４は、これらのデータを用いて、実施形態１で説明した処理フローにしたがってクラスタリングを実施する。データ入出力部９０５は、これらデータの入出力を制御する。
　本実施形態２においては、仮クラスタ内の要素数が３より小さい（２以下）場合に、ステップＳ５０８～Ｓ５０９を実施することとした。生物学的あるいは医学的に特定の遺伝子の発現量（サンプルデータの値）が変動すると知られている場合、その変動幅をＣＲ値として用いてもよい。このような、遺伝子とＣＲ値の対応関係やサンプルデータのカテゴリとＣＲ値の対応関係をデータベース９０７に保存しておき、データ解析装置９０６に入力された遺伝子やサンプルデータのカテゴリのなかでデータベース９０７に保存しておいたものと合致する該当するものがあれば、演算部９０６はそれをデータベース９０７から読み出す。あるいは、実験誤差データに代えてステップＳ５０４における判定閾値をデータベース９０７にあらかじめ格納しておき、データ解析装置９０４に入力されたサンプルデータや遺伝子情報に基づいて適当な値を用いるようにしてもよい。
　本実施形態２では、ステップＳ５０４における判定閾値は、クラスタ内のサンプルデータ数に対応する固定した値としたが、判定閾値を複数設けておき、各閾値をそれぞれ用いてクラスタリングを実施し、もっとも対数尤度の低い（最も尤もらしい）閾値を選択するようにしてもよい。さらには、閾値をランダムに選択するか、あるいは閾値がある確率分布にしたがう前提の下でランダムに選択してもよい。このときの確率分布関数はあらかじめ複数種類をデータベース９０７に格納しておき、サンプルデータの内容に応じて適当な確率分布関数を選択すればよい。
＜実施の形態２：サンプリング結果＞
　以下では、本実施形態２に係るデータ解析装置９０６によるサンプリング結果の例を説明する。サンプルとしてマウスの間葉系幹細胞（Ｃ３Ｈ１０Ｔ１／２）を９２個用いた。すなわち、細胞を取得したマウスの骨の分化誘導状態を知るためにこのような細胞を採取した。もちろん、他の目的のための他の生体（人を含む）から、別の細胞（たとえば、血中の免疫細胞やがんの組織切片など）を採取してもよい。サンプルデータとして用いる遺伝子発現データは、Ｂｇｌａｐ１、Ｃｏｌ１ａ１、Ｐｐａｒｇ、Ｃｏｌ２ａ１、Ｅｅｆ１ｇの５種類の遺伝子についての定量値とした。遺伝子の種類およびその数は１例であり、可能性のある遺伝子すべてを計測対象にしてもかまわない。
　次に、実験誤差を評価するため、５×１０５個程度以上の多数の細胞から抽出したｍＲＮＡを単一細胞相当量である０．５ｐｇずつ９６個分サンプリングし、５×１０５個の細胞の遺伝子発現量の平均値を計測するためｑＰＣＲ装置を用いて定量した。定量値のばらつきはｃＤＮＡライブラリを作製するときのハンドリング誤差や、ｑＰＣＲ定量時の定量誤差から構成される実験誤差の総量である。この誤差は遺伝子ごとに異なる。この誤差を５次元ベクトルとして評価した。
　実験誤差は以下のような方法によって取得することが望ましい。まず、多数の細胞から抽出したｍＲＮＡサンプルを作成するとき、細胞数の異なる複数のサンプルを準備し、これらのサンプルから単一細胞相当量のｍＲＮＡを採取し、ｃＤＮＡを構築後、各細胞数および各遺伝子に関するｑＰＣＲ定量時の誤差を定量する。そして、細胞数無限大の場合の値を外挿によって推定する。実際には、細胞数の逆数を横軸に実験誤差（標準偏差）を縦軸にとって、ｙ切片の値を実験誤差の推定値とする。
　得られた実験誤差に基づき、許容できるＣＲ値を定める。このＣＲ値（ベクトル値）をσとする。本実施形態２においては、実験誤差と最小ＣＲ値が一致するものとする。ただし、特定の遺伝子の発現は、同じ細胞状態であっても時間的に変化しており、実験誤差によるデータのばらつき以外に生物学的にもデータの変動幅が分かっている場合がある。これらの変動をクラスタリングにおいて用いたくない場合は、σの値をサンプルデータに基づいて一部変更しても構わない。具体的には、たとえば、特定の遺伝子について１０倍程度の変動がある場合、この遺伝子のみについて、実験誤差ではなくこの変動幅の数値をＣＲ値として用いてもよい。ただし、実験誤差が生物学的あるいは医学的変動より十分小さい場合に限られる。
　図１０は、本実施形態２において、クラスタ数の尤もらしさを表す対数尤度を計算した結果を示す。実験誤差をσとして、ＣＲ値は１σから２．１σまでについて評価した。すべてのＣＲ値について、クラスタ数ｋ＝１５のとき、安定的に極値が得られた。
　図１１は、階層的クラスタリングの結果と本実施形態２に係るデータ解析装置９０６によるクラスタリング結果を比較する図である。１５個のクラスタは階層的クラスタリング結果に対応するトーナメント図の下に四角の箱で表現している。１５個のクラスタは、５個の有意なクラスタと１０個の例外クラスタから構成されることが分かった。５個の有意なクラスタは階層的クラスタリングの結果と一致する。本実施形態２により、有意なクラスタと例外クラスタの境界が明確になったことが分かる。
　図１２は、各クラスタ（１～５の番号で表示）のＢｇｌａｐ１、ＰＰａｒｇ, Ｃｏｌ２ａ１の遺伝子発現量を示す図である。縦軸は、単一細胞中の各遺伝子の分子数である。横軸はクラスタ番号を示している。黒いバーは分化誘導剤あり、灰色のバーは分化誘導剤なしに対応する細胞についての遺伝子発現量である。分化誘導剤あり・なしの区別は事前に分かっている細胞情報の例である。
　図１２に示すように、分化誘導剤あり・なしは、個々のクラスタのなかでは細胞を区別するための重要な情報でないことが分かる。一方、３つの遺伝子の発現プロファイル、すなわち、（Ｂｇｌａｐ１，Ｐｐａｒｇ，Ｃｏｌ２ａ１）の発現量を＋－で表現すると、クラスタ１は（－，＋，＋）、クラスタ２は（＋，＋，＋）、クラスタ３は（－，－，＋）、クラスタ４は（－，＋，－）、クラスタ５は（＋，－，－）という明確な特徴を持っていることがわかる。同時にクラスタごとの遺伝子発現量の変動量は、上記特徴を識別するには十分でないことがわかる。事前情報との対応は各クラスタに属する細胞数に反映されるので、その意味においても正確なクラスタリングが重要であることがわかる。
　図１３は、図１２で説明した３つの遺伝子について、クラスタリングを実施せず、分化誘導ありの場合となしの場合について発現量を示した図である。図１３に示すようにＢｇｌａｐ１については事前情報に対応した変化が見られるが、ＰｐａｒｇおよびＣｏｌ２ａ１については事前情報に対応した発現量の変化がなく、これら２つの遺伝子に関しては情報が得られない。以上から、クラスタリングによって、より詳細な生体の遺伝子発現に関する情報が得られることが分かる。
　以上のように、本実施形態２によれば、個々の細胞中の遺伝子発現量を計測し、クラスタリングによって、細胞を構成する生体の情報を得ることができる。すなわち本実施形態２に係るデータ解析装置９０６は、生体中にどのような種類（クラスタ）に属する細胞がどの程度の数だけ存在するかを推定することによって生体の状態を推定するための装置である。計測対象となっている生体の健康状態が変化すると、クラスタの種類とクラスタに属する細胞の数が変化する場合において、本実施形態２は有効である。
＜実施の形態３＞
　図１４は、本発明の実施形態３に係るデータ解析装置９０６の機能ブロック図である。本実施形態３では、遺伝子発現の定量方法として、大規模ＤＮＡシーケンサを用いる構成例を説明する。これにより、計測対象とする遺伝子数を大幅に増やすことができる。
　機能ブロック９０１は、クラスタリング対象サンプルをプレート上のポリＴプローブ付ビーズの入った反応容器に１つずつ分注し、反応容器内で細胞を破砕し、ｍＲＮＡをビーズ表面でトラップすることによって抽出する。
　機能ブロック９０２は、実験誤差測定のため、複数の遺伝子について既知量のｍＲＮＡを個々の反応容器に分注する。既知量ＲＮＡを導入した、反応容器に対応するシーケンシングデータからその標準変偏差を算出することによって、サンプル処理とシーケンシングに関する実験誤差を定量することができる。
　機能ブロック９０３は、逆転写反応とｍＲＮＡ分解を実施する。このとき、ｃＤＮＡの末端に一括増幅用プライマを導入し、次にこのプライマを用いて一括増幅を実施する。さらに、断片化し、各反応槽ごとに配列のことなる細胞認識タグをもった増幅プライマ（プライマ配列はすべての容器で共通）を用いて一括増幅を実施し、シーケンシングライブラリを構築する。シーケンシングライブラリの末端には容器ごとすなわち細胞ごとに異なる配列のタグが挿入されているので、以下の処理においてはサンプルを混合することができる。混合したサンプルを大規模シーケンサでシーケンシングするために、エマルジョンＰＣＲやブリッジＰＣＲなどの個別増幅を用いてシーケンシングを実行する。シーケンシングは、蛍光計測を用いる装置、ＦＥＴを用いる装置、ナノポアを用いる装置などいずれの装置を用いてもよい。断片化したサンプルから得られたシーケンシングデータを既知の遺伝子配列にマッピングすることによって、どの遺伝子のどの領域の配列が計測できたかを判定する。その後、データを遺伝子毎に集計し、遺伝子毎の発現量データを算出する。このときの算出アルゴリズムは当該専門家が一般的に用いるアルゴリズムを用いてよい。その結果、細胞毎・遺伝子毎の発現量データが得られる。
　クラスタリング手順については実施形態１～２と同様であるが、計測遺伝子数が数万程度と非常に大きいため、細胞をより詳細に分類することができる。
　同様の処理を用いて、各細胞のゲノムデータを解析することもできる。データ解析装置９０６に入力されるデータとしては、ゲノムを複数の領域に分割し、各分割領域についてカウントした変異数である。この場合の計測目的は、例えば癌の組織切片からの細胞を計測することによって、癌の発生や転移のメカニズムを解明すること、あるいは分子標的薬を選択するための診断、などが考えられる。
　ゲノムデータを解析することによって得られるデータについて補足する。ゲノム全体またはその一部を単一細胞ごとに配列解析し、（ｍＲＮＡの配列解析におけるデータを用いることを想定）たとえば５０ｋ塩基ごとに領域を分けて、変異を計測する。計測対象は、たとえば一塩基置換、欠失、挿入、遺伝子コピー数異常、などである。入力データはそれぞれの変異数である。実験誤差は、変異がないサンプルを人工的に作製し、それを評価することによって、評価することができる。
　ゲノムを直接配列決定するためには、ｍＲＮＡ抽出の替わりにＲＮＡ分解、ＤＮＡ抽出後、断片化、ポリＡテイリングを酵素試薬の添加によって実行する必要がある。その後はｍＲＮＡの処理と同様である。
＜実施の形態４＞
　本発明の実施形態４では、細胞中の遺伝子の発現量（分子数）を定量することによって細胞を特徴付ける代わりに、免疫染色法などで得られた細胞サンプルを蛍光顕微鏡でイメージングして、蛍光強度と分子数の対応データまたは単分子蛍光カウンティングによって細胞中の蛋白質量を定量することにより細胞を分類する構成例を説明する。
　本実施形態４においては、遺伝子の種類が蛋白の種類に対応する。細胞ごとの蛋白量（分子数）を、サンプルデータとしてデータ解析装置９０６に入力する。免疫染色などのサンプル作製時は、蛍光計測時の誤差を評価して実験誤差としてデータ解析装置９０６に入力する。これにより、細胞のクラスタリングを実行することができる。
　図１５は、本実施形態４に係るデータ解析装置９０６の機能ブロック図である。機能ブロック９０１は、細胞サンプルとして組織切片を採取し、脱パラフィン処理後に免疫染色を実施する。このとき、免疫染色と異なる蛍光波長の既知量の蛍光色素を細胞内に導入すると同時に核染色と細胞膜染色を実施する。７色の多色の蛍光計測を実施することによって、ターゲット蛋白（３種類）と実験誤差計測用色素、核、細胞膜に対応する像を得る。ターゲット蛋白の種類を増やす場合には計測色を増やす必要がある。細胞膜と核を認識し、細胞膜中に核が１つある物体を細胞と認識する。個々の細胞と認識された領域のターゲット蛋白に相当する蛍光色の強度、または単分子蛍光時の分子数を定量値として記録する。同時に実験誤差用色素の定量も実行する。実験誤差は、認識した細胞の面積あたりの強度の標準偏差を求め、細胞の平均面積から算出する。
　本発明は上記した実施形態に限定されるものではなく、様々な変形例が含まれる。上記実施形態は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施形態の構成の一部を他の実施形態の構成に置き換えることもできる。また、ある実施形態の構成に他の実施形態の構成を加えることもできる。また、各実施形態の構成の一部について、他の構成を追加・削除・置換することもできる。
　上記各構成、機能、処理部、処理手段等は、それらの一部や全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリ、ハードディスク、ＳＳＤ（Ｓｏｌｉｄ　Ｓｔａｔｅ　Ｄｒｉｖｅ）等の記録装置、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に格納することができる。
　９０６：データ解析装置、９０４：演算部、９０５：データ入出力部、９０８：サンプルデータ入力部、９０９：実験誤差データ入力部。

Claims (11)

  1. 　複数のサンプルデータをクラスタリング解析する装置であって、
    　複数のサンプルデータを受け取るサンプルデータ入力部と、
    　前記サンプルデータの実験誤差についての情報を記述した実験誤差データを受け取る実験誤差データ入力部と、
    　前記複数のサンプルデータをクラスタリング空間内においてクラスタリングする演算部と、
    　前記クラスタリングの結果を出力する出力部と、
    　を備え、
    　前記演算部は、
    　　前記実験誤差データが記述している前記実験誤差の範囲に応じて、前記クラスタリングを実施する際におけるクラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ取得しておき、
    　　仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングし、
    　　前記複数のサンプルデータのうちいずれのクラスタにも属さない例外データについては、
    　　前記クラスタリング空間上で前記例外データからさらに前記クラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、前記例外データがそのクラスタに属するものと判定し、前記領域がいずれのクラスタにも含まれない場合は、前記例外データが独立したクラスタを構成するものと判定する
    　ことを特徴とするデータ解析装置。
  2. 　前記演算部は、
    　　各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属する確率の尤もらしさを表す第１対数尤度と、各前記サンプルデータが前記クラスタリングの結果得られた各クラスタに属さない確率の尤もらしさを表す第２対数尤度とをそれぞれ算出する処理を、前記第１対数尤度と前記第２対数尤度を用いて算出される前記クラスタリング結果の尤もらしさが所定閾値に達するまで繰り返すことにより、最適なクラスタ総数を求め、
    　　得られた最適なクラスタ総数に準じて、前記複数のサンプルデータの最終的なクラスタリング結果を決定する
    　ことを特徴とする請求項１記載のデータ解析装置。
  3. 　前記演算部は、
    　　前記仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングする過程において、前記サンプルデータが仮設定したクラスタに属すると仮定して前記第１対数尤度と前記第２対数尤度を算出し、
    　　前記仮設定したクラスタに属すると仮定した前記サンプルデータについては、前記クラスタリング空間上において、前記仮設定したクラスタの中心からの距離が離れるにしたがって、前記仮設定したクラスタに属する確率を低く評価し、
    　　前記仮設定したクラスタに属さないと仮定した前記サンプルデータについては、前記クラスタリング空間上において、前記仮設定したクラスタの中心からの距離が離れるにしたがって、前記仮設定したクラスタに属さない確率を高く評価する
    　ことを特徴とする請求項２記載のデータ解析装置。
  4. 　前記演算部は、
    　　クラスタ内に属する前記サンプルデータの個数が所定個数以上であるか否かを基準として、前記サンプルデータが前記例外データであるか否かを判定する
    　ことを特徴とする請求項１記載のデータ解析装置。
  5. 　前記演算部は、前記所定個数をランダムに選択する
    　ことを特徴とする請求項４記載のデータ解析装置。
  6. 　前記演算部は、所定の確率分布にしたがって前記所定個数をランダムに選択する
    　ことを特徴とする請求項４記載のデータ解析装置。
  7. 　前記演算部は、
    　　前記クラスタ範囲パラメータを掃引し、各前記クラスタ範囲パラメータの値を採用した場合における前記クラスタリングによって得られたクラスタ総数を取得し、
    　　前記第１対数尤度と前記第２対数尤度に基づき算出した前記クラスタリング結果の尤もらしさの値が極値となるクラスタ総数を、最適なクラスタ総数として採用する
    　ことを特徴とする請求項２記載のデータ解析装置。
  8. 　前記演算部は、
    　　前記クラスタリングを実施する過程において得られる情報を用いて、前記クラスタリングの結果の信頼度指標を算出し、
    　前記出力部は、
    　　前記信頼度指標を前記クラスタリングの結果と併せて出力する
    　ことを特徴とする請求項１記載のデータ解析装置。
  9. 　前記演算部は、
    　　前記第１対数尤度と前記第２対数尤度に基づき算出した前記クラスタリング結果の尤もらしさの値を、前記クラスタリングの結果の信頼度指標として算出する
    　ことを特徴とする請求項８記載のデータ解析装置。
  10. 　前記サンプルデータ入力部と前記実験誤差データ入力部は、細胞の解析結果に関するデータをそれぞれ前記サンプルデータおよび前記実験誤差データとして受け取り、
    　前記演算部は、前記クラスタリングによって前記細胞をグループ化する
    　ことを特徴とする請求項１記載のデータ解析装置。
  11. 　複数のサンプルデータをクラスタリング解析する方法であって、
    　複数のサンプルデータを受け取るサンプルデータ入力ステップ、
    　前記サンプルデータの実験誤差についての情報を記述した実験誤差データを受け取る実験誤差データ入力ステップ、
    　前記複数のサンプルデータをクラスタリング空間内においてクラスタリングする演算ステップ、
    　前記クラスタリングの結果を出力する出力ステップ、
    　を有し、
    　前記演算ステップにおいては、
    　　前記実験誤差データが記述している前記実験誤差の範囲に応じて、前記クラスタリングを実施する際におけるクラスタ境界を緩和するためのクラスタ範囲パラメータをあらかじめ取得しておき、
    　　仮設定したクラスタ総数に準じて前記複数のサンプルデータをクラスタリングし、
    　　前記複数のサンプルデータのうちいずれのクラスタにも属さない例外データについては、
    　　前記クラスタリング空間上で前記例外データからさらに前記クラスタ範囲パラメータによって定められる距離だけ離れた領域がいずれかのクラスタに含まれる場合は、前記例外データがそのクラスタに属するものと判定し、前記領域がいずれのクラスタにも含まれない場合は、前記例外データが独立したクラスタを構成するものと判定する
    　ことを特徴とするデータ解析方法。
PCT/JP2012/080003 2012-11-20 2012-11-20 データ解析装置、データ解析方法 WO2014080447A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US14/443,118 US20150302042A1 (en) 2012-11-20 2012-11-20 Data analysis apparatus and data analysis method
JP2014548349A JP6029683B2 (ja) 2012-11-20 2012-11-20 データ解析装置、データ解析プログラム
PCT/JP2012/080003 WO2014080447A1 (ja) 2012-11-20 2012-11-20 データ解析装置、データ解析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2012/080003 WO2014080447A1 (ja) 2012-11-20 2012-11-20 データ解析装置、データ解析方法

Publications (1)

Publication Number Publication Date
WO2014080447A1 true WO2014080447A1 (ja) 2014-05-30

Family

ID=50775654

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/080003 WO2014080447A1 (ja) 2012-11-20 2012-11-20 データ解析装置、データ解析方法

Country Status (3)

Country Link
US (1) US20150302042A1 (ja)
JP (1) JP6029683B2 (ja)
WO (1) WO2014080447A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016158576A (ja) * 2015-03-03 2016-09-05 富士フイルム株式会社 細胞群検出装置および方法並びにプログラム
JP2018530815A (ja) * 2015-08-17 2018-10-18 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 生体データにおけるパターン認識のマルチレベルアーキテクチャ
WO2022009342A1 (ja) * 2020-07-08 2022-01-13 富士通株式会社 情報処理プログラム、情報処理方法および情報処理装置
US11327470B2 (en) 2017-12-21 2022-05-10 Mitsubishi Power, Ltd. Unit space generating device, plant diagnosing system, unit space generating method, plant diagnosing method, and program

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6274114B2 (ja) * 2013-01-10 2018-02-07 富士通株式会社 制御方法、制御プログラム、および制御装置
EP3361416A1 (en) 2017-02-09 2018-08-15 Tata Consultancy Services Limited Statistical modeling techniques based neural network models for generating intelligence reports
US10846311B2 (en) * 2017-11-01 2020-11-24 Mad Street Den, Inc. Method and system for efficient clustering of combined numeric and qualitative data records
US10635939B2 (en) 2018-07-06 2020-04-28 Capital One Services, Llc System, method, and computer-accessible medium for evaluating multi-dimensional synthetic data using integrated variants analysis
CN109325059A (zh) * 2018-12-03 2019-02-12 枘熠集成电路（上海）有限公司 一种数据比较方法及装置
CN110287456A (zh) * 2019-06-30 2019-09-27 张家港宏昌钢板有限公司 基于数据挖掘的大盘卷轧制表面缺陷分析方法
CN110490229A (zh) * 2019-07-16 2019-11-22 昆明理工大学 一种基于spark和聚类算法的电能表检定误差诊断方法
CN113418885A (zh) * 2021-07-22 2021-09-21 合肥学院 一种用于紫外分光光度计实验数据的分析方法
CN113673582B (zh) * 2021-07-30 2023-05-09 西南交通大学 基于系统聚类分析的铁路动态基准点多层级分群方法
CN116796214B (zh) * 2023-06-07 2024-01-30 南京北极光生物科技有限公司 一种基于差分特征的数据聚类方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004272350A (ja) * 2003-03-05 2004-09-30 Nec Corp クラスタリング装置、クラスタリング方法、クラスタリングプログラム
JP2006163894A (ja) * 2004-12-08 2006-06-22 Hitachi Software Eng Co Ltd クラスタリングシステム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2517286C2 (ru) * 2008-04-25 2014-05-27 Конинклейке Филипс Электроникс Н.В. Классификация данных выборок

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004272350A (ja) * 2003-03-05 2004-09-30 Nec Corp クラスタリング装置、クラスタリング方法、クラスタリングプログラム
JP2006163894A (ja) * 2004-12-08 2006-06-22 Hitachi Software Eng Co Ltd クラスタリングシステム

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
DAVE, R.N.: "Characterization and detection of noise in clustering", PATTERN RECOGNITION LETTERS, vol. 12, no. 11, November 1991 (1991-11-01), pages 657 - 664 *
ENDO, Y.: "Fuzzy c-Means for Data with Tolerance", PROC. 2005 INTERNATIONAL SYMPOSIUM ON NONLINEAR THEORY AND ITS APPLICATIONS (NOLTA2005), 2005, pages 345 - 348 *
HIDETOMO ICHIHASHI: "Fuzziness and Robustness in Fuzzy Clustering : Regularization, Robustification and IRLS", JOURNAL OF JAPAN SOCIETY FOR FUZZY THEORY AND INTELLIGENT INFORMATICS, vol. 17, no. 4, 15 August 2005 (2005-08-15), pages 392 - 397 *
YUKIHIRO HAMASUNA: "On Hard Clustering for Data with Tolerance", JOURNAL OF JAPAN SOCIETY FOR FUZZY THEORY AND INTELLIGENT INFORMATICS, vol. 20, no. 3, 15 June 2008 (2008-06-15), pages 388 - 398 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2016158576A (ja) * 2015-03-03 2016-09-05 富士フイルム株式会社 細胞群検出装置および方法並びにプログラム
JP2018530815A (ja) * 2015-08-17 2018-10-18 コーニンクレッカ フィリップス エヌ ヴェＫｏｎｉｎｋｌｉｊｋｅ Ｐｈｉｌｉｐｓ Ｎ．Ｖ． 生体データにおけるパターン認識のマルチレベルアーキテクチャ
JP7041614B2 (ja) 2015-08-17 2022-03-24 コーニンクレッカ フィリップス エヌ ヴェ 生体データにおけるパターン認識のマルチレベルアーキテクチャ
JP7041614B6 (ja) 2015-08-17 2022-05-31 コーニンクレッカ フィリップス エヌ ヴェ 生体データにおけるパターン認識のマルチレベルアーキテクチャ
US11327470B2 (en) 2017-12-21 2022-05-10 Mitsubishi Power, Ltd. Unit space generating device, plant diagnosing system, unit space generating method, plant diagnosing method, and program
WO2022009342A1 (ja) * 2020-07-08 2022-01-13 富士通株式会社 情報処理プログラム、情報処理方法および情報処理装置

Also Published As

Publication number Publication date
JPWO2014080447A1 (ja) 2017-01-05
JP6029683B2 (ja) 2016-11-24
US20150302042A1 (en) 2015-10-22

Similar Documents

Publication Publication Date Title
JP6029683B2 (ja) データ解析装置、データ解析プログラム
CN105096225B (zh) 辅助疾病诊疗的分析系统、装置及方法
JP6313757B2 (ja) 統合デュアルアンサンブルおよび一般化シミュレーテッドアニーリング技法を用いてバイオマーカシグネチャを生成するためのシステムおよび方法
CN112005306A (zh) 选择、管理和分析高维数据的方法和系统
US20050159896A1 (en) Apparatus and method for analyzing data
US9940383B2 (en) Method, an arrangement and a computer program product for analysing a biological or medical sample
CN109411015A (zh) 基于循环肿瘤dna的肿瘤突变负荷检测装置及存储介质
Cao et al. ROC curves for the statistical analysis of microarray data
US20200227168A1 (en) Machine learning in functional cancer assays
WO2014024142A2 (en) Population classification of genetic data set using tree based spatial data structure
CN112634987B (zh) 一种单样本肿瘤dna拷贝数变异检测的方法和装置
JP5854346B2 (ja) トランスクリプトーム解析方法、疾病判定方法、コンピュータプログラム、記憶媒体、及び解析装置
Raza et al. A novel anticlustering filtering algorithm for the prediction of genes as a drug target
CN107463797B (zh) 高通量测序的生物信息分析方法及装置、设备及存储介质
KR101067352B1 (ko) 생물학적 네트워크 분석을 이용한 마이크로어레이 실험 자료의 작용기작, 실험/처리 조건 특이적 네트워크 생성 및 실험/처리 조건 관계성 해석을 위한 알고리즘을 포함한 시스템 및 방법과 상기 방법을 수행하기 위한 프로그램을 갖는 기록매체
KR102111820B1 (ko) 동적 네트워크 바이오마커의 검출 장치, 검출 방법 및 검출 프로그램
CN101517579A (zh) 蛋白质查找方法和设备
KR20180051333A (ko) 유전체내 암 특이적 진단 마커 검출
WO2012149107A2 (en) Stratifying patient populations through characterization of disease-driving signaling
Wang et al. Learning dynamics by computational integration of single cell genomic and lineage information
JP6517933B2 (ja) 検査システム、検査装置、及び検査方法
CN110751983A (zh) 一种筛选特征mRNA用于诊断早期肺癌的方法
CN110475874A (zh) 脱靶序列在dna分析中的应用
CN116825367A (zh) 一种组织硬度预测模型的建立、组织硬度预测方法及装置
Gulla An integrated systems biology approach to investigate transcriptomic data of thyroid carcinoma

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 12888920

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2014548349

Country of ref document: JP

Kind code of ref document: A

WWE Wipo information: entry into national phase

Ref document number: 14443118

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 12888920

Country of ref document: EP

Kind code of ref document: A1