A method and system for flow cytometric classification and enumeration based on density estimation

Abstract

The application provides a method and system for flow cytometry classification and counting based on density estimation, and the method comprises the following steps: step 1, reading an original flow data file to obtain a data set; step 2, standardizing the data set; step 3, automatically reducing the dimension of the data set by using a UMAP algorithm; step 4, analyzing the cell distribution based on a density method, and automatically clustering a cell population with uniform distribution by using a DBSCAN algorithm; and step 5, performing flow cytometry counting. The method adopted by the application has a faster processing speed and higher calculation accuracy, and has important guiding significance in cell detection.

Description

Technical Field

[0001] This invention belongs to the technical field of medical data processing and flow cytometry data analysis, and particularly relates to a method and system for flow cytometry cell classification and counting based on density estimation. Background Technology

[0002] In routine blood tests, the number of red blood cells and white blood cells plays a crucial role in diagnosing many diseases. An abnormal increase or decrease in red blood cell or white blood cell counts suggests a potential problem with the patient's bodily functions. Doctors prescribe medication based on changes in cell counts, combined with the patient's symptoms and other relevant test results, to make an accurate diagnosis.

[0003] Flow cytometry (FCM), developed in the 1970s, is a technique for precise and rapid multi-parameter quantitative analysis of the physicochemical and biological characteristics of biological cells, as well as the sorting of specific cell populations. A flow cytometer is a typical instrument for flow cytometry. It borrows from fluorescence microscopy and utilizes advancements in fluorescent dyes, laser technology, monoclonal antibody technology, and computer technology. It can rapidly analyze the content and types of chemical substances in individual cells and their populations, and isolate and purify cells containing specified chemical substances. It is widely used in some laboratories and hospitals. Its principle utilizes the focusing principle of fluid dynamics. The cells or particles to be analyzed are arranged in a line and rapidly flow through a detection beam one by one. The multi-angle scattered light and multi-color fluorescence emitted by each cell or particle are collected by an optical system and converted into electrical signals by a photoelectric sensor. After electronic signal processing and sampling, these signals become digital signals, which are then stored and analyzed by a computer. The characteristic data of all cells or particles acquired by the flow cytometer become flow cytometry data. Flow cytometry (FCM) has significantly improved detection speed and statistical accuracy, and can simultaneously measure multiple parameters from a single cell. With its advantages of high speed, high precision, good accuracy, large-scale analysis, and multi-parameter analysis, FCM has become the most advanced quantitative cell analysis technique and an important clinical testing instrument. In recent years, FCM has not only made significant progress in biomedicine and clinical laboratory science, but has also expanded into various fields of biology, including microbiology, medicine, fermentation, and environmental protection. Furthermore, it remains an indispensable tool in cell research.

[0004] Traditionally, the analysis of flow cytometry data relies on experienced personnel projecting the data into a two-dimensional scatter plot, and then using region gating to analyze groups of interest, such as classification and counting; this is known as manual gating. Flow cytometers, as well as flow cytometry-based blood analyzers, urine analyzers, and particle analyzers, all collect and analyze two-dimensional or multi-dimensional particle data. The collected signals generate two-dimensional or three-dimensional scatter plots, which are then divided into multiple regions. Particles whose multi-parameter signals of cells or particles fall within the same region are grouped into the same category, and the number and percentage of particles falling into these categories are counted to analyze the statistical characteristics of the sample, thereby identifying different particles in the liquid and classifying them into different categories.

[0005] As flow cytometry develops towards multi-laser and high-throughput technologies, the amount of data generated by flow cytometers has increased exponentially. Traditional manual data analysis methods can no longer meet the needs of rapid analysis, and rapid automatic data analysis has become the main direction for the future development of flow cytometry.

[0006] Data clustering analysis has become a very active research topic, but few people have studied the use of clustering methods to analyze cell classification data detected by flow cytometry.

[0007] Clustering is an unsupervised classification process that divides a dataset into groups or clusters, ensuring that data objects within the same group are highly similar, while those in different groups are less similar. The key to clustering is grouping similar data together, making it suitable for exploring the internal relationships between samples. When clustering large datasets, the K-means algorithm is also much faster than hierarchical clustering algorithms.

[0008] The steps for cell classification in flow cytometry using K-means involve detecting fluorescently labeled cells using FCM to obtain fluorescence spectral data, preprocessing this data, and then constructing a clustering model using the K-means algorithm to achieve cell classification. The key to studying the application of the K-means algorithm in cell classification in flow cytometry is to use this algorithm to extract effective information from the fluorescence spectra and establish a correct mapping relationship between fluorescence spectra and cell classification. K-means clustering accuracy is poor; it is only suitable for accurately clustering and distinguishing cell populations that exhibit a near-circular distribution after dimensionality reduction. However, it has extremely poor clustering ability for irregularly distributed cell populations after dimensionality reduction. Furthermore, in the K-means algorithm, k is given in advance, and selecting this value is very difficult. Often, the optimal number of categories for a given dataset is unknown beforehand. The choice of initial cluster centers has a significant impact on the clustering results. If the initial values ​​are poorly chosen, effective clustering results may not be obtained, which is an obstacle to automatic clustering, introducing human error and hindering the use of automated engineering and the discovery and mining of novel cell populations. Therefore, there is an urgent need for a rapid and efficient flow cytometry method and system for cell classification and counting based on density estimation. Summary of the Invention

[0009] Purpose of the invention: The technical problem to be solved by the present invention is to provide a method and system for flow cytometry cell classification and counting based on density estimation, which addresses the shortcomings of the prior art.

