Malignant cell separation method and system based on transcriptome sequencing
Through a multi-level threshold and correlation analysis method based on transcriptome sequencing, the accuracy and sensitivity issues of malignant cell identification are solved, and efficient and accurate malignant cell separation is achieved. It is applicable to various tumor single-cell transcriptome data, including primary tumors, metastatic tumors and circulating tumor cells.
Patent Information
- Application Number
- CN202510873980.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-27
- Publication Date
- 2025-10-17
AI Technical Summary
Existing malignant cell identification methods have suboptimal accuracy and sensitivity in single-cell sequencing data analysis, making it difficult to effectively identify malignant tumor cells in the tumor microenvironment, especially in hepatocellular carcinoma (HCC), due to the high heterogeneity and difficulty in processing noisy data.
A transcriptome sequencing-based approach was used to perform copy number variation (CNV) analysis using the R package inferCNV. A multi-level threshold strategy and correlation analysis were combined to identify malignant cells. This approach included data preprocessing, CNV deviation calculation, patient- and sample-level threshold determination, reference correlation analysis, and independence verification, and generated visual charts.
It significantly improves the accuracy of malignant cell identification, adapts to sample heterogeneity, objectively determines thresholds, enhances the credibility and repeatability of results, realizes a highly automated analysis process, and provides intuitive visual results.
Smart Images

