An xgboost-based cell microdroplet identification method
By using XGBoost machine learning and the RankMSE metric, a cell-empty droplet binary classification model was constructed, which solved the problem of distinguishing empty droplets from cells in single-cell RNA sequencing and achieved cell identification with higher accuracy and stability.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- 张浩
- Filing Date
- 2023-10-09
- Publication Date
- 2026-07-21
AI Technical Summary
Existing technologies struggle to effectively distinguish between real cells and empty droplets in single-cell RNA sequencing, leading to data analysis errors, particularly in the selection of UMI counting thresholds, making it difficult to simultaneously retain cells with low RNA content and filter out empty droplets.
We employed the XGBoost machine learning method, combined with the RankMSE index and multiple rounds of downsampling, to construct a cell-empty droplet binary classification model. Through iterative training, we optimized the model and improved the identification accuracy by utilizing cell quality control features and gene expression levels.
Robust cell identification was achieved across different datasets, effectively excluding empty droplets and low-quality cells, improving identification accuracy and stability, and reducing the false positive rate.
Smart Images
Figure CN117316296B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of single-cell RNA sequencing technology, specifically to a cell droplet identification method based on XGBoost. Background Technology
[0002] Single-cell RNA sequencing has gradually developed into a mature technology in recent years. In particular, the emergence of microfluidic droplet-based single-cell sequencing technology has dramatically increased the cell throughput of single-cell experiments. In these droplet-based single-cell capture methods, to ensure that most cells are individually encapsulated, a much higher number of droplets than cells are needed. This results in most droplets being empty droplets without cells, ultimately producing data containing millions of cells mixed with empty droplets. Ideally, these empty droplets do not contain any cellular RNA molecules. However, under the shear stress of the fluid, during apoptosis, cell fragmentation, or secretion, cytoplasmic RNA, mitochondrial RNA, ribosomes, or other cellular debris flow out of the droplets and into the ambient solution. This RNA present throughout the reaction environment is called ambient RNA. After entering the solution, ambient RNA diffuses over time and gradually enters the originally RNA-free empty droplets, resulting in a large number of empty droplet data points containing unique molecular identifiers (UMIs) in the final single-cell data, making them indistinguishable directly from cells. Furthermore, due to factors such as amplification, the UMI count of some empty droplets can even exceed that of some cells with low RNA content, making them even more difficult to distinguish from real cell data. Researchers often choose a UMI count threshold to screen cells; however, this single-dimensional classification presents two problems: 1. An excessively high threshold can filter out most empty droplets, but it will also lose cells with low RNA content; 2. An excessively low threshold, while obtaining most cells, will also introduce too many empty droplets, leading to erroneous analytical conclusions. Therefore, for droplet-based single-cell data, identifying real cells from a large number of empty droplets still presents a computational challenge. Therefore, an XGBoost-based cell droplet identification method is proposed. Summary of the Invention
[0003] To address the shortcomings of existing technologies, this invention provides a cell droplet identification method based on XGBoost, which solves the problems mentioned in the background section.
[0004] To achieve the above objectives, the present invention provides the following technical solution: a method for identifying cell droplets based on XGBoost, comprising the following steps:
[0005] Step 1: Use the RankMSE metric to search for cell and empty droplet boundaries and construct a training set, and solve the class imbalance problem in the training data through multiple rounds of downsampling;
[0006] Step 2: Combining pre-calculated cell quality control features such as cell gene expression levels and cell entropy, construct a cell-empty drop binary classification model suitable for the current data using the machine learning method XGBoost;
[0007] Step 3: The cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data through iterations to obtain an optimized prediction model;
[0008] Step 4: Identify cells using a cell-empty drop binary classification model.
[0009] Optionally, in step one: using the RankMSE index to search for the boundaries between cells and empty droplets and construct a training set, and solving the class imbalance problem in the training data through multiple rounds of downsampling, RankMSE is first calculated based on nCount_rank and nFeature_rank, and all droplets in the data are predefined using the RankMSE curve, mainly defined as four categories: cells, uncertain, empty droplets, and discarded droplets. The predefined standard for cells is all droplets before the first trough, the uncertain is defined as all droplets between the first trough and the peak of the second peak, the empty droplets are defined as all droplets from the peak of the second peak to the second trough, and the remaining droplets are discarded droplets with extremely low RNA content but extremely large data volume.
[0010] Optionally, in step one: using the RankMSE metric to search for the boundaries between cells and empty droplets and construct a training set, and solving the class imbalance problem in the training data through multiple rounds of downsampling, multiple rounds of UMI downsampling are performed on the predefined labeled cells and uncertain droplets to simulate cells and uncertain droplets at the empty droplet sequencing depth, and the originally limited number of droplets is expanded to be equal to or even greater than the number of empty droplets. Finally, these downsampled cell droplets and predefined empty droplet data are used as the first round of training set, and the standardized gene expression matrix is put into XGBoost for training.
[0011] Optionally, in step three: the cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data in an iterative manner to obtain an optimized prediction model. After the first round of training, the model mainly predicts and scores the uncertain microdroplets and empty droplets simulated in step two, and finally classifies them according to the set threshold. The empty droplet data in the training set is readjusted, the empty droplet data in the previous round is filtered and retained, and the data newly predicted as empty droplets from uncertain microdroplets in this round is added to the training set. This iterative training is repeated until the training error reaches the minimum, and finally an optimized model is obtained.
[0012] Optionally, in step three: the cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data through iterations to obtain an optimized prediction model. All data, including simulated droplet data, are put into the iteratively optimized model for the final round of prediction. The prediction scores of cells and droplets are obtained by averaging multiple simulated droplets. Finally, the droplets are comprehensively scored by combining global quality control indicators such as the cell Gini index. The droplet data is classified according to a threshold, which is 0.9 by default. After the classification is completed, the number of false positives is estimated, and cell droplets with lower comprehensive scores are reclassified as uncertain droplets based on the estimated number.
[0013] Optionally, in step three: the cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data in an iterative manner to obtain an optimized prediction model. Before training the model, indicators such as the proportion of mitochondrial content, the proportion of ribosome content, cell entropy, and cell redundancy were calculated for the droplet data and added to the feature set that originally only contained genes.
[0014] This invention provides a method for identifying cell droplets based on XGBoost, which has the following beneficial effects:
[0015] This XGBoost-based cell droplet identification method, compared to the existing four commonly used methods—CellRangerV2, CellRangerV3, EmptyDrops, and zUMIs—has the advantage of being able to robustly identify cells in different datasets and can effectively exclude empty droplets, low-quality cells, or cell debris from the data, thus exhibiting higher accuracy and stability compared to existing technologies. Attached Figure Description
[0016] Figure 1 This is a schematic diagram of the steps and structure of the present invention;
[0017] Figure 2 This is a flowchart of the present invention. Detailed Implementation
[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.
[0019] Please see Figures 1 to 2 This invention provides a technical solution: a method for identifying cell droplets based on XGBoost, comprising the following steps:
[0020] Step 1: Use the RankMSE metric to search for cell and empty droplet boundaries and construct a training set, and solve the class imbalance problem in the training data through multiple rounds of downsampling;
[0021] First, RankMSE is calculated based on nCount_rank and nFeature_rank. Then, all droplets in the data are predefined using the RankMSE curve, primarily categorized into four types: cells, uncertain droplets, empty droplets, and discarded droplets. Cells are defined as all droplets before the first trough; uncertain droplets are defined as all droplets between the first trough and the peak of the second peak; empty droplets are defined as all droplets from the peak of the second peak to the second trough; and the remaining droplets are discarded droplets with extremely low RNA content but a large data volume. dropSplit will not perform any operations on discarded droplets. To prevent incorrect judgments when searching for the first trough, dropSplit first calculates the inflection point ranks of nCount and nFeature, and automatically selects the trough around the smaller of the two inflection point ranks.
[0022] Among them, RankMSE:
[0023] In a droplet-based single-cell experimental dataset, if the UMI count of droplet g is C g The non-zero gene expression count is F. g C g and F g The ranks in all droplets are as follows: and Then define the RankSE of the droplet g:
[0024]
[0025] Due to RankSE g The index assigned to each droplet cannot reflect the continuous trend of its Rank before and after, so we assign all droplets according to R. C After sorting in ascending order, a sliding window method is used to smooth RankSE: a sliding window of length w is selected from the i-th droplet G. i Slide backwards one by one, and calculate the average of the RankSE of the first w droplets as the droplet G. i RankMSE i :
[0026]
[0027] Cell entropy, cell redundancy, and cell Gini index:
[0028] For a droplet, if it has n genes with non-zero UMI counts, the UMI counts of these genes form a vector A = {x1, x2, x3, ..., x...} n If these genes are assigned to UMI, then the probability is given by the vector P = {p1, p2, p3, ..., p}. n},in The cell entropy of this droplet is defined as:
[0029]
[0030] The cell Gini index of this droplet is defined as:
[0031]
[0032] The total UMI count of the droplet is m = ∑x. If these m UMIs are evenly distributed among n genes and the cell entropy is recalculated, the theoretical maximum cell entropy ExpectedMax(CellEntropy) of the droplet can be obtained, and the cell redundancy is defined as:
[0033]
[0034] Based on the vector P after uniform distribution according to UMI, the theoretical maximum cell Gini index ExpectedMax(CellGini) can be obtained by recalculation. We define the standardized cell Gini index as:
[0035] Normlize(CellGini)=CellGini+1-ExpectedMax(CellGini) (6)
[0036] The cell Gini index in this method uses the standardized value by default.
[0037] This method calculates the inflection point for two indicators: total UMI count and total non-zero gene count. Taking nCount as an example, when the total UMI count nCount of a given droplet g is calculated... g After the value, nCount g Sort the data in descending order based on size and calculate their rank. At this point, we can set the function log(nCount) = f(log(R)). C To make the curve sufficiently smooth, the `smooth.spline` function in R is used, with the degrees of freedom parameter set to 10, to fit the curve. Since inflection point detection on a logarithmic scale is easily affected by the scale level of the inflection point, the fitted smooth curve is exponentially transformed back to a normal scale. At this point, the function is `fitted(nCount) = f(fitted(R...)`. CBy calculating the first derivative f' and second derivative f'” at each point under the function curve, the curvature at each point is finally obtained:
[0038]
[0039] The inflection point of the curve is where the curvature is greater than 0 and the absolute value is the largest.
[0040] In the above method, the method for detecting the peak or valley of the curve and the method for peak selection include: given an arbitrary sequence S and the required left and right shoulder widths l and r for detecting the peak, the position i of the peak in the sequence needs to satisfy: After determining the peak position i in the sequence S, return the peak position i. If the number of positions n that satisfy the peak condition is greater than 1, then return all positions I = {i1, i2, ..., i...}. n If the negative number is taken before the given sequence, it becomes a process of detecting valley positions in the sequence, and the method is the same as above.
[0041] In a single detection, if multiple peaks exist, the first round of peak selection begins. From front to back, the difference *a* between two peaks is detected in pairs, and the difference *b* between the minimum value of the sequence between the two peaks and the peak value of the second highest peak is detected in pairs. Given a tolerance value *t*, if... If the second highest peak of the two peaks is retained, then only the highest peak of the two peaks is retained. If the number of peaks is still greater than 1 after all peak positions have been retained and filtered, then a second round of peak selection is performed on the remaining peaks until only one peak exists. The position of this peak in the sequence is the final selected peak position.
[0042] Multiple rounds of UMI downsampling were performed on predefined data labeled as cells and uncertain droplets to simulate cells and uncertain droplets at the empty droplet sequencing depth. The originally limited number of droplets was expanded to be equal to or even greater than the number of empty droplets. Finally, these downsampled cell droplets and predefined empty droplet data were used as the first training set. The standardized gene expression matrix was put into XGBoost for training. The significance of downsampling here is to correct the sequencing depth between droplets so that they are truly comparable in the feature space of genes. Multiple rounds of downsampling can simulate the distribution of cells and uncertain droplets at low sequencing depth and can solve the class imbalance problem in the training data.
[0043] Step 2: Combining pre-calculated cell quality control features such as cell gene expression levels and cell entropy, construct a cell-empty drop binary classification model suitable for the current data using the machine learning method XGBoost;
[0044] Step 3: The cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data through iterations to obtain an optimized prediction model;
[0045] After the first round of training, the model primarily predicts and scores the uncertain droplets and empty droplets simulated in the previous step, and finally classifies them according to a set threshold. Then, the empty droplet data in the training set is readjusted; the empty droplet data from the previous round is filtered and retained, and data newly predicted as empty droplets from uncertain droplets in the current round is added to the training set. This iterative training process is repeated until the training error is minimized, ultimately yielding an optimized model. In essence, this model fits most of the empty droplet data that can be distinguished from cell data in the original empty droplet dataset. By gradually incorporating similar empty droplet data from uncertain droplets, the final model achieves optimal classification ability for high-quality empty droplet data in uncertain regions. Furthermore, since the model is binary, its ability to classify cells in uncertain regions is also correspondingly improved.
[0046] All data, including simulated droplet data, are fed into the iteratively optimized model for the final round of prediction. The prediction scores for cells and droplets are obtained by averaging multiple simulated droplets. Finally, global quality control indicators such as the cell Gini index are combined to give the droplets a comprehensive score. The droplet data are classified according to a default threshold of 0.9. In order to control false positives caused by excessive model training errors or excessively low threshold settings, the number of false positives is estimated after classification. Based on the estimated number, cell droplets with lower comprehensive scores are reclassified as uncertain droplets.
[0047] The number of false positives is estimated as follows: Generally, the higher the dropSplitScore, the more likely the droplets are cells. If the classification error rate during model training is e∈(0,0.5), and the dropSplitScore threshold is T∈(0.5,1), the false positive rate r in the final classification result is estimated as follows:
[0048] r = 1 - (1 - e) × T (8)
[0049] Step 4: Identify cells using a cell-empty drop binary classification model;
[0050] Before training the model, indicators such as the proportion of mitochondrial content, the proportion of ribosome content, cell entropy, and cell redundancy were calculated for the droplet data and added to the feature set that originally only contained genes. The calculation of these indicators helps to improve the accuracy of the model. In addition, this method provides visualization options for relevant quality control indicators during the calculation and training process, so as to monitor the training process and make timely adjustments.
[0051] verify:
[0052] In multiple simulated datasets, the high-quality PBMC dataset, and the ESC cell line dataset with high background noise, the cell identification capabilities of our method dropSplit, along with four commonly used methods—CellRangerV2, CellRangerV3, EmptyDrops, and zUMIs—were tested and evaluated. The results show that dropSplit can robustly identify cells in different datasets and effectively exclude empty droplets, low-quality cells, or cell debris from the data, demonstrating higher accuracy and stability compared to other methods.
[0053] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.
Claims
1. A method for identifying cell droplets based on XGBoost, characterized in that: Includes the following steps: Step 1: Use the RankMSE metric to search for cell and empty droplet boundaries and construct a training set, and solve the class imbalance problem in the training data through multiple rounds of downsampling; Step 2: Combining the cell gene expression level and cell entropy pre-calculated cell quality control features, construct a cell-empty drop binary classification model suitable for the current data using the machine learning method XGBoost; Step 3: The cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data through iterations to obtain an optimized prediction model; Step 4: Identify cells using a cell-empty drop binary classification model; Step 1: Use the RankMSE index to search for the boundaries between cells and empty droplets and construct a training set. Then, solve the class imbalance problem in the training data through multiple rounds of downsampling. First, calculate RankMSE based on nCount_rank and nFeature_rank. Then, predefine all droplets in the data using the RankMSE curve and define them as four categories: cells, uncertain, empty droplets, and discarded droplets. The predefined criteria for cells are all droplets before the first trough. The uncertain droplets are defined as all droplets between the first trough and the peak of the second peak. The empty droplets are defined as all droplets from the peak of the second peak to the second trough. The remaining droplets are discarded droplets with extremely low RNA content but a large amount of data. dropSplit will not perform any operations on the discarded droplets; to prevent incorrect judgments when searching for the first valley, dropSplit first calculates the inflection point ranks of nCount and nFeature, and automatically selects the valley around the smaller inflection point rank. Among them, RankMSE: In a single-cell experimental dataset based on droplets, if droplets The UMI count is The non-zero gene expression count is , and The ranks in all droplets are as follows: and Then define a droplet RankSE: （1） because The indicators assigned to each droplet cannot reflect the continuous trend of its rank change. Therefore, all droplets are classified according to... After sorting in ascending order, the sliding window method is used to process the data. Perform smoothing: Select a length of The sliding window from the first A tiny droplet Slide backwards one by one and calculate the previous one. A tiny droplet The average value of the droplets of : （2） Cell entropy, cell redundancy, and cell Gini index: For a droplet, if its UMI count is non-zero, then the number of genes present is... These genes' UMI counts constitute a vector. The probability that these genes are assigned to UMI is then a vector. ,in The cell entropy of this droplet is defined as: （3） The cell Gini index of this droplet is defined as: （4） The total UMI count of the droplet is If these m UMIs are evenly distributed to By recalculating the cell entropy from each gene, the theoretical maximum cell entropy of the droplet can be obtained. And define cell redundancy as: （5） Based on the vector after uniform distribution in UMI Recalculation yields the theoretical maximum cellular Gini index. We define the standardized cell Gini index as: （6） The cell Gini index in this method uses the standardized value by default. This method calculates the inflection point for two indicators: total UMI count and total non-zero gene count. Taking nCount as an example, when a given droplet... Total UMI count After the value, Sort the data in descending order based on size and calculate their rank. At this point, a function can be set. To make the curve sufficiently smooth, the `smooth.spline` function in R is used, with the degrees of freedom parameter set to 10, to fit the curve. Since inflection point detection on a logarithmic scale is easily affected by the scale level of the inflection point, the fitted smooth curve is exponentially transformed back to a normal scale. At this point, the function is... By calculating the first derivative at each point under the curve of this function and second derivative Finally, the curvature of each point is obtained: （7） The inflection point of the curve is where the curvature is greater than 0 and the absolute value is the largest. In the above methods, the methods for detecting curve peaks or valleys and selecting peaks include: given a set of arbitrary numbers The required left and right shoulder widths for the detection peak and The location of the peak in the sequence Then it needs to satisfy: After evaluating all elements in the sequence, return the sequence. Peak position If the number of positions that meet the peak condition Then return all positions. If the negative number is taken before the given sequence, it becomes a process of detecting valley positions in the sequence, and the method is the same as above. In a single detection, if multiple peaks are present, the first round of peak selection begins, sequentially detecting the difference between two peaks in pairs, from front to back. And the difference between the minimum value of the sequence between the two peak positions and the peak value of the second peak. When a tolerance value is given ,like If the peaks are selected, the second highest peak of the two peaks is retained; otherwise, only the highest peak of the two peaks is retained. If the number of peaks is still greater than 1 after all peak positions have been selected, a second round of peak selection is performed on the remaining peaks until only one peak exists. The position of this peak in the sequence is the final selected peak position. Multiple rounds of UMI downsampling were performed on predefined labeled cells and uncertain droplets to simulate cells and uncertain droplets at empty droplet sequencing depth. The originally limited number of droplets was expanded to be equal to or even greater than the number of empty droplets. These downsampled cell droplets and predefined empty droplet data were used as the first round of training set. The standardized gene expression matrix was put into XGBoost for training. The significance of downsampling is to correct the sequencing depth between droplets so that they are truly comparable in the feature space of genes. Multiple rounds of downsampling can simulate the distribution of cells and uncertain droplets at low sequencing depth and can solve the class imbalance problem in the training data.
2. The method for identifying cell droplets based on XGBoost according to claim 1, characterized in that: Step 3: The cell-empty droplet binary classification model is retrained by iteratively adding newly predicted empty droplet data to obtain an optimized prediction model. After the first round of training, the uncertain microdroplets and empty droplets simulated in Step 2 are predicted and scored, and finally classified according to the set threshold. The empty droplet data in the training set is readjusted, the empty droplet data in the previous round is filtered and retained, and the data newly predicted as empty droplets from uncertain microdroplets in this round is added to the training set. This iterative training is repeated until the training error reaches the minimum, and finally an optimized model is obtained.
3. The method for identifying cell droplets based on XGBoost according to claim 1, characterized in that: Step three: The cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data through iterations to obtain an optimized prediction model. All data, including simulated droplet data, are put into the iteratively optimized model for the final round of prediction. The prediction scores of cells and droplets are obtained by averaging multiple simulated droplets. Finally, the droplets are comprehensively scored by combining the cell Gini index global quality control index. The droplet data is classified according to a threshold, which is 0.9 by default. After the classification is completed, the number of false positives is estimated, and cell droplets with lower comprehensive scores are reclassified as uncertain droplets based on the estimated number.
4. The method for identifying cell droplets based on XGBoost according to claim 1, characterized in that: Step 3: The cell-empty droplet binary classification model is retrained by adding newly predicted empty droplet data in an iterative manner to obtain an optimized prediction model. Before training the model, the proportion of mitochondrial content, proportion of ribosome content, cell entropy, and cell redundancy index of microdroplet data are calculated and added to the feature set that originally only contained genes.