[0010] The method includes the following steps: Step 1, read the original streaming data file to obtain the dataset;

[0011] Step 2: Standardize the dataset;

[0012] Step 3: Use the UMAP (Uniform Manifold Approximation and Projection for Dimension Reduction) algorithm to automatically reduce the dimensionality of the dataset;

[0013] Step 4: Analyze cell distribution based on density method. For uniformly distributed cell populations, use the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) algorithm for automatic clustering.

[0014] Step 5: Perform flow cytometry cell counting.

[0015] Step 2 includes: performing min-max normalization on the dataset, transforming the original streaming data in the dataset using the following formula:

[0016]

[0017] Where i is a variable, X[i] represents the i-th data in the original streaming data, X′[i] is the data after X[i] is standardized, i = 1, 2, ..., N, N represents the number of data in the original streaming data, min{X[i]} represents the minimum value among the N data, and max{X[i]} represents the maximum value among the N data.

[0018] z-score standardization, also known as standard deviation standardization, is used to assess the distance of a sample point from the population mean. It is suitable when the maximum and minimum values ​​of attribute A are unknown, or when there are outliers outside the range of values. Specifically:

[0019]

[0020] in, Represents the sample mean. S represents the standard deviation of the sample. X[i] represents the original data in the i-th dimension of the original streaming data, and X′[i] represents the standardized data.

[0021] Scaling the data using z-score normalization or min-max normalization makes it easier to find the optimal cluster radius in step 4.2. In this invention, min-max normalization is used to cluster the data within the range of 0-1.

[0022] Step 4 includes:

[0023] Step 4-1: Use SPSS (Statistical Product and Service Solutions) chi-square test to determine whether the data in the dataset follows a uniform distribution.

[0024] Step 4-2: When the data distribution is uniform, the DBSCAN algorithm is used to automatically cluster the cell distribution.

[0025] Step 4-3: When the data distribution is uneven, mean offset clustering is used to automatically cluster the cell distribution.

[0026] Step 4-2 includes:

[0027] Step 4-2-1, determine the radius r1 and the minimum number of points minPoints: For the data points in the dataset after dimensionality reduction in Step 3, starting from an arbitrarily unvisited data point, determine whether the number of points contained within a circle centered on the data point with a radius of r1 is greater than or equal to the minimum number of points minPoints within the region. If it is greater than or equal to the minimum number of points minPoints within the region, the data point is marked as a core point; otherwise, proceed to Step 4-2-2;

[0028] Among them, the value of r1 is crucial. If it is set very small, it means that no point is a core sample, which may cause all points to be marked as noise. If r1 is set very large, it may cause all points to form a single cluster. Since the data has been normalized to numbers between 0 and 1 in Step 2, the setting of r1 is generally less than 0.5. In the present invention, r1 = 0.3;

[0029] The value range of the minimum number of points minPoints is 3 < minPoints < 10. Such a setting can filter out outliers while retaining the points that are abnormal in the cell population. When minPoints is set too large, the abnormal points will be judged as noise points. When minPoints is set too small, the noise points will also be detected as core points.

[0030] Step 4-2-2, if a data point exists within a circle centered on a core point with a radius of r, the data point is marked as an edge point; otherwise, it is a noise point. Repeat Step 4-2-1 to Step 4-2-2 until all data points have been visited.

[0031] Step 4-3 includes:

[0032] Step 4-3-1, determine the sliding window radius r2, randomly select a data point in the dataset after dimensionality reduction in Step 3 as the center of the circle, and start sliding a circular sliding window with a radius of r2;

[0033] Among them, the value of r2 is crucial. If it is set very small, it means that no point is a core sample, which may cause all points to be marked as noise. If r2 is set very large, it may cause all points to form a single cluster. Since the data has been normalized to numbers between 0 and 1 in Step 2, the setting of r2 is generally less than 0.5. In the present invention, r2 = 0.3.

[0034] Step 4-3-2, every time it slides to a new region, calculate the mean value within the sliding window as the center point, and the number of points within the sliding window is the density within the window; in each movement, the window will move towards the region with higher density;

[0035] Step 4-3-3: Move the window while calculating the center point and density within the window, until the density inside the circle no longer increases.

[0036] Steps 4-3-4 and 4-3-1 to 4-3-3 will generate more than two sliding windows. When more than two sliding windows overlap, the window containing the most points is retained. Then, clustering is performed according to the sliding window where the data points are located. After clustering, different colors are used to label cells of different categories to obtain a clustered color image. There are a total of N cell categories after clustering.

[0037] Step 5 includes:

[0038] Step 5-1: Convert the clustered color image into a cell binarized image:

[0039] For each type of cell, count the cells. First, for the first type of cells after clustering, binarize the first type of cells into white in the clustered color image, and binarize the remaining points other than the first type of cells into black to obtain the cell binarized image.

[0040] Step 5-2: Perform distance transformation on the binarized cell image to obtain the distance-transformed binary cell image;

[0041] Step 5-3: Use the findContours() contour extraction algorithm in OpenCV to extract all contours of the distance-transformed binary cell image and obtain the number of contours.

[0042] Step 5-4: Use the contourArea algorithm in OpenCV to calculate the area of ​​all contours;

[0043] Step 5-5: If the area of ​​a contour is greater than the threshold th1, proceed to step 5-6; otherwise, directly return the number of contours as the counting result for each type of cell.

[0044] Steps 5-6: Cell counting based on density images;

[0045] Steps 5-7: For cells in other clusters after clustering, repeat steps 5-1 to 5-6 until the second, third, ..., jth, ... Nth clusters are counted.