Figure CN120808871A_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application relates to a malignant cell separation method and system based on transcriptome sequencing, belonging to the cross field of bioinformatics, tumor biology and computational biology. BACKGROUND
[0002] With the rapid development of single-cell sequencing technology, researchers can analyze the gene expression profiles of various cell types in tissues at the single-cell level. In the field of tumor research, single-cell sequencing technology provides an unprecedented opportunity to understand tumor heterogeneity, tumor microenvironment and tumor evolution. Hepatocellular carcinoma (HCC) is the fourth leading cause of cancer death worldwide, and its tumor microenvironment is complex and diverse, containing malignant tumor cells, immune cells, stromal cells and other cell types.
[0003] However, in single-cell sequencing data analysis, accurate identification and separation of malignant tumor cells still face many challenges: first, the differences between tumor cells and normal cells at the transcriptome level are often not obvious; second, tumor cells themselves are highly heterogeneous; third, existing cell type annotation methods are mainly based on known marker genes, and have limited effect on identifying malignant cells with complex genomic variations.
[0004] Currently, the commonly used malignant cell recognition methods in the industry include: recognition based on the expression pattern of specific tumor markers, inference methods based on copy number variation (CNV), and detection methods based on gene mutations. Among them, the CNV-based method is widely used because it can capture the genomic instability characteristics of malignant cells, but existing methods still have limitations in handling high-noise data, setting reasonable thresholds, and integrating multi-dimensional information, resulting in suboptimal accuracy and sensitivity of malignant cell recognition.
[0005] With the development of precision medicine for tumors, there is an increasing demand for accurate separation of malignant cells in the tumor microenvironment. Accurate identification of malignant cells not only helps to reveal the mechanisms of tumor evolution, but also provides important basis for the development of targeted therapy and immunotherapy strategies. Therefore, developing an efficient and accurate malignant cell separation method has important scientific value and clinical application prospects. SUMMARY
[0006] The purpose of the present application is to provide a malignant cell separation method and system based on transcriptome sequencing, which uses single-cell transcriptome data to perform copy number variation (CNV) analysis and malignant cell identification through the R package inferCNV.
[0007] To achieve the above-mentioned purpose, the technical solution adopted by the present application is as follows:
[0008] On the one hand, the present application provides a malignant cell separation method based on transcriptome sequencing,
[0009] Step 1, reading and preprocessing transcriptome sequencing data;
[0010] Step 2, calculating window CNV based on gene position;
[0011] Step 3, calculating patient-level and sample-level CNV thresholds respectively;
[0012] Step 4, performing reference correlation analysis, including: reference cell determination, cell correlation calculation, correlation threshold determination;
[0013] Step 5, integrating CNV threshold and correlation for cell classification;
[0014] Step 6, independence verification;
[0015] Step 7, generating analysis results and visualization charts.
[0016] Further, the step 1 includes:
[0017] (1) reading the inferCNV prediction matrix I ∈ Rm×n through the data interface, m is the number of genes ≥ 20000 to ensure genome coverage, n is the number of cells ≥ 1000 to ensure statistical reliability, and reading cell annotation information A containing cell type and tissue source label, and gene position information G containing gene coordinate mapping;
[0018] (2) the gene position information G contains the chromosome coordinate mapping of each gene g:
[0019] Map(g i )→(Chr k ,Pos start ,Pos end )
[0020] Where Pos is in base pairs bp, representing the physical start and end positions of the gene on the chromosome;
[0021] (3) sorting the genes in the inferCNV matrix according to the chromosome physical position, establishing the gene position index Map(g)→P, where P=(Chr, Pos) represents the physical position of gene g on the chromosome;
[0022] (4) calculating the deviation D ij of each gene in each cell relative to the normal expression value:
[0023] D ij =I ij -median(I i,normal )
[0024] Where median(I_i,normal) represents the median expression value of gene i in normal cells. This step realizes the normalization of CNV signal based on normal cell reference.
[0025] Furthermore, the step 2 includes:
[0026] (1) Setting a fixed-size gene window W;
[0027] (2) Slide the window along the chromosome position and calculate the CNV deviation mean D of all genes in each window w , where D w Reflects the degree of copy number variation in the w-th window:
[0028]
[0029] where i∈window w represents all genes in the window w, and c represents cells;
[0030] (3) Calculate the CNV score of cell c to quantify the overall intensity of copy number variation across the genome:
[0031]
[0032] Where N is the total number of windows. The square operation amplifies the significant variation signal and suppresses the background noise, so that the CNV score can distinguish malignant cells from non-malignant cells.
[0033] Furthermore, in step 3, the patient-level CNV threshold is calculated as follows:
[0034] (1) Integrate the CNV score data of all tissue samples of the patient to generate a patient-level CNV score vector Vp∈R^{np}, where np≥5000 cells;
[0035] (2) Calculate the skewness S of Vp p :
[0036]
[0037] Where μ and σ are the mean and standard deviation of Vp respectively;
[0038] (3) According to the skewness value S p Determine distribution characteristics;
[0039] (4) Apply the inflection point detection algorithm to determine the patient-level CNV score threshold T p :
[0040] (4.1) Construct the normalized cumulative distribution curve Γ(x);
[0041] (4.2) Calculate the vertical distance d from each point on the curve to the line L connecting the beginning and end, which satisfies:
[0042]
[0043] where (xi, yi) is the coordinate of the i-th point on the curve, (x1, y1) and (xn, yn) are the coordinates of the beginning and end points, respectively;
[0044] (4.3) Take argmax(d) as the threshold point T p ;
[0045] (5) Based on the threshold T p , identify the high CNV cell set:
[0046] H p = {c | Score CNV (c) > T p}.
[0047] Further, in step 3, the sample-level CNV threshold is calculated as follows:
[0048] (1) Only using the CNV score data of the current sample, generate a sample-level CNV score vector Vs ∈ R^{ns}, where ns is the number of cells in the sample;
[0049] (2) Calculate the skewness S s of Vs, and use the same skewness calculation formula as in claim 4 to determine the distribution characteristics;
[0050] (3) Apply the same inflection point detection algorithm as the patient-level to determine the sample-level CNV score threshold Ts, and refer to the algorithm steps in claim 4 for step
[0051] (4) Based on the threshold T s , identify the high CNV cell set within the sample:
[0052] H s = {c | Score CNV (c) > T s}.
[0053] Further, in step 4, the reference cell determination includes:
[0054] (1) Integrate the CNV data of all samples of the patient to generate a global CNV score vector Vg ∈ R^{ng}, where ng is the total number of cells in the patient;
[0055] (2) Calculate the P-th percentile of Vg as the reference threshold T r :
[0056] T r= percentile(V g , 100 - P)
[0057] where P is the percentile parameter, and the default value is 0.5;
[0058] (3) Select cells with CNV scores higher than T r as the reference cell set R:
[0059] R = {c | Score CNV (c) > T r};
[0060] The cell correlation calculation includes:
[0061] (1) Calculate the average CNV pattern vector M R of the reference cells:
[0062]
[0063] where Mc represents the CNV pattern vector of cell c, and |R| represents the number of reference cells;
[0064] (2) Calculate the Spearman correlation coefficient p R of each cell c with M c :
[0065]
[0066] where d i represents the rank difference of the CNV pattern of cell c and M R in the ith window, and N is the total number of windows;
[0067] (3) Generate the correlation score vector CorrS ∈ R^n for all cells:
[0068] The correlation threshold determination includes:
[0069] (1) Select cells with CNV scores lower than min(Ts, Tp) as the benchmark population B:
[0070] B = {c | Score CNV (c) < min(T s , T p )}
[0071] These cells represent the non-malignant background, which is used to determine the correlation background distribution;
[0072] (2) Perform multi-method normal distribution fitting on the correlation distribution of cells in B, and calculate the mean μB and standard deviation σB;
[0073] (3) Calculate the correlation threshold Tcorr as the criterion for separating malignant cells:
[0074] T corr =μ B +k·σ B
[0075] where k = 2.57 corresponds to a 99% confidence interval;
[0076] (4) Identify the high-correlation cell set Hcorr based on the threshold Tcorr:
[0077] H corr ={c|ρ c >T corr}.
[0078] Furthermore, the step 5 includes:
[0079] (1) Identify the high CNV cell set Hp based on the patient-level CNV threshold Tp;
[0080] (2) Identify the high-correlation cell set Hcorr based on the correlation threshold Tcorr;
[0081] (3) Take the intersection of the two to determine the malignant cell set M:
[0082] M=H p ∩H corr
[0083] That is, cells that simultaneously meet the requirements of high CNV score and high correlation are judged as malignant cells;
[0084] (4) Calculation of the malignant cell ratio P M :
[0085]
[0086] Where |M| represents the number of malignant cells, and n is the total number of cells in the current sample;
[0087] (5) Generate cell classification result data table C, which contains information such as cell ID, cell type label (malignant / non-malignant), CNV score, correlation score, and tissue source, providing basic data for downstream analysis.
[0088] Furthermore, the independence verification in step 6 includes distribution independence analysis, density peak analysis and comprehensive judgment, wherein:
[0089] The distribution independence analysis includes:
[0090] (1) Extract the correlation score vectors CorrS_M and CorrS_N of the malignant cell set M and the non-malignant cell set N respectively, where CorrS_M = {c1 M ,c2 M ,...,c m M} represents the correlation scores of m malignant cells, and CorrS_N = {c1 N ,c2 N ,...,c n N} represents the correlation scores of n non-malignant cells;
[0091] (2) Calculate Hartigan's Dip Test statistic D and p-value:
[0092] [D,p] = dip.test(Corr_S)
[0093] where the D statistic measures the maximum deviation of the empirical distribution function from unimodal distribution, and when p < 0.05, the unimodal hypothesis is rejected at a 95% confidence level, and it is determined as a traditional multimodal distribution;
[0094] (3) Calculate the normalized mean difference NMD:
[0095]
[0096] where mean(CorrS_M) and mean(CorrS_N) represent the mean of the correlation scores of malignant and non-malignant cells respectively, and sd(CorrS) represents the overall standard deviation of all cell correlation scores, used to standardize the mean difference between different samples;
[0097] (4) Calculate the overlap coefficient OVL of the two types of cell distribution by kernel density estimation:
[0098] OVL = ∫min(f M (x),f N (x))dx × 100%
[0099] where fM and fN represent the kernel density estimation functions of malignant and non-malignant cells respectively, the integral interval is the joint support set of the two distributions, and the OVL value range is [0, 1], 0 represents complete non-overlap, and 1 represents complete overlap;
[0100] (5) Calculate the Bhattacharyya distance BD:
[0101] BD = 0.25 × ln(0.25 × ((σ 2 M / σ 2 N)+(σ 2 N / σ2 M)+2))+0.25x((μM-μN) 2 / (σ 2 M+σ 2 N))
[0102] where μMand μNare the mean values of the malignant and non-malignant cell related scores, σ 2 Mand σ 2 Nare the corresponding variances, the larger the BD value, the more significant the difference between the two distributions;
[0103] (6) Fitting unimodal and bimodal Gaussian mixture models, and calculating the BIC difference ΔBIC:
[0104] ΔBIC = BIC1 - BIC2
[0105] where BIC1 represents the Bayesian information criterion of the unimodal Gaussian model, and BIC2 represents the Bayesian information criterion of the bimodal Gaussian mixture model, ΔBIC > 10 strongly supports the bimodal model, and ΔBIC > 5 moderately supports the bimodal model;
[0106] (7) Calculate the relative density ratio RDR:
[0107] RDR = fM(μM) / max(fN(μM), 1x10-10)
[0108] where fM(μM) represents the nuclear density value of malignant cells at their mean value μM, fN(μM) represents the nuclear density value of non-malignant cells at the mean value μM of malignant cells, and 1x10-10 in the denominator is a small constant to prevent zero, RDR > 2 indicates that malignant cells form a relatively dense distribution in their aggregation area;
[0109] (8) Calculate the local background ratio LBR:
[0110] LBR = ∫[Q5, Q 95 ]fM(x)dx / max(∫[Q5, Q 95 ]fN(x)dx, 1x10 -10 )
[0111] where Q5and Q 95 represent the 5% and 95% quantiles of the malignant cell distribution, respectively, defining the main distribution interval of malignant cells, and the numerator and denominator represent the density integrals of malignant and non-malignant cells in this interval, respectively, LBR > 2 indicates that malignant cells are significantly enriched in their main distribution area relative to background noise;
[0112] (9) Peak separation degree calculation:
[0113] Identify all local maxima of the density curve using kernel density estimation, get the array of peak positions peaks_positions, and calculate the minimum of the distance between adjacent peaks:
[0114] Peak_Separation = min(diff(peaks_positions)) / sd(CorrS_all)
[0115] where diff(peaks_positions) calculates the array of differences between adjacent peak positions, diff(peaks_positions) = [peak2-peak1, peak3-peak2,..., peak N -peak n-1 ], Peak_Separation measures the degree of separation between peaks in units of standard deviation, >1 indicates that there is a clear separation between peaks;
[0116] (10) Enhanced multimodality score:
[0117] Adopt the conditional cumulative scoring mechanism of multi-dimensional evidence fusion:
[0118] S_enhanced = min((∑score_conditions / 15)×10, 10)
[0119] where 15 is the theoretical maximum sum of each condition score, and 10 is the upper limit of the final score, and the rules for each condition score are:
[0120] BIC evidence: ΔBIC>50 plus 3 points; 20<ΔBIC≤50 plus 2 points; 10<ΔBIC≤20 plus 1 point
[0121] Peak separation evidence: >2σ plus 3 points; 1σ<separation≤2σ plus 2 points; 0.5σ<separation≤1σ plus 1 point, where σ is sd(CorrS_all)
[0122] Relative density evidence: RDR>10 plus 3 points; 5<RDR≤10 plus 2 points; 2<RDR≤5 plus 1 point
[0123] Local background evidence: LBR>10 plus 3 points; 5<LBR≤10 plus 2 points; 2<LBR≤5 plus 1 point
[0124] Traditional multimodality evidence: Hartigan's Dip Test passes (p<0.05) plus 3 points
[0125] When S_enhanced≥5, it is determined to be an enhanced multimodal distribution;
[0126] The density peak analysis comprises:
[0127] (1) Calculate the kernel density estimation f(x) of the correlation distribution:
[0128]
[0129] where K is the standard Gaussian kernel function K(u)=(1 / √2π)exp(-u 2 / 2), n is the total number of observations, h is the bandwidth parameter, the base bandwidth is determined using the Silverman rule, then an adjustment factor adjust=0.8 is set to improve sensitivity to small peaks, the final bandwidth h_final=h×adjust, the number of sampling points is set to 1024 to ensure high-resolution density estimation, xi is the ith observation;
[0130] (2) Implement a dynamic threshold adjustment strategy:
[0131] relative_position=(x-x_min) / (x_max-x_min)
[0132] position_gain=exp(2×relative_position) / exp(2)
[0133] adaptive_threshold=base_threshold / position_gain
[0134] final_threshold=max(adaptive_threshold,min_threshold_limit)
[0135] x_min and x_max are the minimum and maximum values of the correlation score, respectively, relative_position∈[0,1] represents the relative position of the position in the distribution range, exp(2)≈7.389 is a normalization constant, base_threshold=0.03×max(f(x)) is the base height threshold, min_threshold_limit=0.01×max(f(x)) is the lower limit protection to prevent the threshold from being too low, this strategy ensures that the high correlation area obtains a lower detection threshold to improve sensitivity; (3) Implement additional threshold optimization for high correlation areas:
[0136] high_corr_bonus=1+(x-mean(CorrS_all)) / range(CorrS_all)
[0137] optimized_threshold = adaptive_threshold / high_corr_bonus
[0138] where range(CorrS_all) = max(CorrS_all) - min(CorrS_all) is the range of the correlation scores, high_corr_bonus > 1 is the gain factor for high correlation regions, and mean(CorrS_all) is the mean of all cell correlation scores. This optimization provides more sensitive detection of high correlation regions where malignant cells tend to cluster.
[0139] (4) Detect local maximum points p that satisfy the following conditions as peaks:
[0140] f(p) > f(p-1) AND f(p) > f(p+1) AND f(p) > tfinai(p) AND Prom(p) > min_prominence(p)
[0141] where min_prominence(p) is the minimum prominence threshold, with a base value of max(f(x)) x 0.01, and is dynamically adjusted according to the position gain function.
[0142] (5) Calculate the peak prominence Prom(p):
[0143] Prom(p) = peak_height - max(min_left, min_right)
[0144] where peak_height is the density value at the peak, min_left and min_right are the minimum density values found to the left and right of the peak, respectively. Prominence reflects the significance of the peak relative to the surrounding "peaks and valleys". Prom(p) > 0.01 x max(f(x)) indicates that the peak is sufficiently prominent.
[0145] (6) Calculate the peak-to-valley ratio Dist:
[0146]
[0147] where p1 and p2 are the position coordinates of two peaks, respectively, and σ represents the standard deviation of the overall distribution. Dist is expressed in units of standard deviation, facilitating comparison between different samples. Dist > 0.5 indicates that the peaks are distinguishably separated.
[0148] (7) Calculate the peak-to-valley ratio PVR:
[0149]
[0150] where valley_min represents the lowest density value between two peaks, PVR > 1.3 means the peak is well separated from the background;
[0151] (8) Calculate long-distance micro-peak reward DB:
[0152] distance_score = |secondary_peak - main_peak| / sd(CorrS_all)
[0153] height_ratio = secondary_peak_height / main_peak_height
[0154] where main_peak is the position of the highest density peak, secondary_peak is the position of the secondary peak, main_peak_height and secondary_peak_height are the corresponding density heights, distance_score represents the normalized peak distance, height_ratio represents the height ratio of the secondary peak to the main peak, and the scoring rules are as follows: when distance_score > 2, height_ratio > 0.5, 4 points, 0.3-0.5, 3 points, 0.1-0.3, 2 points, ≤0.1, 1 point; when 1 < distance_score ≤ 2, reduce 1 point according to the above rules;
[0155] (9) High correlation region peak reward:
[0156] corr_position_bonus = min(2, (x - mean(CorrS_all)) / range(CorrS_all) x 5)
[0157] where x is the peak position, mean(CorrS_all) is the overall mean, and range(CorrS_all) is the value range. The peak located in the high correlation region (x > mean(CorrS_all)) and on the right side of the main peak can obtain an additional corr_position_bonus score reward, up to 2 points;
[0158] (10) Bimodality comprehensive score:
[0159] Bimodality_Score = min(∑feature_scores, 10)
[0160] where second_to_main_prominence_ratio=prominence_secondary / prominence_main represents the ratio of the second highest peak to the highest peak, and feature_scores include inter-peak distance score (2 points for normalized_distance>0.8, 1 point for >0.5), peak-valley separation score (2 points for valley_ratio>1.3, 1 point for >1.1), long-distance peak reward (calculated according to distance_score and height_ratio, maximum 4 points), high correlation region reward (maximum 2 points), secondary peak prominence score (3 points for second_to_main_prominence_ratio>0.3, 2 points for >0.2, 1 point for >0.1), multiple peak number reward (1 point for detecting >2 peaks), distribution asymmetry reward (1 point for |skewness|>0.5), and the presence of a clear second peak is determined when Bimodality_Score≥4;
[0161] The comprehensive judgment adopts a hierarchical priority decision mechanism:
[0162] (1) Highest priority decision:
[0163] When Bimodality_Score≥4, the determination result is "obvious second peak exists", indicating that the density peak value analysis detects a significant second peak;
[0164] (2) Second highest priority decision:
[0165] When S_enhanced≥5 and simultaneously satisfy NMD>1.5, OVL<0.3, ΔBIC>10, and BD>0.5, the determination result is "distribution analysis confirms independence";
[0166] (3) Medium priority decision:
[0167] When simultaneously satisfy NMD>1.0, OVL<0.5, ΔBIC>5, and BD>0.3, the determination result is "possible independence";
[0168] (4) Default decision:
[0169] The determination result in other cases is "unimodal distribution", indicating that no evidence of independence of malignant cells is detected;
[0170] (5) Confidence quantitative evaluation:
[0171] Confidence_Score = min(10, max(0, 0.5 x (Bimodality_Score / 10) x 5 + 0.2 x (S_enhanced / 10) x 2 + 0.3 x min(ABIC / 10, 3)))
[0172] In the formula, the peak analysis contributes 50% weight (coefficient 0.5 x 5), the enhanced multimodality analysis contributes 20% weight (coefficient 0.2 x 2), and the BIC model comparison contributes 30% weight (coefficient 0.3, and the upper limit of the contribution of ABIC is 3 points), and the final score range is [0, 10];
[0173] (6) Verification report generation:
[0174] The output includes a complete verification report containing the final determination category, confidence score, key statistical indicators and decision basis, providing quantitative reliability evaluation for downstream analysis.
[0175] Further, the step 7 comprises:
[0176] (1) Generating a cell classification summary statistics table S, containing the number and proportion of each type of cell;
[0177] (2) Generating a CNV heat map H to show the CNV patterns of different types of cells, the behavior window, and the cells as columns, and the color representing the CNV deviation intensity;
[0178] (3) Generating a correlation distribution graph D to show the distribution difference between malignant and non-malignant cells, including the nuclear density curve, threshold line and peak value label;
[0179] (4) Generating a verification result visualization chart V, containing the multimodality analysis and peak detection results;
[0180] (5) Saving the analysis results as a structured data file R, containing cell classification results, statistical indicators and verification conclusions.
[0181] On the other hand, the present application also provides a malignant cell separation system based on transcriptome sequencing, which is used to realize the above-mentioned method, and the system comprises:
[0182] A data input interface module for receiving single cell inferCNV prediction matrix, cell annotation information and gene position information;
[0183] A preprocessing engine for performing gene ordering and CNV deviation calculation;
[0184] A CNV calculation unit for implementing sliding window calculation and CNV quantity quantification;
[0185] A multi-level threshold processor, comprising a patient-level threshold calculation subunit and a sample-level threshold calculation subunit;
[0186] a correlation analysis module, comprising a reference cell selector, a correlation calculator, and a threshold determiner;
[0187] Cell classification decision unit, which performs dual-standard judgment based on CNV score and correlation;
[0188] Independence verification system, including distribution analyzer and peak detector;
[0189] Result output and visualization module, used to generate analysis reports and charts.
[0190] The present invention has the following beneficial effects:
[0191] (1) Improved classification accuracy: Through a multi-level threshold strategy and dual-standard judgment, the accuracy of malignant cell identification is significantly improved by considering both CNV scores and correlation.
[0192] (2) Adaptation to sample heterogeneity: The double-layer threshold design at the patient level and sample level can simultaneously maintain consistency in cross-sample analysis and adaptability to the characteristics of individual samples, effectively solving the technical problem of high heterogeneity of tumor samples.
[0193] (3) Objective threshold determination: The automatic threshold calculation method based on inflection point detection and normal distribution fitting avoids the deviation caused by subjective settings and improves the repeatability and reliability of the analysis results.
[0194] (4) Comprehensive verification mechanism: Through two complementary methods, distribution independence analysis and density peak analysis, the independence and real existence of malignant cell populations are comprehensively verified, which enhances the credibility of the results.
[0195] (5) High degree of automation: The entire analysis process is highly automated, with no need for human intervention from data input to result output, which greatly improves analysis efficiency and reduces operational complexity.
[0196] (6) Visual presentation of results: Provides a variety of visual charts to intuitively display the analysis process and results, making it easier for researchers to interpret and evaluate.
[0197] (7) Wide range of applications: This method is applicable to single-cell transcriptome data of various types of tumors, including primary tumors, metastatic tumors, and circulating tumor cells, and has broad application prospects. BRIEF DESCRIPTION OF THE DRAWINGS
[0198] Figure 1 This is a flow chart of the malignant cell separation method based on multi-level thresholds of the present invention;
[0199] Figure 2 Flowchart for multi-tier CNV score threshold determination;
[0200] Figure 3 Flowchart for reference correlation analysis;
[0201] Figure 4 Flowchart for cell classification decision;
[0202] Figure 5 Flowchart for independence verification system;
[0203] Figure 6 Scatter plot for malignant and non-malignant cell CNV score-correlation distribution;
[0204] Figure 7 Heat map for chromosome CNV pattern of single-sample malignant cell separation result;
[0205] Figure 8 Distribution independence analysis and peak density detection result plot for independence verification system. DETAILED DESCRIPTION
[0206] The present application will be described in detail below in conjunction with the accompanying drawings and specific embodiments.
[0207] Embodiment One
[0208] A malignant cell separation method based on transcriptome sequencing, comprising the following physical implementation steps:
[0209] S1: Data reading and preprocessing
[0210] To solve the technical problems of high noise and discontinuous gene positions in the inferCNV result matrix of single-cell transcriptome data, this step is implemented through the following physical processing operations:
[0211] (1) Read the inferCNV prediction matrix (I ∈ Rm×n, m is the number of genes ≥ 20000, n is the number of cells ≥ 1000), cell annotation information (A), and gene position information (G) through the data interface;
[0212] (2) The gene position information G contains the chromosome coordinate mapping of each gene g:
[0213] Map(g i )→(Chr k ,Pos start ,Pos end )
[0214] Where Pos is in base pairs (bp) and represents the physical start and end positions of the gene on the chromosome;
[0215] (3) Sort genes in the inferCNV matrix according to their chromosomal physical positions, and build a gene position index Map(g)→P, where P=(Chr, Pos) represents the physical position of gene g on the chromosome;
[0216] (4) Calculate the deviation of each gene in each cell relative to the normal expression value Dij(unit: relative expression deviation value):
[0217] D ij = I ij -median(I i,normal )
[0218] where median(I_i,normal) represents the median expression value of gene i in normal cells, and this step realizes the CNV signal standardization based on normal cell reference.
[0219] S2: Window CNV calculation
[0220] To solve the problem of gene-level noise interference and chromosome region CNV detection, the following sub-steps are implemented:
[0221] (1) Set a fixed size of gene window (W=100 genes), where W represents the number of continuous genes covering the chromosomal region, and the size corresponds to a physical distance of about 10-15 Mb (1.0-1.5×10^7 base pairs) in the human genome, which is sufficient to capture common cancer-related CNV events (such as EGFR amplification of about 0.5-2 Mb) and greater than the noise-related length scale (<5 Mb) of transcriptome sequencing technology;
[0222] (2) Slide the window along the chromosomal position, and calculate the CNV deviation mean (Dw, unit: average relative expression deviation value) of all genes in each window, where Dw reflects the copy number variation degree of the wth window:
[0223]
[0224] where i∈window w represents all genes in window w, and c represents the cell;
[0225] (3) Calculate the CNV score (Scor eCNV (c), unit: square deviation sum) of cell c, which quantifies the overall intensity of the whole genome copy number variation:
[0226]
[0227] Where N is the total number of windows, square operation amplifies significant variant signals, suppresses background noise, and enables CNV score to effectively distinguish between malignant cells (high score) and non-malignant cells (low score).
[0228] S3: Multi-level threshold determination
[0229] To solve the problems of sample heterogeneity and threshold objective determination, this step is realized through a double threshold strategy:
[0230] S3.1: Patient-level CNV threshold calculation
[0231] To establish a consistent classification standard across samples, this sub-step is realized by the following operations:
[0232] (1) Integrate the CNV score data of all tissue samples of the patient to generate a patient-level CNV score vector Vp∈R^{n_p}, where np≥5,000 cells;
[0233] (2) Calculate the skewness (Skewness) Sp (dimensionless statistic) of Vp:
[0234]
[0235] Where μ and σ are the mean and standard deviation of Vp, respectively;
[0236] (3) Determine the distribution characteristics according to the skewness value Sp: |Sp|>2 indicates a highly skewed distribution that needs special processing;
[0237] (4) Apply the inflection point detection algorithm to determine the patient-level CNV score threshold T_p (unit: same as CNV score):
[0238] ① Construct a normalized cumulative distribution curve Γ(x);
[0239] ② Calculate the vertical distance d of each point on the curve to the connecting line L between the beginning and end, which satisfies:
[0240]
[0241] ③ Take argmax(d) as the threshold point Tp;
[0242] (5) Identify the high CNV cell set based on the threshold Tp:
[0243] H p = {c | Score CNV (c) > T p}
[0244] S3.2: Sample-level CNV threshold calculation
[0245] To adapt to the single sample-specific distribution, this sub-step is implemented by the following operations:
[0246] (1) Using only the CNV score data of the current sample, generate a sample-level CNV score vector Vs∈R^{n_s}, where ns is the number of cells in the sample;
[0247] (2) Calculate the skewness Ss of Vs to determine the distribution characteristics;
[0248] (3) Apply the same inflection point detection algorithm as the patient level to determine the sample-level CNV score threshold Ts (unit: same as CNV score);
[0249] Based on the threshold Ts, identify the high CNV cell set within the sample:
[0250] H s ={c∣Score CNV (c)>T s}。
[0251] S4: Reference correlation analysis
[0252] To solve the problem of malignant cell CNV pattern recognition, this step is implemented by the following physical processing operations:
[0253] S4.1: Reference cell determination
[0254] To establish a reliable reference CNV pattern, this sub-step is implemented by the following operations:
[0255] (1) Integrate the CNV data of all samples of the patient to generate a global CNV score vector Vg∈R^{n_g}, where ng is the total number of cells in the patient;
[0256] (2) Calculate the P percentile (default P = 0.5%) of Vg as the reference threshold Tr (unit: same as CNV score):
[0257] T r =percentile(V g ,100-P)
[0258] This percentage is selected based on the biological characteristics that highly malignant cells usually account for <1% observed in tumor heterogeneity research;
[0259] (3) Select cells with CNV scores higher than Tr as the reference cell set R:
[0260] R={c∣Score CNV (c)>T r}
[0261] These cells represent the most significant cell population for CNV pattern, to establish the reference pattern of malignant cells.
[0262] S4.2: Cell correlation calculation
[0263] To quantify the similarity of each cell to the malignant pattern, this sub-step is implemented by:
[0264] (1) Calculate the average CNV pattern vector M of reference cells R ∈R^N (unit: average relative expression deviation value):
[0265]
[0266] Where Mc represents the CNV pattern vector of cell c, |R| represents the number of reference cells, and N is the total number of windows;
[0267] (2) Calculate the Spearman correlation coefficient p_c (dimensionless, range [-1, 1]) of each cell c with M R
[0268]
[0269] Where d_i represents the rank difference of the CNV pattern of cell c and M_R at the i-th window, and N is the total number of windows;
[0270] (3) Generate the correlation score vector Corr_S ∈ R^n for all cells.
[0271] S4.3: Correlation threshold determination
[0272] To objectively determine the correlation classification boundary, this sub-step is implemented by:
[0273] (1) Select cells with CNV scores lower than min(T_s, T_p) as the reference population B:
[0274] B = {c | Score CNV (c) < min(T s , T p )}
[0275] These cells represent the non-malignant background, used to determine the correlation background distribution;
[0276] (2) Perform multi-method normal distribution fitting on the correlation distribution of cells in B, and calculate the mean μB and standard deviation σB:
[0277] μ1, σ1 = mean(Corr_S_B), sd(Corr_S_B)
[0278] μ2,σ2 = fitdistr(Corr_S_B, "normal")
[0279] μ3,σ3 = fitdist(Corr_S_B, "norm")
[0280] μ_B = median(μ1,μ2,μ3)
[0281] σ_B = median(σ1,σ2,σ3)
[0282] (3) Calculate correlation threshold T_corr (dimensionless, range [-1, 1]):
[0283] T corr = μ B +k·σ B
[0284] where k = 2.57 corresponding to 99% confidence interval (statistical significance level a = 0.01), which selection is based on strict requirement for false positive control in tumor biology;
[0285] (4) Identify high correlation cell set H_corr based on threshold T_corr:
[0286] H corr = {c | p c > T corr}.
[0287] S5: Cell classification
[0288] To achieve robust malignant cell identification, this step is achieved by multi-criteria intersection determination:
[0289] (1) Identify high CNV cell set Hp according to patient-level CNV threshold Tp;
[0290] (2) Identify high correlation cell set Hcorr according to correlation threshold Tcorr;
[0291] (3) Take the intersection of the two to determine the malignant cell set M:
[0292] M = H p ∩H corr
[0293] That is, cells that meet both high CNV score and high correlation are determined as malignant cells;
[0294] (4) Calculate the proportion of malignant cells P M (percent):
[0295]
[0296] Where n is the total number of cells;
[0297] (5) Generate cell classification result data table C, which contains the type label and related statistical data of each cell, providing a basis for downstream analysis.
[0298] S6: Independence Verification
[0299] To evaluate the independence and true existence of the inferred malignant cell population, this step is implemented through three sub-steps: distribution independence analysis, density peak analysis, and comprehensive judgment:
[0300] S6.1: Distributional Independence Analysis
[0301] To quantify the degree of separation between malignant and non-malignant cell distribution, this substep is achieved by the following operations:
[0302] (1) Data preprocessing:
[0303] Extract the correlation score vectors CorrS_M and CorrS_N of the malignant cell set M and the non-malignant cell set N respectively, where CorrS_M={c1 M ,c2 M ,...,c m M} represents the correlation score of m malignant cells, CorrS_N={c1 N ,c2 N ,...,c n N} represents the correlation score of n non-malignant cells;
[0304] (2) Hartigan's Dip multimodality test:
[0305] Calculate Hartigan's Dip Test statistic D and p-value, where the D statistic measures the maximum degree to which the empirical distribution function deviates from the unimodal distribution. When p < 0.05, the unimodal hypothesis is rejected at the 95% confidence level, and the distribution is judged to be a traditional multimodal distribution.
[0306] (3) Calculation of standardized mean difference:
[0307] NMD=|mean(CorrS_M)-mean(CorrS_N)| / sd(CorrS_all)
[0308] where mean(CorrS_M) and mean(CorrS_N) represent the mean of correlation scores of malignant and non-malignant cells respectively, sd(CorrS_all) represents the population standard deviation of correlation scores of all cells, for normalizing the mean difference between different samples, NMD > 1.5 means the difference of mean correlation between two types of cells is more than 1.5 standard deviations, with significant separation;
[0309] (4) Distribution overlap coefficient calculation:
[0310] OVL = ∫min(fM(x), fN(x))dx
[0311] where fM and fN represent the kernel density estimation functions of malignant and non-malignant cells respectively, the integral interval is the joint support set of two distributions, the value range of OVL is [0, 1], 0 represents completely no overlap, 1 represents completely overlap, OVL < 0.3 means the overlap degree of two types of cells is low, with good separation;
[0312] (5) Bhattacharyya distance calculation:
[0313] BD = 0.25 x ln(0.25 x ((σ 2 M / σ 2 N)+(σ 2 N / σ 2 M)+2))+0.25 x ((μM-
[0314] μN) 2 / (σ 2 M+σ 2 N))
[0315] where μM and μN are the mean of correlation scores of malignant and non-malignant cells respectively, σ 2 M and σ 2 N are the corresponding variances, the larger the BD value represents the more significant difference between two distributions, BD > 0.5 means two distributions have significant difference;
[0316] (6) Mixed Gaussian model comparison:
[0317] ΔBIC = BIC_single - BIC_double
[0318] where BIC_single represents the Bayesian information criterion of single peak Gaussian model, BIC_double represents the Bayesian information criterion of double peak Gaussian mixed model, ΔBIC > 10 strongly supports the double peak model, ΔBIC > 5 moderately supports the double peak model;
[0319] (7) Relative density ratio calculation:
[0320] RDR = fM(μM) / max(fN(μM), 1 x 10-10)
[0321] where fM(μM) represents the nuclear density value of malignant cells at their mean μM, fN(μM) represents the nuclear density value of non-malignant cells at the mean μM of malignant cells, 1 x 10-10 in the denominator is a tiny constant to prevent zero, and RDR > 2 indicates that malignant cells form a relatively dense distribution in their aggregation area;
[0322] (8) Local background ratio calculation:
[0323] LBR = ∫[Q5, Q 95 ] fM(x) dx / max(∫[Q5, Q 95 ] fN(x) dx, 1 x 10-10) -10
[0324] where Q5 and Q 95 represent the 5th and 95th percentiles of the malignant cell distribution, respectively, defining the main distribution interval of malignant cells, and the numerator and denominator represent the density integrals of malignant and non-malignant cells in this interval, respectively, and LBR > 2 indicates that malignant cells are significantly enriched in their main distribution area relative to background noise;
[0325] (9) Peak separation calculation:
[0326] Use kernel density estimation to identify all local maxima of the density curve to obtain the peak position array peaks_positions, and calculate the minimum value of the distance between adjacent peaks:
[0327] Peak_Separation = min(diff(peaks_positions)) / sd(CorrS_all)
[0328] where diff(peaks_positions) calculates the difference array of adjacent peak positions, Peak_Separation measures the separation degree between peaks in standard deviation units, and > 1 indicates that there is a clear separation between peaks;
[0329] (10) Enhanced multi-peak score:
[0330] Adopt a conditional cumulative scoring mechanism for multi-dimensional evidence fusion:
[0331] S_enhanced = min((∑score_conditions / 15) x 10, 10)
[0332] where 15 is the theoretical maximum sum of condition scores, 10 is the upper limit of the final score, and the rules for each condition score are: BIC evidence (ΔBIC > 50 plus 3 points; 20 < ΔBIC ≤ 50 plus 2 points; 10 < ΔBIC ≤ 20 plus 1 point), peak separation evidence (> 2σ plus 3 points; 1σ < separation ≤ 2σ plus 2 points; 0.5σ < separation ≤ 1σ plus 1 point), relative density evidence (RDR > 10 plus 3 points; 5 < RDR ≤ 10 plus 2 points; 2 < RDR ≤ 5 plus 1 point), local background evidence (LBR > 10 plus 3 points; 5 < LBR ≤ 10 plus 2 points; 2 < LBR ≤ 5 plus 1 point), traditional multimodality evidence (Hartigan's Dip Test passed plus 3 points), and when S_enhanced ≥ 5, it is determined to be an enhanced multimodal distribution;
[0333] S6.2: Density peak analysis
[0334] To detect the obvious peaks in the correlation distribution, this sub-step is implemented by the following operations:
[0335] (1) High-resolution kernel density estimation:
[0336] f(x) = (1 / nh) x ∑ i K((x-xi) / h)
[0337] where K is the standard Gaussian kernel function K(u) = (1 / √2π)exp(-u 2 / 2), n is the total number of observations, h is the bandwidth parameter, the base bandwidth is determined using the Silverman rule, then an adjustment factor adjust = 0.8 is set to improve the sensitivity to small peaks, and the final bandwidth h_final = h x adjust, the number of sampling points is set to 1024 to ensure high-resolution density estimation, and xi is the ith observation;
[0338] (2) Dynamic threshold adjustment strategy:
[0339] relative_position = (x-x_min) / (x_max-x_min)
[0340] position_gain = exp(2 x relative_position) / exp(2)
[0341] adaptive_threshold = base_threshold / position_gain
[0342] final_threshold = max(adaptive_threshold, min_threshold_limit)
[0343] where x_min and x_max are the minimum and maximum of the correlation score respectively, relative_position ∈ [0, 1] represents the relative position of the location within the distribution range, exp(2) ≈ 7.389 is the normalization constant, base_threshold = 0.03 x max(f(x)) is the base height threshold, min_threshold_limit = 0.01 x max(f(x)) is the lower limit protection to prevent the threshold from being too low, and this strategy ensures that the high correlation area obtains a lower detection threshold to improve sensitivity;
[0344] (3) High correlation area threshold optimization:
[0345] high_corr_bonus = 1 + (x - mean(CorrS_all)) / range(CorrS_all)
[0346] optimized_threshold = adaptive_threshold / high_corr_bonus
[0347] where range(CorrS_all) = max(CorrS_all) - min(CorrS_all) is the range of the correlation score, high_corr_bonus ≥ 1 is the gain coefficient of the high correlation area, and mean(CorrS_all) is the mean of all cell correlation scores. This optimization provides more sensitive detection capability for the high correlation area where malignant cells tend to gather;
[0348] (4) Peak identification conditions:
[0349] Detect local maximum points p that satisfy the following conditions as peaks:
[0350] f(p) > f(p-1) AND f(p) > f(p+1) AND f(p) ≥ τfinal(p) AND Prom(p) ≥ min_prominence(p)
[0351] where min_prominence(p) is the minimum prominence threshold, the base value is set to max(f(x)) x 0.01, and is dynamically adjusted according to the position gain function;
[0352] (5) Peak prominence calculation:
[0353] Prom(p) = peak_height - max(min_left, min_right)
[0354] where peak_height is the density value at the peak, min_left and min_right are the lowest density values searched on the left and right of the peak, prominence reflects the prominence of the peak relative to the surrounding "valleys", Prom(p) > 0.01 x max(f(x)) means that the peak is prominent enough;
[0355] (6) Peak distance calculation:
[0356] Dist = |p1-p2| / σ
[0357] where p1 and p2 are the position coordinates of the two peaks, and σ represents the standard deviation of the overall distribution, and the peak distance is expressed in standard deviation units, which is convenient for comparison between different samples, and Dist > 0.5 means that the peaks are distinguishably separated;
[0358] (7) Peak-valley separation ratio calculation:
[0359] PVR = min(peak1_height, peak2_height) / valley_min_height
[0360] where valley_min represents the lowest density value between the two peaks, and PVR > 1.3 means that the peaks are well separated from the background;
[0361] (8) Detection of remote micro-peak:
[0362] distance_score = |secondary_peak-main_peak| / sd(CorrS_all)
[0363] height_ratio = secondary_peak_height / main_peak_height
[0364] In the formula, main_peak is the position of the peak with the highest density, secondary_peak is the position of the secondary peak, main_peak_height and secondary_peak_height are the corresponding density heights, distance_score represents the normalized peak distance, and height_ratio represents the height ratio of the secondary peak to the main peak. The scoring rules are as follows: when distance_score > 2, height_ratio > 0.5 is 4 points, 0.3 < height_ratio ≤ 0.5 is 3 points, 0.1 < height_ratio ≤ 0.3 is 2 points, and ≤ 0.1 is 1 point; when 1 < distance_score ≤ 2, reduce 1 point according to the above rules;
[0365] (9) Peak reward in high correlation region:
[0366] corr_position_bonus = min(2, (x - mean(CorrS_all)) / range(CorrS_all) x 5)
[0367] where x is the peak position, mean(CorrS_all) is the population mean, and range(CorrS_all) is the range of values. Peaks located in the high correlation region (x > mean(CorrS_all)) and to the right of the main peak can receive an additional corr_position_bonus bonus of up to 2 points.
[0368] (10) Bimodality Score:
[0369] Bimodality_Score = min(∑feature_scores, 10)
[0370] where second_to_main_prominence_ratio = prominence_secondary / prominence_main represents the ratio of the prominence of the secondary peak to the main peak, and feature_scores include the inter-peak distance score (2 points for normalized_distance > 0.8, 1 point for > 0.5), the peak-valley separation score (2 points for valley_ratio > 1.3, 1 point for > 1.1), the distant peak bonus (calculated based on distance_score and height_ratio, up to 4 points), the high correlation region bonus (up to 2 points), the secondary peak prominence score (3 points for second_to_main_prominence_ratio > 0.3, 2 points for > 0.2, 1 point for > 0.1), the multiple peak number bonus (1 point for detecting > 2 peaks), the distribution asymmetry bonus (1 point for |skewness| > 0.5), and the bimodality score (up to 10 points). When Bimodality_Score ≥ 4, it is determined that there is a clear second peak.
[0371] S6.3: Comprehensive Judgment
[0372] To provide the final verification conclusion, this sub-step uses a hierarchical priority decision mechanism:
[0373] (1) Highest Priority Decision:
[0374] When Bimodality_Score ≥ 4, the decision result is "obvious second peak exists", indicating that the density peak value analysis detects a significant second peak.
[0375] (2) Second Highest Priority Decision:
[0376] When S_enhanced≥5 and simultaneously satisfy NMD>1.5, OVL<0.3, ΔBIC>10, BD>0.5, the determination result is "Distribution analysis confirms independence";
[0377] (3) Medium priority determination:
[0378] When simultaneously satisfy NMD>1.0, OVL<0.5, ΔBIC>5, BD>0.3, the determination result is "Possible independence";
[0379] (4) Default determination:
[0380] In other cases, the determination result is "Unimodal distribution", indicating that no evidence of independence of malignant cells is detected;
[0381] (5) Confidence quantitative evaluation:
[0382] Confidence_Score=min(10,max(0,0.5×(Bimodality_Score / 10)×5+0.2×(S_enhanced / 10)×2+0.3×min(ΔBIC / 10,3)))
[0383] In the formula, the peak analysis contributes 50% weight (coefficient 0.5×5), the enhanced multimodality analysis contributes 20% weight (coefficient 0.2×2), and the BIC model comparison contributes 30% weight (coefficient 0.3, and the upper limit of ΔBIC contribution is 3 points), and the final score range is [0, 10];
[0384] (6) Verification report generation:
[0385] Output a complete verification report containing the final determination category, confidence score, key statistical indicators and decision basis, providing quantitative reliability evaluation for downstream analysis.
[0386] S7: Result output and visualization
[0387] In order to facilitate result interpretation and analysis, this step is realized by the following operations:
[0388] (1) Generate cell classification summary statistics table S:
[0389] Including the number and proportion of each type of cell, the detection rate of malignant cells, and the summary of key statistical indicators;
[0390] (2) Generate CNV heat map H:
[0391] Show the CNV pattern of different types of cells, the behavior sliding window, the column as the cell, and the color represents the CNV deviation intensity, the malignant cell area is obviously distinguished from the non-malignant cell;
[0392] (3) Generate correlation distribution map D:
[0393] Show the distribution difference between malignant and non-malignant cells, including nuclear density curve, dynamic threshold line, detected peak label and overlapping area annotation;
[0394] (4) Generate verification result visualization chart V:
[0395] Comprehensive display of multimodal analysis results, peak detection details, confidence scores and decision basis;
[0396] (5) Save analysis results:
[0397] Save the analysis results as a structured data file R, including cell classification results, all statistical indicators, verification conclusions and confidence evaluation, for subsequent analysis and result reproduction.
[0398] The inflection point detection algorithm involved in the implementation:
[0399] To solve the threshold determination problem caused by the skew distribution of single cell data, the inflection point detection algorithm is based on the following principles:
[0400] (1) Linear transformation is used for data normalization:
[0401]
[0402] This transformation eliminates the dimensional effect and preserves the distribution form, making CNV score data of different scales comparable;
[0403] (2) LOESS algorithm is used for smoothing:
[0404] y smooth = loess(y', x', span = 0.2)
[0405] Where span = 0.2 means using 20% of the data points for local weighted regression, which effectively filters out noise while preserving the inflection point characteristics;
[0406] (3) Geometric distance calculation uses vector projection method:
[0407] Let the connecting line be vector AB, and the distance d of point P to AB satisfy:
[0408]
[0409] This distance quantifies the convexity of the curve, and the maximum distance point is the inflection point;
[0410] (4) Sensitivity parameterized threshold screening:
[0411] dthreshold =max(d)×sensitivitycandidates={P|d(P)>d threshold}
[0412] The default sensitivity is 0.9, which controls the strictness of candidate point screening;
[0413] (5) Derivative analysis aids judgment:
[0414] best_knee=argmax P∈candidates |f′(P left )-f′(P right )|
[0415] Select the point where the curvature change (derivative change) is most significant as the final inflection point;
[0416] (6) Multi-level fallback mechanism:
[0417] if|candidates|<1:
[0418] sensitivity = sensitivity - 0.1
[0419] goto step 4
[0420] When there are insufficient candidate points, the sensitivity is automatically reduced to ensure that valid results are always returned.
[0421] Reference cell selection in this implementation:
[0422] To solve the problem of establishing a reference CNV pattern, this reference cell selection scheme is based on the following principles:
[0423] (1) Global percentile threshold calculation:
[0424] T r =percentile(V g ,100-P)
[0425] The default value of P is 0.5%, which means that the top 0.5% cells with the highest CNV scores are selected as references;
[0426] (2) Tissue-specific analysis:
[0427]
[0428] Where Rt represents the reference cell set in tissue t, and Pt represents the proportion of reference cells in the tissue.
[0429] The correlation threshold calculation in this implementation is:
[0430] To solve the problem of objectively determining the boundaries of relevance classification, this solution is based on the following principles:
[0431] (1) Multi-method normal distribution fitting:
[0432] μ1,σ1=mean(Corr_S_B),sd(Corr_S_B)
[0433] μ2,σ2=fitdistr(Corr_S_B,"normal")
[0434] μ3,σ3=fitdist(Corr_S_B,"norm")
[0435] μ_B=median(μ1,μ2,μ3)
[0436] σ_B=median(σ1,σ2,σ3)
[0437] The results of the three methods are combined to improve the robustness of the fitting parameters;
[0438] (2)Configured standard deviation multiplier:
[0439] T corr =μ+k·σ
[0440] The default value of k is 2.57, which corresponds to a 99% confidence interval (α = 0.01), and can be adjusted to 1.96 (95%) or 3.09 (99.9%) as needed;
[0441] (3) Normality verification mechanism:
[0442] p AD =anderson.darling.test(Corr_S B )p KS =kolmogorov.smirnov.test(Corr_S B )p SW
[0443] =shapiro.wilk.test(Corr_S B )
[0444] When the p values were all > 0.05, the distribution was judged to meet the normality assumption;
[0445] (4) Shape metric analysis:
[0446]
[0447] Where |skewness|<0.5 and |kurtosis|<0.5 indicate that the distribution is close to normal;
[0448] (5) Q-Q plot visual verification:
[0449] qqplot(Corr_S_B, rnorm(n, mu, sigma))
[0450] Visual assessment of the fit quality by comparing quantile relationship of actual distribution and theoretical normal distribution.
[0451] Density peak analysis in this implementation:
[0452] To solve the problem of malignant cell population independence verification, this scheme is based on the following principles:
[0453] (1) Position adaptive threshold strategy:
[0454] relative_position = (x - x_min) / (x_max - x_min)
[0455] position_gain = exp(2 × relative_position) / exp(2)
[0456] adaptive_threshold = base_threshold / position_gain
[0457] This strategy makes the threshold decrease exponentially with the correlation position index. When the relative position is close to 1, the threshold can be reduced to about 14% (1 / 7.389) of the base threshold, ensuring high sensitivity detection in the malignant cell aggregation area.
[0458] (2) Double condition peak verification:
[0459] f(p) > f(p-1) AND f(p) > f(p+1) AND f(p) ≥ τfinal(p) AND Prom(p) ≥ min_prominence(p)
[0460] Simultaneously consider the local maximum condition, dynamic height threshold and minimum prominence requirement to ensure that the detected peak has both statistical significance and biological significance, where min_prominence is set to max(f(x)) × 0.01;
[0461] (3) Bidirectional peak valley ratio evaluation:
[0462] PVR = min(peak1_height, peak2_height) / valley_min_height
[0463] Search for local minimum density points from both sides of the peak, calculate the ratio of lower peak to valley to assess the degree of peak separation, PVR>1.3 indicates that the peak is well separated from the background, >1.1 indicates moderate separation;
[0464] (4) Distance peak detection enhancement:
[0465] distance_score=|secondary_peak-main_peak| / sd(CorrS_all)
[0466] height_ratio=secondary_peak_height / main_peak_height
[0467] According to the standard deviation unit distance and height ratio of the main peak, give the far small peak additional score, when distance_score>2 can get the highest 4 points reward, when 1<distance_score≤2, according to the proportion of decreasing, this mechanism is particularly beneficial to detect the separation of malignant cell peak relative to the main background peak;
[0468] (5) Multi-dimensional feature fusion scoring system:
[0469] Bimodality_Score=min(∑feature_scores,10)
[0470] The scoring system integrates 7 core indicators, including: peak distance score (normalized_distance>0.8 2 points), valley separation score (valley_ratio>1.3 2 points), distance peak reward (highest 4 points), high correlation area reward (highest 2 points), secondary peak prominence score (second_to_main_prominence_ratio>0.3 3 points), multi-peak number reward (>2 peaks 1 point), distribution asymmetry reward (|skewness|>0.5 1 point), when Bimodality_Score≥4 is judged as there is obvious second peak.
[0471] Distribution independence analysis in this implementation:
[0472] In order to comprehensively evaluate the statistical independence of malignant cell population, this scheme is based on the following principles:
[0473] (1) Hartigan's Dip Test multimodality test:
[0474] D=sup_x|F_n(x)-F_unimodal(x)|
[0475] The test is a non-parametric method to assess unimodality of a distribution, D statistic measures the maximum deviation of empirical distribution function from unimodal distribution, p<0.05 indicates significant deviation from unimodal hypothesis at 95% confidence level, contributes 3 points for traditional multimodality evidence;
[0476] (2) Normalized mean difference:
[0477] NMD = |mean(CorrS_M) - mean(CorrS_N)| / sd(CorrS_all)
[0478] Normalized mean difference considers the overall data dispersion, eliminates the scale difference between samples, NMD>1.5 indicates strong significant difference between the two types of cells, >1.0 indicates moderate significant difference, this index provides key evidence for strengthening multimodality score;
[0479] (3) Distribution overlap coefficient:
[0480] OVL = ∫min(fM(x), fN(x))dx
[0481] The overlap of the two types of cell distribution is accurately quantified by kernel density estimation, the integral of the minimum value of the two density functions is calculated, OVL<0.3 indicates good separation, <0.5 indicates moderate separation, this index directly reflects the separability of the two types of cells;
[0482] (4) Bayesian information criterion model comparison:
[0483] ABIC = BIC_single - BIC_double
[0484] ABIC>10 indicates that the double peak mixed model is significantly better than the single peak model, strongly supports the independent existence of malignant cell population, >5 indicates moderate support, this index provides the highest 3 points for BIC evidence score;
[0485] (5) Relative density ratio:
[0486] RDR = fM(μM) / max(fN(μM),1×10 -10 )
[0487] This index is particularly suitable for identifying small and dense populations, compares the density of the two types of cells at the mean position of malignant cells, RDR>10 can obtain 3 points, 5
[0488] (6) Local background ratio:
[0489] LBR = ∫[Q5,Q 95 ]fM(x)dx / max(∫[Q5,Q95 ]fN(x)dx,1×10 -10 )
[0490] The density integral ratio of the two types of cells in the main distribution area of malignant cells (5%-95% percentile range) is evaluated, LBR>10 is 3 points, 5
[0491] (7) Bhattacharyya distance:
[0492] BD=0.25*ln(0.25*((sigma 2 M / sigma 2 N)+(sigma 2 N / sigma 2 M)+2))+0.25*((mu
[0493] M-mu 2 N) / (sigma 2 M+sigma 2 N)
[0494] Considering the mean difference and variance difference of the two distributions, BD>0.5 indicates that the two distributions have significant differences, and >0.3 indicates moderate differences. This distance index provides an important quantitative basis for comprehensive judgment;
[0495] (8) Multidimensional evidence fusion mechanism:
[0496] S_enhanced=min((∑score_conditions / 15)*10,10)
[0497] The condition cumulative score mechanism is used to integrate BIC evidence, peak separation evidence, relative density evidence, local background evidence and traditional multi-peak evidence. The theoretical maximum total score of 15 points is mapped to a 10-point scoring system. S_enhanced>=5 is determined as a strong multi-peak distribution, which is the second highest priority decision for comprehensive judgment and provides a key basis.
[0498] Example two
[0499] A malignant cell separation system based on multi-level threshold, the system implementation architecture is as follows
[0500] Data input interface module, for receiving single cell inferCNV prediction matrix, cell annotation information and gene location information;
[0501] Preprocessing engine, performing gene ordering and CNV bias calculation;
[0502] CNV calculation unit, implementing sliding window calculation and CNV quantification;
[0503] The multi-level threshold processor includes a patient-level threshold subunit and a sample-level threshold subunit;
[0504] The correlation analysis module includes a reference cell selector, a correlation calculator, and a threshold determinator;
[0505] The cell classification decision unit performs a double standard judgment based on the CNV score and the correlation;
[0506] The independence verification system includes a distribution analyzer and a peak detector;
[0507] The result output and visualization module generates an analysis report and a chart.
[0508] The following experimental verification and data analysis are performed for the above-mentioned scheme:
[0509] We used single-cell sequencing data from 7 liver cancer patients from a public database for verification, each patient including at least one normal tissue-derived sample and one tumor tissue-derived sample. The research process is as follows:
[0510] The raw data is input into inferCNV for preliminary CNV prediction;
[0511] The inferCNV results are input into our improved algorithm for processing;
[0512] The algorithm automatically calculates the CNV score and the correlation score;
[0513] The dynamic threshold is determined based on patient-specific data;
[0514] The malignant cells and non-malignant cells are divided according to the threshold;
[0515] The verification algorithm is applied to confirm the separation effect.
[0516] The verification results show that:
[0517] Our method can automatically determine the best threshold suitable for each patient's sample ( Figure 6 );
[0518] The identified malignant cell population has similar CNV expression patterns ( Figure 7 );
[0519] These cells express significantly high levels of epithelial and hepatocyte markers;
[0520] Through the verification algorithm, the labeled malignant cells and other cells show a significant bimodal distribution ( Figure 8 ), confirming that the algorithm effectively separates the malignant cells.
[0521] These experimental results demonstrate that our improved method can more accurately deal with the tumor heterogeneity problem, achieving more accurate separation of malignant cells, especially in complex tumor environments where traditional fixed threshold methods fail.
[0522] The above shows and describes the basic principles, main features and advantages of the present application. Those skilled in the art should understand that the above examples do not limit the protection scope of the present application in any form, and any technical solutions obtained by equivalent replacement or the like fall within the protection scope of the present application. The parts not involved in the present application are the same as or can be realized by using the prior art.
Claims
1. A method for isolating malignant cells based on transcriptome sequencing, characterized in that: The steps include: Step 1: Read and preprocess transcriptome sequencing data; Step 2: Calculate window CNV based on gene position; Step 3: Calculate the CNV thresholds at the patient level and sample level respectively; Step 4: perform reference correlation analysis, including: determining reference cells, calculating cell correlation, and determining correlation thresholds; Step 5: Cell classification based on CNV threshold and correlation; Step 6: Independence verification; Step 7: Generate analysis results and visualization charts.
2. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: The step 1 comprises: (1) Read the inferCNV prediction matrix I∈Rm×n through the data interface, where m is the number of genes ≥20,000 to ensure genome coverage and n is the number of cells ≥1,000 to ensure statistical reliability. At the same time, read the cell annotation information A including cell type and tissue source labels, and the gene location information G including gene coordinate mapping; (2) Gene location information G contains the chromosome coordinate mapping of each gene g: Map(g i )→(Chr k ,Pos start ,Pos end ) Pos is in base pairs bp, which indicates the physical start and end position of the gene on the chromosome; (3) Sort the genes in the inferCNV matrix according to their physical positions on the chromosome and establish a gene position index Map(g)→P, where P = (Chr, Pos) represents the physical position of gene g on the chromosome; (4) Calculate the deviation D of each gene relative to the normal expression value in each cell ij : D ij =I ij -median(I i,normal ) Where median(I_i,normal) represents the median expression value of gene i in normal cells. This step realizes the normalization of CNV signal based on normal cell reference.
3. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: The step 2 includes: (1) Setting a fixed-size gene window W; (2) Slide the window along the chromosome position and calculate the CNV deviation mean D of all genes in each window w , where D w Reflects the copy number variation degree of the w-th window: where i∈window w represents all genes in the window w, and c represents cells; (3) Calculate the CNV score of cell c to quantify the overall intensity of copy number variation across the genome: Where N is the total number of windows. The square operation amplifies the significant variation signal and suppresses the background noise, so that the CNV score can distinguish malignant cells from non-malignant cells.
4. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: In step 3, the patient-level CNV threshold is calculated as follows: (1) Integrate the CNV score data of all tissue samples of the patient to generate a patient-level CNV score vector Vp∈R^{np}, where np≥5000 cells; (2) Calculate the skewness S of Vp p : Where μ and σ are the mean and standard deviation of Vp respectively; (3) According to the skewness value S p Determine distribution characteristics; (4) Apply the inflection point detection algorithm to determine the patient-level CNV score threshold T p : (4.1) Construct the normalized cumulative distribution curve Γ(x); (4.2) Calculate the vertical distance d from each point on the curve to the line L connecting the end and the end to satisfy: Where (xi,yi) is the coordinate of the i-th point on the curve, (x1,y1) and (xn,yn) are the coordinates of the first and last points respectively; (4.3) Take argmax(d) as the threshold point T p ; (5) Based on threshold T p Identify high CNV cell populations: H p ={c∣Score CNV (c)>T p }。 5. The method for isolating malignant cells based on transcriptome sequencing according to claim 4, characterized in that: In step 3, the sample-level CNV threshold is calculated as follows: (1) Using only the CNV score data of the current sample, generate a sample-level CNV score vector Vs∈R^{ns}, where ns is the number of cells in the sample; (2) Calculate the skewness S of Vs s , determine the distribution characteristics; (3) Determine the sample-level CNV score threshold Ts using the same inflection point detection algorithm as that used at the patient level; (4) Based on threshold T s Identify high CNV cell populations within a sample: H s ={c∣Score CNV (c)>T s }。 6. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: In step 4, determining the reference cells includes: (1) Integrate the CNV data of all samples of the patient to generate a global CNV score vector Vg∈R^{ng}, where ng is the number of all cells in the patient; (2) Calculate the P percentile of Vg as the reference threshold T r : T r =percentile(V g ,100-P) Where P is the percentile parameter, the default value is 0.5; (3) Select CNV scores higher than T r The cells are used as the reference cell set R: R={c∣Score CNV (c)>T r }; The cell correlation calculation includes: (1) Calculate the average CNV pattern vector M of the reference cell R : Where Mc represents the CNV pattern vector of cell c, |R| represents the number of reference cells; (2) Calculate c and M for each cell R Spearman correlation coefficient ρ c : where d i The CNV pattern of cell c is different from that of M R The rank difference in the i-th window, N is the total number of windows; (3) Generate the correlation score vector CorrS∈R^n of all cells; The correlation threshold determination includes: (1) Select cells with CNV scores lower than min(Ts,Tp) as the reference population B: B={c∣Score CNV (c)<min(T s ,T p )} These cells represent the nonmalignant background and are used to determine the background distribution of correlations; (2) Perform multi-method normal distribution fitting on the correlation distribution of cells in B and calculate the mean μB and standard deviation σB; (3) Calculate the correlation threshold Tcorr as the criterion for separating malignant cells: T corr =μ B +k·s B where k = 2.57 corresponds to a 99% confidence interval; (4) Identify the high-correlation cell set Hcorr based on the threshold Tcorr: H corr ={c∣ρ c >T corr }。 7. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: The step 5 comprises: (1) Identify the high CNV cell set Hp based on the patient-level CNV threshold Tp; (2) Identify the high-correlation cell set Hcorr based on the correlation threshold Tcorr; (3) Take the intersection of the two to determine the malignant cell set M: M=H p ∩H corr That is, cells that simultaneously meet the requirements of high CNV score and high correlation are judged as malignant cells; (4) Calculation of the malignant cell ratio P M : Where |M| represents the number of malignant cells, and n is the total number of cells in the current sample; (5) Generate cell classification result data table C, which contains cell ID, cell type label, CNV score, correlation score and tissue source information, providing basic data for downstream analysis.
8. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: The independence verification in step 6 includes distribution independence analysis, density peak analysis and comprehensive judgment, wherein: The distribution independence analysis includes: (1) Extract the correlation score vectors CorrS_M and CorrS_N of the malignant cell set M and the non-malignant cell set N respectively, where CorrS_M = {c1 M ,c2 M ,...,c m M } represents the correlation score of m malignant cells, CorrS_N={c1 N ,c2 N ,...,c n N } represents the correlation score of n non-malignant cells; (2) Calculate Hartigan's Dip Test statistic D and p-value: [D,p]=dip.test(Corr_S) The D statistic measures the maximum degree to which the empirical distribution function deviates from the unimodal distribution. When p < 0.05, the unimodal hypothesis is rejected at the 95% confidence level and the distribution is determined to be a traditional multimodal distribution. (3) Calculate the standardized mean difference (NMD): Where mean(CorrS_M) and mean(CorrS_N) represent the means of the correlation scores of malignant and non-malignant cells, respectively, and sd(CorrS) represents the overall standard deviation of all cell correlation scores, which is used to standardize the mean differences between different samples; (4) Calculate the overlap coefficient OVL of the two types of cell distributions by kernel density estimation: OVL=∫min(f M (x),f N (x))dx×100% Where fM and fN represent the kernel density estimation functions of malignant and non-malignant cells, respectively. The integration interval is the joint support set of the two distributions. The OVL value range is [0,1], where 0 indicates no overlap and 1 indicates complete overlap. (5) Calculate the Bhattacharyya distance BD: BD=0.25×ln(0.25×((σ 2 M / s 2 N)+(σ 2 N / a 2 M)+2))+0.25×((μM-μN) 2 / (s 2 M+s 2 N)) where μM and μN are the mean of the correlation scores of malignant and non-malignant cells, respectively, and σ 2 M and σ 2 N is the corresponding variance, and the larger the BD value is, the more significant the difference between the two distributions is; (6) Fit the unimodal and bimodal Gaussian mixture models and calculate the BIC difference ΔBIC: ΔBIC=BIC1-BIC2 Where BIC1 represents the Bayesian information criterion for the unimodal Gaussian model, BIC2 represents the Bayesian information criterion for the bimodal Gaussian mixture model, ΔBIC>10 strongly supports the bimodal model, and ΔBIC>5 moderately supports the bimodal model; (7) Calculate the relative density ratio RDR: RDR=fM(μM) / max(fN(μM),1×10 -10 ) Where fM (μM) represents the nuclear density value of malignant cells at their mean μM, fN (μM) represents the nuclear density value of non-malignant cells at the mean μM of malignant cells, and the 1×10 -10 A very small constant to prevent division by zero; (8) Calculate the local background ratio LBR: LBR=∫[Q5,Q 95 ]fM(x)dx / max(∫[Q5,Q 95 ]fN(x)dx,1×10 -10 ) Among them, Q5 and Q 95 They represent the 5% and 95% quantiles of the distribution of malignant cells, respectively, defining the main distribution interval of malignant cells. The numerator and denominator represent the density integrals of malignant and non-malignant cells in this interval, respectively; (9) Peak separation calculation: Identify all local maxima of the density curve using kernel density estimation to obtain the array of peak positions peaks_positions, and calculate the minimum value of the distances between adjacent peaks: Peak_Separation = min(diff(peaks_positions)) / sd(CorrS_all) Where diff(peaks_positions) calculates the difference array of adjacent peak positions, diff(peaks_positions) = [peak2-peak1, peak3-peak2, ..., peak n -peak n -1], Peak_Separation measures the degree of separation between peaks in units of standard deviation; (10) Enhance the multimodality score: Adopt a conditional cumulative scoring mechanism for multi-dimensional evidence fusion: S_enhanced = min((∑score_conditions / 15) × 10, 10) Where 15 is the theoretical maximum total of each condition score, 10 is the upper limit of the final score, and the scoring rules for each condition are: BIC evidence: ΔBIC > 50 adds 3 points; 20 < ΔBIC ≤ 50 adds 2 points; 10 < ΔBIC ≤ 20 adds 1 point; Peak separation evidence: > 2σ adds 3 points; 1σ < separation ≤ 2σ adds 2 points; 0.5σ < separation ≤ 1σ adds 1 point, where σ is sd(CorrS_all); Relative density evidence: RDR > 10 adds 3 points; 5 < RDR ≤ 10 adds 2 points; 2 < RDR ≤ 5 adds 1 point; Local background evidence: LBR > 10 adds 3 points; 5 < LBR ≤ 10 adds 2 points; 2 < LBR ≤ 5 adds 1 point; Traditional multimodality evidence: Hartigan's Dip Test passes (p < 0.05) adds 3 points; When S_enhanced ≥ 5, it is determined as an enhanced multimodal distribution; The density peak analysis includes: (1) Calculate the kernel density estimate f(x) of the correlation distribution: Where K is the standard Gaussian kernel function K(u) = (1 / √2π)exp(-u 2 / 2), n is the total number of observations, h is the bandwidth parameter, the Silverman rule is used to determine the basic bandwidth, and then the adjustment factor adjust = 0.8 is set to improve the sensitivity to small peaks. The final bandwidth h_final = h × adjust, the number of sampling points is set to 1024 to ensure high-resolution density estimation, and xi is the i-th observation; (2) Implement a dynamic threshold adjustment strategy: relative_position = (x - x_min) / (x_max - x_min) position_gain = exp(2 × relative_position) / exp(2) adaptive_threshold = base_threshold / position_gain final_threshold = max(adaptive_threshold, min_threshold_limit) x_min and x_max are the minimum and maximum values of the correlation scores respectively, relative_position ∈ [0, 1] represents the relative position of the position within the distribution range, exp(2) is the normalization constant, base_threshold = 0.03 × max(f(x)) is the base height threshold, min_threshold_limit = 0.01 × max(f(x)) is the lower limit protection to prevent the threshold from being too low, and this strategy ensures that the high-correlation region obtains a lower detection threshold to improve sensitivity; (3) Implement additional threshold optimization for the high-correlation region: high_corr_bonus = 1 + (x - mean(CorrS_all)) / range(CorrS_all) optimized_threshold=adaptive_threshold / high_corr_bonus Where range(CorrS_all)=max(CorrS_all)-min(CorrS_all) is the range of the correlation score, high_corr_bonus≥1 is the gain coefficient for the high correlation region, and mean(CorrS_all) is the mean of the correlation scores of all cells. This optimization provides more sensitive detection capabilities for high correlation regions where malignant cells are prone to aggregate. (4) Detect the local maximum point p that meets the following conditions as the peak value: f(p)>f(p-1)ANDf(p)>f(p+1)ANDf(p)≥τfinal(p)ANDProm(p)≥min_prominence(p) Where min_prominence(p) is the minimum prominence threshold, the base value is set to max(f(x))×0.01, and is dynamically adjusted according to the position gain function; (5) Calculate the peak prominence Prom(p): Prom(p)=peak_height-max(min_left,min_right) Where peak_height is the density value at the peak, min_left and min_right are the lowest density values found on the left and right sides of the peak, respectively. Prominence reflects the significance of the peak relative to the surrounding valleys. Prom(p) > 0.01 × max(f(x)) indicates that the peak is sufficiently prominent. (6) Calculate the peak-to-peak distance Dist: Where p1 and p2 are the position coordinates of the two peaks, σ represents the standard deviation of the overall distribution, and Dist represents the distance between peaks in units of standard deviation, which is convenient for comparison between different samples; (7) Calculate the peak-to-valley ratio (PVR): Where valley_min represents the lowest density value between the two peaks, and PVR>1.3 indicates that the peak is well separated from the background; (8) Calculate the long-distance micro-peak reward DB: distance_score=|secondary_peak-main_peak| / sd(CorrS_all) height_ratio=secondary_peak_height / main_peak_height Among them, main_peak is the peak position with the highest density, secondary_peak is the secondary peak position, main_peak_height and secondary_peak_height are the corresponding density heights respectively, distance_score represents the normalized distance between peaks, and height_ratio represents the height ratio of the secondary peak to the main peak. The scoring rules are as follows: when distance_score > 2, if height_ratio > 0.5, it gets 4 points; if 0.3 ≤ height_ratio ≤ 0.5, it gets 3 points; if 0.1 ≤ height_ratio ≤ 0.3, it gets 2 points; if height_ratio ≤ 0.1, it gets 1 point; when 1 < distance_score ≤ 2, 1 point is deducted according to the above rules; (9) Peak reward in high correlation region: corr_position_bonus = min(2, (x - mean(CorrS_all)) / range(CorrS_all) × 5). Among them, x is the peak position, mean(CorrS_all) is the overall mean, and range(CorrS_all) is the value range. Peaks located in the high correlation region and on the right side of the main peak can obtain an additional bonus score of corr_position_bonus, with a maximum of 2 points; (10) Comprehensive score for bimodality: Bimodality_Score = min(∑feature_scores, 10) Among them, second_to_main_prominence_ratio = prominence_secondary / prominence_main represents the ratio of the prominence of the secondary peak to the main peak. feature_scores include the distance between peaks score, valley separation score, far - distance peak reward, high correlation region reward, secondary peak prominence score, multi - peak quantity reward, distribution asymmetry reward. When Bimodality_Score ≥ 4, it is determined that there is an obvious second peak; The comprehensive judgment adopts a hierarchical priority decision - making mechanism: (1) Highest - priority judgment: When Bimodality_Score ≥ 4, the judgment result is "There is an obvious second peak", indicating that the density peak analysis detects a significant second peak; (2) Second - highest - priority judgment: When S_enhanced ≥ 5 and at the same time satisfy NMD > 1.5, OVL < 0.3, ΔBIC > 10, BD > 0.5, the judgment result is "Distribution analysis confirms independence"; (3) Medium - priority judgment: When at the same time satisfy NMD > 1.0, OVL < 0.5, ΔBIC > 5, BD > Confidence_Score=min(10,max(0,0.5×(Bimodality_Score / 10)×5+0.2× (S_enhanced / 10)×2+0.3×min(ΔBIC / 10,3))) In the formula, peak analysis contributes 50% weight, enhanced multimodality analysis contributes 20% weight, BIC model comparison contributes 30% weight, and the final score range is [0,10]; (6) Verification report generation: The output includes a complete validation report containing the final judgment category, confidence score, key statistical indicators, and decision basis, providing a quantitative reliability assessment for downstream analysis.
9. The method for isolating malignant cells based on transcriptome sequencing according to claim 1, characterized in that: The step 7 comprises: (1) Generate a cell classification summary statistics table S, which includes the number and proportion of each type of cells; (2) Generate a CNV heatmap H, which shows the CNV patterns of different cell types, the behavior window is listed as cells, and the color indicates the CNV deviation intensity; (3) Generate a correlation distribution map D to show the distribution difference between malignant and non-malignant cells, including a nuclear density curve, a threshold line, and a peak marker; (4) Generate a visualization chart V of the verification results, including multimodality analysis and peak detection results; (5) The analysis results are saved as a structured data file R, which contains cell classification results, statistical indicators, and verification conclusions.
10. A malignant cell separation system based on transcriptome sequencing, characterized in that: The system is used to implement the method according to any one of claims 1 to 9, and the system comprises: Data input interface module, used to receive single-cell inferCNV prediction matrix, cell annotation information and gene location information; Preprocessing engine for performing gene sequencing and CNV deviation calculation; CNV calculation unit, used to implement sliding window calculation and CNV score quantification; A multi-level threshold processor, comprising a patient-level threshold calculation subunit and a sample-level threshold calculation subunit; a correlation analysis module, comprising a reference cell selector, a correlation calculator, and a threshold determiner; Cell classification decision unit, which performs dual-standard judgment based on CNV score and correlation; Independence verification system, including distribution analyzer and peak detector; Result output and visualization module, used to generate analysis reports and charts.