[0046] Where j is a variable, representing which cell category is being counted, and N is the total number of cell clusters obtained in step 4.

[0047] In step 5-1, the cell binarization image is obtained using the following formula:

[0048]

[0049] Where j is a variable, j = 1, 2, ..., N, B j [m,n]=255 means that the pixel value at the pixel with x-coordinate m and y-coordinate n in the binarized image of the j-th cell is 255. m is a variable representing the m-th column in the binarized image, and n is a variable representing the n-th row in the binarized image. m=1,2,...,W, n=1,2,...,H, W is the width of the binarized image, and H is the height of the binarized image.

[0050] Step 5-2 includes: The distance transformation is to calculate the distance from a non-zero pixel in an image to the nearest zero pixel, as shown in the following formula:

[0051]

[0052] Among them, B j [m,n] represents the non-zero pixels in the binary image of the j-th cell type, B j [p,q] represents the zero pixel in the binarized image of the j-th cell class, where p is a variable representing the p-th column in the binarized image, and q is a variable representing the q-th row in the binarized image. j [p, q] = 0 indicates that the pixel value at the x-coordinate p and y-coordinate q in the binary image of the j-th cell type is 0;

[0053] D(B j [m, n] = 255, B j [p, q] = 0) represents the non-zero pixel B in the binary image of the j-th cell type. j [m, n] and zero pixel B j The distance between [p, q], P j [m, n] represents the pixel value at the pixel with x-coordinate m and y-coordinate n in the binary image after distance transformation; min(·) represents the minimum value function, and D(P, Q) represents the Euclidean distance between points P and Q.

[0054] th2 is the distance threshold between non-zero pixels and zero pixels in the binary image of cell j. The general value of th2 is less than 10. The setting of th is mainly to separate adhered or close cells. If it is set too high, too many cells will be filtered out, making the counting results inaccurate.

[0055] In step 5-5, th1 is the area threshold for determining severe cell-cell adhesion, where th1 < 50. In this invention, the specific value of th1 is determined based on the specific cell type.

[0056] Steps 5-6 include:

[0057] Step 5-6-1: Binarize the image P of the j-th cell type. j [m, n] are convolved with the Gaussian filtered template image G[u, v] to obtain the cell density image D of the j-th type of cells. j [m, n]:

[0058]

[0059]

[0060] Where u represents the x-coordinate of the Gaussian filter template image, v represents the y-coordinate of the Gaussian filter template image, σ is the standard deviation of the Gaussian image, and e represents the base of the exponential function, e = 2.71828183, which is the natural constant. This is a convolution operation;

[0061] In this invention, when the scale of the Gaussian filter template image is set to 5, u = -2, -1, ..., 2, v = -2, -1, ..., 2, σ < 2.

[0062] Gaussian filtering works by performing convolution calculations on each point of the input array with the input Gaussian filter template, and then combining these results into a filtered output array. In simple terms, Gaussian filtering is a weighted average process for the entire image, where the value of each pixel is obtained by weighting its own value and the values ​​of other pixels in its neighborhood.

[0063] Step 5-6-2: Convert the density image of cell type j to D j The summation of the values ​​at all pixels in the range [m,n] gives the number of cells in the j-th class. This process can be repeated to obtain the number of cells in the 1st, 2nd, ..., Nth classes.

[0064] The present invention also provides a system for flow cytometry cell classification and counting based on density estimation, including a data acquisition unit, a preprocessing unit, a dimensionality reduction unit, a judgment unit, a clustering unit, and a counting unit;

[0065] The data acquisition unit is used to acquire flow cytometry data: when cells enter the flow cytometer, after passing through multicolor fluorescence, they are collected by the optical system and converted into electrical signals by photoelectric sensors. The electrical signals are converted into digital signals by A / D (analog-to-digital converter). The digital signals are called feature data. The data acquisition unit stores all the feature data of cells or particles acquired by the flow cytometer, which is called the raw flow cytometry data file. The raw flow cytometry data file is stored in the dataset.

[0066] The preprocessing unit is used to standardize the dataset, and then send the processed dataset into the dimensionality reduction unit.

[0067] The dimensionality reduction unit uses the UMAP algorithm to automatically reduce the dimensionality of the dataset;

[0068] The judgment unit uses the SPSS chi-square test to determine whether the data in the dataset follows a uniform distribution. When the data distribution is uniform, the clustering unit uses the DBSCAN algorithm to automatically cluster the cell distribution; otherwise, the clustering unit uses mean-shift clustering to automatically cluster the cell distribution.

[0069] The counting unit is used to count the cells in each clustered cell image separately.

[0070] The beneficial effects of this invention are:

[0071] 1. It can automatically cluster and count cell populations of different densities, and can effectively eliminate noise interference and non-specific signals.

[0072] 2. The method used in this invention has a faster processing speed and higher calculation accuracy, and has important guiding significance in cell detection. Attached Figure Description

[0073] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments, and the advantages of the present invention in the above and / or other aspects will become clearer.

[0074] Figure 1 This is a flowchart of the method of the present invention.

[0075] Figure 2 This is a counting flowchart.

[0076] Figure 3 This is a schematic diagram of a 3×3 two-dimensional Gaussian template.

[0077] Figure 4 This is a schematic diagram showing the binary image before and after convolution with a template image. Detailed Implementation

[0078] This invention provides a flow cytometry cell classification and counting system based on density estimation, comprising a data acquisition unit, a preprocessing unit, a dimensionality reduction unit, a judgment unit, a clustering unit, and a counting unit. Details are as follows:

[0079] (1) Data acquisition unit

[0080] After cells enter the flow cytometer, they undergo multicolor fluorescence, are collected by an optical system, and converted into electrical signals by a photoelectric sensor. After electronic signal processing and sampling, they become digital signals, which are then stored and analyzed by a computer. All characteristic data of cells or particles acquired by the flow cytometer are called flow cytometry data.

[0081] (2) Preprocessing unit

[0082] The streaming data is standardized and then fed into the dimensionality reduction unit.

[0083] (3) Dimensionality reduction unit

[0084] If the feature dimensionality of the dataset is too high, it is difficult to construct accurate clusters. High dimensionality does not necessarily mean hundreds or thousands of features; even 10-dimensional features can cause accuracy problems. After standardization, streaming data undergoes dimensionality reduction to transform high-dimensional data into low-dimensional data.

[0085] (4) Judgment Unit

[0086] Use the SPSS chi-square test to test whether the data follows a uniform distribution.

[0087] (5) Clustering unit

[0088] When the distribution follows a uniform distribution, the DBSCAN algorithm is used to automatically cluster the cell distribution; otherwise, mean offset clustering is used to automatically cluster the cell distribution.

[0089] (6) Counting unit

[0090] Count the cells in each cluster image separately.

[0091] Therefore, the present invention also provides a method for flow cytometry cell classification and counting based on density estimation, the flowchart of which is shown below. Figure 1 , Figure 2 As shown, the specific steps include the following:

[0092] Step 1: Read the raw streaming data file to obtain the dataset;

[0093] Step 2: Standardize the dataset;

[0094] Because the features in the dataset are not within the same range, the entire dataset needs to be standardized. In other words, each feature in the dataset of this invention has a unique size and range for its data.

[0095] Step 3: Use the UMAP algorithm to automatically reduce the dimensionality of the dataset;

[0096] Step 4: Analyze cell distribution based on density methods. For uniformly distributed cell populations, the DBSCAN algorithm is used for automatic clustering; when clusters have significantly different densities, mean-shift clustering is used. Specifically, this includes:

[0097] SPSS chi-square test is used to test whether the data in the dataset follows a uniform distribution. If the data distribution is uniform, the DBSCAN algorithm is used to automatically cluster the cell distribution; otherwise, mean shift clustering is used to automatically cluster the cell distribution.

[0098] Step 5: Perform flow cytometry cell counting. This specifically includes:

[0099] Step 5-1: For the first type of cells after clustering, binarize this type of cell into white in the clustered color image, and binarize the remaining points into black to obtain the cell binarized image. The calculation formula is:

[0100]

[0101] Where j is a variable, j = 1, 2, ..., N, B j [m, n] = 255 means that the pixel value at the pixel with x-coordinate m and y-coordinate n in the binarized image of the j-th cell is 255. m is a variable representing the m-th column in the binarized image, and n is a variable representing the n-th row in the binarized image. m = 1, 2, ..., W, n = 1, 2, ..., H, W is the width of the binarized image, and H is the height of the binarized image.

[0102] Step 5-2: Perform distance transformation on the binarized cell image;

[0103]

[0104] Among them, B j [m, n] represents the non-zero pixels in the binary image of the j-th cell type, B j [p, q] represents the zero pixel in the binarized image of the j-th cell, where p is a variable representing the p-th column in the binarized image and q is a variable representing the q-th row in the binarized image. That is, the value is 0 at the pixel with horizontal coordinate p and vertical coordinate q in the binarized image.

[0105] D(B j [m, n] = 255, B j [p, q] = 0) represents the non-zero pixel B in the binary image of the j-th cell type. j [m, n] and zero pixel B j The distance between [p, q], P j [m, n] represents the pixel value at the pixel with x-coordinate m and y-coordinate n in the binary image after distance transformation; min(·) represents the minimum value function, and D(P, Q) represents the Euclidean distance between points P and Q.

[0106] th2 is the distance threshold between non-zero pixels and zero pixels in the binary image of cell j. The general value of th2 is less than 10. The setting of th is mainly to separate adhered or close cells. If it is set too large, too many cells will be filtered out, making the counting results inaccurate. In the embodiment of the present invention, th2 = 3.

[0107] Step 5-3: Use the findContours() contour extraction algorithm in OpenCV to extract all contours of the image and obtain the number of contours;

[0108] Step 5-4: Use the contourArea algorithm in OpenCV to calculate the area of ​​all contours;

[0109] Step 5-5: If the area of ​​the contour is greater than th1, proceed to step 5-6; otherwise, directly return the number of contours as the counting result of this type of cell.

[0110] th1 is the area threshold for determining severe cell-cell adhesion, where th1 < 50. In this invention, the specific value of th1 is determined based on the specific cell type.

[0111] Steps 5-6: Cell counting based on density images;

[0112] (1) Obtain the density image D of the j-th cell type by performing Gaussian convolution on the binary image of the j-th cell type. j [m, n];

[0113] Gaussian filtering works by convolving each point of the input array with the input Gaussian filter template, and then combining these results into a filtered output array. In simpler terms, Gaussian filtering is a weighted average process for the entire image; the value of each pixel is obtained by weighting its own value and the values ​​of its neighboring pixels. The binarized image B of the j-th cell type... j [m, n], the cell density image of the j-th class is obtained by convolving it with the Gaussian filtered template image. The calculation formula is:

[0114]

[0115] in, This is a convolution operation;

[0116] Advanced convolutional images are represented by the following formula:

[0117]

[0118] Where G[u, v] represents the Gaussian filter template image, u represents the x-coordinate u in the Gaussian filter template image, v represents the y-coordinate v in the Gaussian filter template image, σ is the standard deviation of the Gaussian filter template image, in this invention σ=1.5, e represents the base of the exponential function, e=2.71828183, which is the natural constant.

[0119] The specific operation of Gaussian filtering is as follows: scan each pixel in the image with a template (or convolution, mask), and replace the value of the center pixel of the template with the weighted average gray value of the pixels in the neighborhood determined by the template.

[0120] Taking a 3×3 Gaussian template as an example (to prevent overflow, the image elements in the leftmost column, rightmost column, topmost row, and bottommost row are not processed), the basic steps of Gaussian filtering are as follows:

[0121] 1) Scan the image from top to bottom and from left to right;

[0122] 2) Move the Gaussian template to a non-zero point in the image, perform a convolution operation between the coefficients of the Gaussian template and the area it covers, and replace the value at that point with the weighted average value. Repeat this process until all points in the image have been traversed.

[0123] A 3×3 two-dimensional Gaussian template, such as Figure 3 As shown, Figure 4 The image shown is a schematic diagram of the density image after a binary image is convolved with a template image.

[0124] (2) Sum the probabilities at all locations in the density image to obtain the number of cells of each class.

[0125] The number of cells in each cell type is obtained by summing up the values ​​of pixels at all locations in the density image.

[0126] Example 1

[0127] This embodiment provides a flow cytometry cell classification and counting system based on density estimation, including a data acquisition unit, a preprocessing unit, a dimensionality reduction unit, a judgment unit, a clustering unit, and a counting unit. Details are as follows:

[0128] (1) Data acquisition unit

[0129] After cells enter the flow cytometer, they undergo multicolor fluorescence, are collected by an optical system, and converted into electrical signals by a photoelectric sensor. After electronic signal processing and sampling, they become digital signals, which are then stored and analyzed by a computer. All characteristic data of cells or particles acquired by the flow cytometer are called flow cytometry data.

[0130] (2) Preprocessing unit

[0131] The streaming data is standardized and then fed into the dimensionality reduction unit.

[0132] (3) Dimensionality reduction unit

[0133] If the feature dimensionality of the dataset is too high, it is difficult to construct accurate clusters. High dimensionality does not necessarily mean hundreds or thousands of features; even 10-dimensional features can cause accuracy problems. After standardization, streaming data undergoes dimensionality reduction to transform high-dimensional data into low-dimensional data.

[0134] (4) Judgment Unit

[0135] Use the SPSS chi-square test to test whether the data follows a uniform distribution.

[0136] (5) Clustering unit

[0137] When the distribution follows a uniform distribution, the DBSCAN algorithm is used to automatically cluster the cell distribution; otherwise, mean offset clustering is used to automatically cluster the cell distribution.

[0138] (6) Counting unit

[0139] Count the cells in each cluster image separately.

[0140] Therefore, the present invention also relates to a flow cytometry method for cell classification and counting based on density estimation, the flowchart of which is shown below. Figure 1 Specifically, it is shown as follows:

[0141] Step 1: Read the raw streaming data file to obtain the dataset;

[0142] Step 2: Standardize the dataset;

[0143] Because the features in the dataset are not in the same range, the entire dataset needs to be standardized. In other words, each feature in the dataset has a unique size and range for its data.

[0144] Step 3: Use the UMAP algorithm to automatically reduce the dimensionality of the dataset;

[0145] In some algorithms, it becomes difficult to accurately construct clusters if the feature dimensionality of the dataset is too high. High dimensionality does not necessarily mean hundreds or thousands of features; even 10-dimensional features can cause accuracy problems.

[0146] The theory behind feature reduction is to transform the original feature set into fewer artificially derived features that still retain most of the information contained in the original features.

[0147] In this invention, the UMAP algorithm is used to automatically reduce the dimensionality of the data.

[0148] Step 4: Analyze cell distribution based on density method. For uniformly distributed cell populations, use the DBSCAN algorithm for automatic clustering; specifically including:

[0149] Step 4-1: Use the SPSS chi-square test to test whether the data in the dataset follows a uniform distribution.

[0150] Step 4-2: When the data distribution is uniform, the DBSCAN algorithm is used to automatically cluster the cell distribution.

[0151] DBSCAN is a density-based clustering algorithm. This type of density clustering algorithm generally assumes that clusters can be determined by the density of sample distribution. Samples of the same cluster are closely connected; that is, any sample of that cluster will have other samples of the same cluster within a short distance of it.

[0152] By grouping closely connected samples into one class, we obtain a cluster category. By further dividing all closely connected samples into distinct categories, we obtain the final cluster category results.

[0153] The DBSCAN algorithm defines a cluster simply as the set of samples connected by the maximum density derived from the density reachability relation, which is the final cluster.

[0154] The DBSCAN algorithm allows one or more core points within a cluster. If there is only one core point, all other non-core point samples in the cluster are within the Eps neighborhood of this core point. If there are multiple core points, then any core point in the Eps neighborhood of any given core point must contain another core point; otherwise, these two core points cannot be density-reached. The set of all samples within the Eps neighborhood of these core points forms a DBSCAN cluster.

[0155] The DBSCAN processing flow is as follows.

[0156] Step 4-2-1, Determine the radius r1 and the minimum number of points minPoints: For the data points in the dimensionality-reduced dataset in Step 3, starting from any unvisited data point, determine whether the number of points contained within a circle with the data point as the center and r1 as the radius is greater than or equal to the minimum number of points minPoints in the region. If it is greater than or equal to the minimum number of points minPoints in the region, the data point is marked as a core point; otherwise, proceed to Step 4-2-2.

[0157] Among them, the value of r1 is crucial. If it is set very small, it means that no point is a core sample, which may cause all points to be marked as noise. If r1 is set very large, it may cause all points to form a single cluster. Since the data has been normalized to numbers between 0 and 1 in step 2, the setting of r1 is generally less than 0.5. In the present invention, r1 = 0.3;

[0158] The value range of the minimum number of points minPoints is 3 < minPoints < 10. Such a setting can filter outlier points while retaining the points with abnormalities in the cell population. When minPoints is set too large, the points with abnormalities will be judged as noise points. When minPoints is set too small, the noise points will also be detected as core points.

[0159] Step 4-2-2, if a data point exists within a circle centered on a core point with a radius of r, the data point is marked as an edge point, otherwise it is a noise point. Repeat steps 4-2-1 to step 4-2-2 until all data points have been visited.

[0160] Step 4-3, when the data distribution is uneven, mean shift clustering is used to automatically cluster the cell distribution;

[0161] Mean shift clustering is an algorithm based on a sliding window to find the dense regions of data points. This is a centroid-based algorithm, which is completed by updating the candidate points of the center point to the mean of the points within the sliding window to locate the center points of each group / class. Then, similar windows are removed from these candidate windows, and finally a set of center points and corresponding groupings are formed. The specific steps include:

[0162] Step 4-3-1, determine the sliding window radius r2, randomly select a data point from the dataset after dimensionality reduction in step 3 as the center of the circle, and start sliding a circular sliding window with a radius of r2;

[0163] Among them, the value of r2 is crucial. If it is set very small, it means that no point is a core sample, which may cause all points to be marked as noise. If r2 is set very large, it may cause all points to form a single cluster. Since the data has been normalized to numbers between 0 and 1 in step 2, the setting of r2 is generally less than 0.5. In the present invention, r2 = 0.3.

[0164] Step 4-3-2, every time it slides to a new area, calculate the mean within the sliding window as the center point, and the number of points within the sliding window is the density within the window; in each movement, the window will move towards the area with higher density;

[0165] Step 4-3-3: Move the window while calculating the center point and density within the window, until the density inside the circle no longer increases.

[0166] Steps 4-3-4 and 4-3-1 to 4-3-3 will generate more than two sliding windows. When more than two sliding windows overlap, the window containing the most points is retained. Then, clustering is performed according to the sliding window where the data points are located. After clustering, different colors are used to label cells of different categories to obtain a clustered color image. There are a total of N cell categories after clustering.

[0167] Step 5: Perform flow cytometry cell counting.

[0168] Step 5-1: Convert the clustered color image into a cell binarized image: For the j-th cell after clustering, binarize this cell type as white in the clustered color image, and binarize the remaining points as black, to obtain the binarized image of the j-th cell type. The calculation formula is as follows:

[0169]

[0170] Where j is a variable, j = 1, 2, ..., N, B j [m, n] = 255 means that the pixel value at the pixel with x-coordinate m and y-coordinate n in the binarized image of the j-th cell is 255. m is a variable representing the m-th column in the binarized image, and n is a variable representing the n-th row in the binarized image. m = 1, 2, ..., W, n = 1, 2, ..., H, W is the width of the binarized image, and H is the height of the binarized image.

[0171] Step 5-2: Perform distance transformation on the binarized cell image;

[0172]

[0173] Among them, B j [m, n] represents the non-zero pixels in the binary image of the j-th cell type, B j [p, q] represents the zero pixel in the binarized image of the j-th cell, where p is a variable representing the p-th column in the binarized image and q is a variable representing the q-th row in the binarized image. That is, the value is 0 at the pixel with horizontal coordinate p and vertical coordinate q in the binarized image.

[0174] D(B j [m, n] = 255, B j [p, q] = 0) represents the non-zero pixel B in the binary image of the j-th cell type. j [m, n] and zero pixel B j The distance between [p, q], P j[m, n] represents the pixel value at the pixel with x-coordinate m and y-coordinate n in the binary image after distance transformation; min(·) represents the minimum value function, and D(P, Q) represents the Euclidean distance between points P and Q.

[0175] th2 is the distance threshold between non-zero pixels and zero pixels in the binarized image of the j-th type of cell. The general value of th2 is 10, and in this invention, th2 = 3.

[0176] Step 5-3: Use the findContours() contour extraction algorithm in OpenCV to extract all contours of the image and obtain the number of contours;

[0177] Step 5-4: Use the contourArea algorithm in OpenCV to calculate the area of ​​all contours;

[0178] Step 5-5: If the area of ​​the contour is greater than th1, proceed to step 5-6; otherwise, directly return the number of contours as the counting result of this type of cell.

[0179] th1 is the area threshold for determining severe cell-cell adhesion, where th1 < 50. In this invention, the specific value of th1 is determined based on the specific cell type.

[0180] Steps 5-6: Cell counting based on density images;

[0181] Step 5-6-1: Binarize the image P of the j-th cell type. j [m, n] are convolved with the Gaussian filtered template image G[u, v] to obtain the cell density image D of the j-th type of cells. j [m, n]:

[0182] Gaussian filtering works by convolving each point of the input array with the input Gaussian filter template, and then combining these results into a filtered output array. In simpler terms, Gaussian filtering is a weighted average process for the entire image; the value of each pixel is obtained by weighting its own value and the values ​​of its neighboring pixels. The binarized image B of the j-th cell type... j [m, n], the cell density image of the j-th class is obtained by convolving it with the Gaussian filtered template image. The calculation formula is:

[0183]

[0184]

[0185] Where u represents the horizontal coordinate u in the Gaussian filter template image, v represents the vertical coordinate v in the Gaussian template image, σ is the standard deviation of the Gaussian image, in this invention σ=1.5, e represents the base of the exponential function, e=2.71828183, which is a natural constant.

[0186] Taking a 3×3 Gaussian template with σ=1.5 as an example, the corresponding Gaussian filter template images are shown in Table 1:

[0187] Table 1

[0188] 0.0453542 0.0566406 0.0453542 0.0566406 0.0707355 0.0566406 0.043542 0.0566406 0.0453542

[0189] To ensure that the sum of the weights of the nine values ​​in Table 1 equals 1, they must be normalized, resulting in Table 2:

[0190] Table 2

[0191] 0.0947416 0.118318 0.0947416 0.118318 0.147761 0.118318 0.0947416 0.118318 0.0947416

[0192] The specific operation of Gaussian filtering is as follows: scan each pixel in the image with a template (or convolution, mask), and replace the value of the center pixel of the template with the weighted average gray value of the pixels in the neighborhood determined by the template.

[0193] Taking a 3×3 Gaussian template as an example (to prevent overflow, the image elements in the leftmost column, rightmost column, topmost row, and bottommost row are not processed), the basic steps of Gaussian filtering are as follows:

[0194] 1) Scan the image from top to bottom and from left to right;

[0195] 2) Move the Gaussian filter template image to a non-zero point in the image, perform a convolution operation between the coefficients of the Gaussian filter template and the area it covers, and replace the value at that point with the weighted average value. Repeat this process until all points in the image have been traversed.

[0196] Step 5-6-2: Convert the density image of cell type j to D j The summation of the values ​​at all pixels in the range [m,n] gives the number of cells in the j-th class. This process can be repeated to obtain the number of cells in the 1st, 2nd, ..., Nth classes.

[0197] Example 2:

[0198] In this embodiment, as Figure 3 As shown, a 3×3 two-dimensional Gaussian filter template image is used as an example for illustration. Figure 4 The image shown is a schematic diagram of the density image after a binary image is convolved with a Gaussian filter template image.

[0199] The number of cells in each cell type is obtained by summing the probabilities at all locations in the density image: specifically, by accumulating the values ​​of pixels at all locations in the density image.

[0200] like Figure 4 As shown, Figure 4 The leftmost image is the binary image, the middle image is a 3×3 two-dimensional Gaussian filtered template image, and the rightmost image is the density image. The final cell count is:

[0201] This invention uses the OpenCV open-source FindContours contour extraction function to find contours in a binary image. It works by systematically scanning the image until it encounters a point in a connected region, using that point as the starting point to trace its contour, marking the elements on the boundary. When the contour is completely closed, the scan returns to the previous position until a new component is found.

[0202] In its specific implementation, this application provides a computer storage medium and a corresponding data processing unit. The computer storage medium is capable of storing a computer program, which, when executed by the data processing unit, can run the invention's content regarding a flow cytometry method for classification and counting based on density estimation, as well as some or all of the steps in various embodiments. The storage medium can be a magnetic disk, optical disk, read-only memory (ROM), or random access memory (RAM), etc.

[0203] Those skilled in the art will clearly understand that the technical solutions in the embodiments of the present invention can be implemented using computer programs and their corresponding general-purpose hardware platforms. Based on this understanding, the technical solutions in the embodiments of the present invention, or the parts that contribute to the prior art, can be embodied in the form of computer programs, i.e., software products. These computer program software products can be stored in a storage medium and include several instructions to cause a device containing a data processing unit (which may be a personal computer, server, microcontroller, MUU, or network device, etc.) to execute the methods described in various embodiments or certain parts of the embodiments of the present invention.

[0204] This invention provides a method and system for flow cytometry cell classification and counting based on density estimation. Many methods and approaches exist for implementing this technical solution; the above description is merely a preferred embodiment of the invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications should also be considered within the scope of protection of this invention. All components not explicitly stated in this embodiment can be implemented using existing technologies.

Claims

1. A method for flow cytometry cell classification and counting based on density estimation, characterized in that, Includes the following steps: Step 1: Read the raw streaming data file to obtain the dataset; Step 2: Standardize the dataset; Step 3: Use the UMAP algorithm to automatically reduce the dimensionality of the dataset; Step 4: Analyze cell distribution based on density method, and use DBSCAN algorithm to automatically cluster uniformly distributed cell populations; Step 5: Perform flow cytometry cell counting; Step 4 includes: Step 4-1: Use SPSS chi-square test to determine whether the data in the dataset follows a uniform distribution; Step 4-2: When the data distribution is uniform, the DBSCAN algorithm is used to automatically cluster the cell distribution. Step 4-3: When the data distribution is uneven, mean shift clustering is used to automatically cluster the cell distribution. Step 4-2 includes: Step 4-2-1, Determine the radius r1 and the minimum number of points minPoints: For the data points in the dimensionality-reduced dataset in Step 3, starting from any unvisited data point, determine whether the number of points contained within a circle with the data point as the center and r1 as the radius is greater than or equal to the minimum number of points minPoints in the region. If it is greater than or equal to the minimum number of points minPoints in the region, the data point is marked as a core point; otherwise, proceed to Step 4-2-2. Step 4-2-2: If a data point exists within a circle centered on the core point with radius r, the data point is marked as an edge point; otherwise, it is marked as a noise point. Step 4-2-3: Repeat steps 4-2-1 to 4-2-2 until all data points have been accessed; Step 4-3 includes: Step 4-3-1: Determine the sliding window radius r2, randomly select a data point in the dimensionality-reduced dataset from Step 3 as the center of the circle, and start sliding with a circular sliding window of radius r2; Step 4-3-2: Each time the window slides to a new region, the mean value within the sliding window is calculated and used as the center point. The number of points within the sliding window is the density within the window. In each move, the window will move towards a region with higher density. Step 4-3-3: Move the window while calculating the center point and density within the window, until the density inside the circle no longer increases. Steps 4-3-4 and 4-3-1 to 4-3-3 will generate more than two sliding windows. When more than two sliding windows overlap, the window containing the most points is retained. Then, clustering is performed according to the sliding window where the data points are located. After clustering, different colors are used to label cells of different categories to obtain a clustered color image. There are a total of N cell categories after clustering.

2. The method according to claim 1, characterized in that, Step 2 includes: performing min-max normalization on the dataset, transforming the original streaming data in the dataset using the following formula: Where i is a variable, X[i] represents the i-th data in the original streaming data, X′[i] is the data after standardization of X[i], i = 1, 2, ..., N, N represents the number of data in the original streaming data, min{X[i]} represents the minimum value among the N data, and max{X[i]} represents the maximum value among the N data.

3. The method according to claim 2, characterized in that, Step 5 includes: Step 5-1: Convert the clustered color image into a cell binarized image: Count the cells for each class. First, for the first class of cells after clustering, binarize the first class of cells in the clustered color image as white, and binarize the remaining points other than the first class of cells as black to obtain the cell binarized image. Step 5-2: Perform distance transformation on the binarized cell image to obtain the distance-transformed binary cell image; Step 5-3: Use the findContours() contour extraction algorithm in OpenCV to extract all contours of the distance-transformed binary cell image and obtain the number of contours. Step 5-4: Use the contourArea algorithm in OpenCV to calculate the area of ​​all contours; Step 5-5: If the area of ​​a contour is greater than the threshold th1, proceed to step 5-6; otherwise, directly return the number of contours as the counting result for each type of cell. Steps 5-6: Cell counting based on density images; Steps 5-7: For cells in other clusters after clustering, repeat steps 5-1 to 5-6 until all cells in all clusters have been counted.

4. The method according to claim 3, characterized in that, In step 5-1, the cell binarization image is obtained using the following formula: Where j = 1, 2, ..., N, B j [m,n]=255 means that the pixel value at the pixel with x-coordinate m and y-coordinate n in the binarized image of the j-th cell is 255, where m=1,2,...,W, n=1,2,...,H, W is the width of the binarized cell image, and H is the height of the binarized cell image.

5. The method according to claim 4, characterized in that, Step 5-2 includes: The distance transformation is to calculate the distance from a non-zero pixel in an image to the nearest zero pixel, as shown in the following formula: Among them, B j [m,n] represents the non-zero pixels in the binary image of the j-th cell type, B j [p,q] represents the zero pixel in the binary image of the j-th cell type, B j [p,q]=0 indicates that the pixel value at the x-coordinate p and y-coordinate q in the binary image of the j-th cell type is 0; D(B j [m,n]=255,B j [p,q]=0)] represents the non-zero pixel B in the binary image of the j-th cell type. j [m,n] and zero pixel B j The distance between [p, q], P j [m,n] represents the pixel value at the pixel with x-coordinate m and y-coordinate n in the binary image after distance transformation; min(·) represents the minimum value function, and D(P,Q) represents the Euclidean distance between points P and Q. th2 is the distance threshold between non-zero pixels and zero pixels in the binarized image of the j-th type of cell.

6. The method according to claim 5, characterized in that, Steps 5-6 include: Step 5-6-1: Binarize the image P of the j-th cell type. j [m,n] is convolved with the Gaussian filtered template image G[u,v] to obtain the cell density image D of the j-th type of cells. j [m,n]: Where u represents the x-coordinate of the Gaussian filter template image, v represents the y-coordinate of the Gaussian filter template image, σ is the standard deviation of the Gaussian image, and e represents the natural constant. This is a convolution operation; Step 5-6-2: Convert the density image of cell type j to D j The sum of the values ​​at all pixels in [m,n] gives the number of cells of type j.

7. A flow cytometry system for cell classification and counting based on density estimation, characterized in that, The system is used to implement the method as described in any one of claims 1 to 6, and the system includes a data acquisition unit, a preprocessing unit, a dimensionality reduction unit, a judgment unit, a clustering unit, and a counting unit; The data acquisition unit is used to acquire flow cytometry data: when cells enter the flow cytometer, after passing through multicolor fluorescence, they are collected by an optical system and converted into electrical signals by a photoelectric sensor. The electrical signals are converted into digital signals by an A / D analog-to-digital converter. The digital signals are called feature data. The data acquisition unit stores all the feature data of cells or particles acquired by the flow cytometer, which is called the raw flow cytometry data file. The raw flow cytometry data file is stored in the dataset. The preprocessing unit is used to standardize the dataset, and then send the processed dataset into the dimensionality reduction unit. The dimensionality reduction unit uses the UMAP algorithm to automatically reduce the dimensionality of the dataset; The judgment unit uses the SPSS chi-square test to determine whether the data in the dataset follows a uniform distribution. When the data distribution is uniform, the clustering unit uses the DBSCAN algorithm to automatically cluster the cell distribution; otherwise, the clustering unit uses mean-shift clustering to automatically cluster the cell distribution. The counting unit is used to count the cells in each clustered cell image separately.

Images