Construction method of tumor characteristic atlas based on ctc enrichment and multi-omics analysis

Abstract

The present application relates to the technical field of image processing, and further relates to a tumor feature atlas construction method based on CTC enrichment and multi-omics analysis. The method comprises the following steps: step one: sequentially performing spiral inertia focusing and deterministic lateral displacement sorting on a to-be-tested whole blood sample, and completing single-cell microcavity positioning capture of circulating tumor cells in the single-cell microcavity capture array region; step two: obtaining an oxygen molecule concentration value-time curve for each positioned and captured circulating tumor cell, and obtaining a metabolic function characteristic parameter group; step three: combining the metabolic function characteristic parameter group to form a multi-omics joint feature vector; and step four: constructing a circulating tumor cell subpopulation based on the multi-omics joint feature vector, extracting a marker molecular feature, and integrating and outputting a tumor feature atlas. The present application can provide more stable technical support for tumor heterogeneity analysis, metastasis potential evaluation, and efficacy monitoring.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, and particularly relates to image generation, specifically a method for constructing tumor feature maps based on CTC enrichment and multi-omics analysis. Background Technology

[0002] Circulating tumor cells (CTCs) serve as an important carrier for liquid biopsies, reflecting information such as tumor burden, metastatic risk, and drug resistance evolution without relying on tissue sampling. Therefore, they have attracted attention in tumor classification, staging, and efficacy assessment. Current technologies for CTC enrichment typically follow two approaches: one is immunoaffinity capture based on surface markers, such as magnetic bead capture or microstructure coating capture targeting epithelial-related markers. This approach offers high specificity, but CTCs may downregulate related markers during epithelial-mesenchymal transition, leading to missed detections and enrichment bias. The other approach is label-free sorting based on physical properties, such as microfiltration membrane pore size, inertial microfluidic size focusing, deterministic lateral displacement critical diameter separation, acoustic sorting, and dielectrophoretic sorting. This reduces marker dependence, but often results in insufficient enrichment purity or unstable recovery rates when faced with fluctuations in blood sample viscosity, differences in cell deformability, and overlapping leukocyte volume distribution. To improve the reliability of downstream analysis, some methods perform staining imaging or flow cytometry verification after enrichment. However, this can easily introduce fixation, dye loading, or phototoxicity that affects cell viability, thereby weakening the feasibility of subsequent functional and omics assays.

[0003] At the single-cell level, transcriptome sequencing is relatively mature, and single-cell proteomics and metabolomics detection are rapidly developing. Multi-omics data fusion often employs typical alignment and dimensionality reduction frameworks, such as alignment based on correlation analysis, batch correction using nearest neighbor matching, and joint representation learning of latent variable models, to map different modalities to a shared space and perform clustering or trajectory analysis. However, existing multi-omics methods still face two key contradictions in practical applications: First, multi-omics data often come from different cells or different batches of samples, and the interplay of intercellular heterogeneity and batch effects leads to uncertainty in the subpopulation boundaries and biological interpretations obtained after fusion; second, even when multi-omics data are obtained from the same batch, real-time characterization of cell functional status is often lacking, especially metabolic phenotypes closely related to tumor invasion and drug resistance, which are difficult to stably quantify at the single-cell scale. In existing metabolic phenotyping assays, external throughput analyzers can output oxygen consumption rate and glycolysis-related indicators, but they are usually measured on a per-well plate or per-cell population basis. Single-cell resolution is insufficient and it is difficult to achieve a one-to-one correspondence with single-cell transcriptomics, proteomics, and metabolomics. Optical methods based on fluorescent or phosphorescent oxygen probes can achieve higher spatial resolution, but within a small volume, factors such as background light, optical path drift, probe membrane thickness inhomogeneity, photobleaching, and environmental oxygen diffusion can significantly affect signal stability, easily causing reading fluctuations. As a result, the oxygen molecule concentration change curve over time is difficult to use for reliable rate estimation. Summary of the Invention

[0004] The main objective of this invention is to provide a method for constructing tumor feature maps based on CTC enrichment and multi-omics analysis, which reduces enrichment bias caused by dependence on immune markers and reduces cell damage and identity mismatch caused by multiple centrifugation and transport, thereby improving the reliability and reproducibility of subsequent single-cell analysis from the source.

[0005] To address the aforementioned problems, the technical solution of this invention is as follows: a method for constructing a tumor feature atlas based on CTC enrichment and multi-omics analysis, the method comprising: Step 1: The whole blood sample to be tested is sequentially subjected to helical inertial focusing and deterministic lateral displacement sorting, and the single-cell microcavity localization and capture of circulating tumor cells is completed in the single-cell microcavity capture array area; Step 2: Obtain the oxygen molecule concentration value over time for each captured circulating tumor cell and obtain a set of metabolic function characteristic parameters; Step 3: Release circulating tumor cells and obtain gene expression profile data, protein expression profile data and metabolite profile data respectively, and combine them with metabolic function characteristic parameter groups to form a multi-omics joint feature vector; Step 4: Construct circulating tumor cell subpopulations based on multi-omics joint feature vectors, infer the subpopulation evolution directed graph, extract landmark molecular features, and integrate and output tumor feature atlas.

[0006] Furthermore, the spiral inertial focusing in step one includes: injecting the whole blood sample to be tested into the spiral inertial focusing channel. The cross-section of the spiral inertial focusing channel is rectangular and the width remains constant along the spiral direction. When the whole blood sample to be tested flows in the spiral inertial focusing channel, it is simultaneously subjected to inertial lift and Dean drag force. Cells with a diameter greater than a preset size threshold form a focusing band at the equilibrium position of the two forces and move along the inner wall of the channel. Red blood cells and platelets migrate to the outside of the channel with the secondary Dean flow and are discharged through the waste liquid outlet, thereby obtaining a primary enriched suspension containing white blood cells and circulating tumor cells.

[0007] Furthermore, the deterministic lateral displacement sorting in step one includes: introducing the primary enriched suspension into the deterministic lateral displacement sorting area, setting up a periodically arranged array of micropillars in the deterministic lateral displacement sorting area, wherein adjacent columns of micropillars in the micropillar array have a fixed lateral offset along the flow direction, when the cell diameter is greater than the critical separation diameter, the cell undergoes lateral displacement along the micropillar offset direction after contacting the micropillar and enters the circulating tumor cell collection channel, when the cell diameter is less than the critical separation diameter, the cell passes through the micropillar gap along the streamline direction and enters the leukocyte discharge channel, thereby obtaining the secondary enriched suspension of circulating tumor cells.

[0008] Furthermore, the single-cell microcavity localization and capture in step one includes: guiding a secondary enrichment suspension of circulating tumor cells into the single-cell microcavity capture array region; the bottom surface of the single-cell microcavity capture array region is provided with a matrix of cylindrical microcavities, the diameter of each cylindrical microcavity being used to accommodate a single circulating tumor cell; a micropore is opened at the center of the bottom of each cylindrical microcavity and the micropore is connected to a negative pressure chamber, so that the circulating tumor cell is drawn into the cylindrical microcavity under negative pressure and the micropore is blocked; after blocking, the negative pressure in the cylindrical microcavity disappears, thereby forming a single-cell occupancy; an oxygen-sensitive phosphorescent sensing membrane is pre-deposited on the inner wall surface of each cylindrical microcavity, the oxygen-sensitive phosphorescent sensing membrane being formed by uniformly dispersing platinum porphyrin luminescent molecules in a polystyrene matrix.

[0009] Further, step two includes: using a direct digital frequency synthesis chip to generate a sinusoidal modulated voltage signal with a fixed frequency; the sinusoidal modulated voltage signal drives an ultraviolet light-emitting diode to generate modulated excitation light through a voltage-to-current conversion circuit; the modulated excitation light is collimated by a collimating lens and irradiates the single-cell microcavity trapping array region to excite the oxygen-sensitive phosphorescence sensing membrane to generate red-band phosphorescence emission; the phosphorescence emission is collected by an objective lens and filtered by a red bandpass filter to remove excitation light scattering before being input into a photomultiplier tube; the photocurrent signal output by the photomultiplier tube is converted into a voltage signal by a transimpedance amplifier and digitally sampled by an analog-to-digital converter; the digitally sampled phosphorescence signal is subjected to phase-locked demodulation processing in a field-programmable gate array; the phase-locked demodulation processing includes generating a reference sine signal and a cosine signal obtained by shifting the reference sine signal by 90 degrees; the digitally sampled phosphorescence signal is multiplied point-by-point with the reference sine signal and the cosine signal respectively, and then passed through a finite impulse response low-pass filter to obtain a DC in-phase component signal and a DC quadrature component signal respectively.

[0010] Furthermore, step two also includes: determining the phase delay of phosphorescence emission relative to the modulated excitation light based on the DC in-phase component signal and the DC quadrature component signal, and converting the phase delay into the oxygen molecule concentration value in the cylindrical microcavity based on the pre-calibrated quenching correspondence; continuously collecting the oxygen molecule concentration value in each cylindrical microcavity and recording the change curve of the oxygen molecule concentration value over time, and determining the decrease in the oxygen molecule concentration value per unit time as the basal oxygen consumption rate of the corresponding circulating tumor cells.

[0011] Furthermore, step two also includes: after the basal oxygen consumption rate is measured, metabolic regulation reagents are sequentially perfused into the single-cell microcavity capture array region through the reagent injection channel of the microfluidic chip. The metabolic regulation reagents include ATP synthase inhibitor solution, mitochondrial uncoupling agent solution, and electron transport chain complex inhibitor solution. The basal oxygen consumption rate before perfusion of ATP synthase inhibitor solution minus the oxygen consumption rate after perfusion of ATP synthase inhibitor solution is taken as ATP synthesis-related oxygen consumption. The peak oxygen consumption rate after perfusion of mitochondrial uncoupling agent solution is taken as maximum respiratory capacity. The residual oxygen consumption after perfusion of electron transport chain complex inhibitor solution is taken as non-mitochondrial oxygen consumption, thereby forming a set of metabolic function characteristic parameters, which includes basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption.

[0012] Furthermore, step three includes: introducing positive pressure gas into the negative pressure chamber to release the negative pressure adsorption at the micropores and perfusing the single-cell microcavity capture array area with cell release buffer, so that circulating tumor cells are released from the cylindrical microcavities and flow out to the single-cell sorting and collection area; in the single-cell sorting and collection area, circulating tumor cells are individually aspirated by micromanipulation capillaries and transferred to independent reaction tubes; single-cell transcriptome sequencing is performed on each circulating tumor cell in each independent reaction tube to obtain gene expression profile data, single-cell proteome mass spectrometry is performed to obtain protein expression profile data, and single-cell metabolome mass spectrometry is performed to obtain metabolite profile data.

[0013] Furthermore, step three also includes: calculating the median basal oxygen consumption rate of all circulating tumor cells as the basal oxygen consumption rate cutoff value; calculating the median maximum respiratory capacity of all circulating tumor cells as the maximum respiratory capacity cutoff value; calculating the difference between the maximum respiratory capacity and the basal oxygen consumption rate of all circulating tumor cells and taking the median as the respiratory reserve cutoff value; and classifying all circulating tumor cells into metabolic phenotypes based on metabolic function characteristic parameter groups. Circulating tumor cells with a basal oxygen consumption rate greater than the basal oxygen consumption rate cutoff value and a maximum respiratory capacity greater than the maximum respiratory capacity cutoff value are classified into oxidized phosphoric acid cells. Circulating tumor cells (CTCs) with a basal oxygen consumption rate (BAC) less than the BAC threshold and a quotient of ATP synthesis-related oxygen consumption and BAC less than a preset threshold are classified into the glycolysis-dependent layer. CTCs with a maximum respiratory capacity greater than the difference between their BAC and respiratory reserve are classified into the metabolic reserve-sufficient layer. CTCs belong to multiple metabolic phenotype layers. Within each metabolic phenotype layer, gene expression profiles, protein expression profiles, and metabolite profiles are standardized and concatenated in gene-protein-metabolite order to form a multi-omics joint feature vector.

[0014] The tumor feature map construction method based on CTC enrichment and multi-omics analysis of the present invention has the following beneficial effects: In terms of metabolic phenotype acquisition, the present invention adopts a phase-sensitive detection approach based on lock-in amplification, which transforms phosphorescence lifetime information into robust readouts of oxygen molecule concentration changes. It can maintain a high signal-to-noise ratio under complex background light and optical path perturbations, making the oxygen molecule concentration value change curve over time more suitable for rate estimation and drug perturbation response analysis. This forms a set of metabolic function feature parameters corresponding to a single cell, and enables this set of parameters to be accurately aligned with gene expression profile data, protein expression profile data, and metabolite profile data at the cell numbering level, constructing a multi-omics joint feature vector, which significantly reduces batch drift and cross-cell splicing errors commonly found in multimodal data fusion.

[0015] In the map construction stage, this invention introduces metabolic continuity-constrained community detection into the dimensionality-reduced neighborhood structure, ensuring that the classification of circulating tumor cell subpopulations simultaneously satisfies molecular similarity and functional compatibility. This reduces over-clustering or erroneous merging caused solely by noise or sparsity. Furthermore, an evolutionary directed graph is established at the subpopulation level based on differences in metabolic activity, providing more interpretable directional information for typing, staging, and evolutionary mechanism studies. Simultaneously, by comparing molecular indicators within and outside subpopulations to extract sets of positive and negative biomarkers, the map output not only includes structured subpopulation relationships but also provides candidate molecular clues for target screening and validation. Overall, this invention integrates enrichment, capture, metabolic dynamic monitoring, and multi-omics acquisition into a closed-loop process, balancing high throughput and single-cell accuracy. This improves the completeness, relevance, and reproducibility of circulating tumor cell phenotypic characterization, providing more stable technical support for tumor heterogeneity analysis, metastatic potential assessment, and efficacy monitoring. Attached Figure Description

[0016] Figure 1 A schematic diagram illustrating the core principle of the circulating tumor cell subpopulation boundary identification method based on nearest-neighbor undirected graphs provided in this embodiment of the invention; Figure 2 This is a standardized numerical distribution heatmap of metabolic functional characteristic parameters of various circulating tumor cell subsets provided in the embodiments of the present invention. Detailed Implementation

[0017] A method for constructing tumor feature atlases based on CTC enrichment and multi-omics analysis, including: Step 1: The whole blood sample to be tested is sequentially subjected to helical inertial focusing and deterministic lateral displacement sorting, and the single-cell microcavity localization and capture of circulating tumor cells is completed in the single-cell microcavity capture array area.

[0018] In a feasible implementation, the whole blood sample to be tested is first treated with anticoagulation at the collection end to maintain stable cell morphology and activity. For example, after collection using dipotassium ethylenediaminetetraacetate anticoagulation tubes, the sample is processed within 2 hours. To reduce the impact of blood viscosity fluctuations on the flow pattern within the helical inertial focusing channel, the whole blood sample to be tested can be mixed with an equal volume of phosphate-buffered saline to form the working solution. Before injection, the microfluidic chip is pre-filled and defoamed: phosphate-buffered saline is used to flush from the inlet end at a flow rate of 0.5 mL / min for 2 minutes, followed by flushing with a buffer solution containing 0.1% polyethylene glycol surfactant for 1 minute. This forms a hydrophilic thin layer on the channel wall, reducing the probability of non-specific cell adsorption in the helical inertial focusing channel and the deterministic lateral displacement sorting region, while preventing air bubbles from entering the single-cell microcavity capture array region and causing local flow field abrupt changes. Subsequently, the whole blood sample to be tested is injected into the inlet of the microfluidic chip at a constant flow rate. The constant flow rate can be selected from any value between 1.2 ml / min and 2.0 ml / min. In this example, 1.5 ml / min is used, so that the whole blood sample to be tested passes through the spiral inertial focusing channel, the deterministic lateral displacement sorting area and enters the single-cell microcavity capture array area in sequence.

[0019] The helical inertial focusing channel is used to push larger-diameter cells to a stable lateral equilibrium position under continuous flow conditions, thereby separating red blood cells and platelets from the major cell band. The cross-section of the helical inertial focusing channel is rectangular, with a constant width along the helical direction. Example dimensions include a width of 300 micrometers and a height of 80 micrometers, with 6 helical turns and a helical radius gradually transitioning from 3 millimeters to 12 millimeters. This geometry allows the whole blood sample to form a repeatable secondary flow structure within the helical inertial focusing channel, where cells exhibit significant size-dependent migration under the combined effects of inertial lift and Dean's drag force. To make "stable focusing" a controllable goal, a common practice is to confine the flow state to a range where the inertial effect is significant but does not cause severe shear damage. This can be characterized by the Reynolds number, written as Reynolds number. ,in This indicates the working solution density of the whole blood sample to be tested. This represents the average flow velocity across the cross-section within the spiral inertial focusing channel. This indicates the hydraulic diameter of the helical inertial focusing channel. This represents the dynamic viscosity of the working solution in the whole blood sample being tested. In the example, we take... 1050 kg per cubic meter 0.002 Pa second Approximately 133 micrometers At approximately 0.8 meters per second, At the tens of magnitude, the inertial lift is sufficient to "lift" cells with diameters larger than a preset size threshold from the streamlines and push them to a lateral equilibrium position. Meanwhile, the Dean's secondary flow introduced by the helical curvature can be characterized by the Dean number, written as... ,in This represents the local radius of curvature of the helical inertial focusing channel. With... Gradually increases along the spiral. The spatial distribution tends to be smooth, and the intensity of the secondary flow will not suddenly increase in a certain section, thereby reducing the probability of cells being subjected to abnormal impact in a certain local area. Based on this mechanical pattern, cells with a diameter larger than the preset size threshold are more likely to form a focusing zone at the equilibrium position of the two forces and move along the inner wall of the channel, while red blood cells and platelets, due to their smaller size and weaker inertial lift, are more likely to migrate to the outer side of the channel with the Dean secondary flow and be discharged through the waste liquid outlet. The preset size threshold can be selected from any value between 10 micrometers and 12 micrometers, so as to retain most white blood cells and circulating tumor cells in the primary enrichment suspension, while preferentially removing red blood cells and platelets at the fluid level. In order to improve the repeatability of the primary enrichment suspension, a shunt geometry can also be arranged at the outlet end of the spiral inertial focusing channel, so that the focusing zone near the inner wall enters the main outlet, while the cells and blood components near the outer wall enter the waste liquid outlet; the shunt ratio can be set to the main outlet accounting for 20% to 40% of the total flow, and 30% is used in the example. In this way, the primary enrichment suspension obtained by the main outlet is more concentrated, and the burden on the subsequent deterministic lateral displacement sorting area is lighter.

[0020] After the primary enrichment suspension enters the deterministic lateral displacement sorting region, a micropillar array is used to perform more refined size threshold screening of cells, thereby increasing the relative enrichment of circulating tumor cells in the secondary enrichment suspension from the leukocyte background. A periodically arranged micropillar array is set up within the deterministic lateral displacement sorting region. In this example, the micropillar diameter is 25 micrometers, the micropillar spacing is 18 micrometers, adjacent columns of micropillars have a fixed lateral offset of 2 micrometers along the flow direction, the array column spacing is 45 micrometers, and the array row spacing is 30 micrometers. This array geometry corresponds to a critical separation diameter range, which can be obtained through calibration in engineering: channel tests are conducted using polystyrene microspheres of known diameters of 8 micrometers, 10 micrometers, and 12 micrometers, observing the final outlet distribution and determining that the critical separation diameter falls within approximately 10 micrometers. The advantage of deterministic lateral displacement sorting is that whether a cell undergoes lateral displacement is not a "probabilistic deflection," but a "cumulative offset" determined by the geometry of the micropillar array. When the cell diameter is larger than the critical separation diameter, the repulsive displacement of the cell after each contact with the micropillar is continuously accumulated by the fixed lateral offset of the array, eventually leading to entry into the circulating tumor cell collection channel. When the cell diameter is smaller than the critical separation diameter, the cell is more likely to pass through the micropillar gaps along the streamline direction and enter the leukocyte discharge channel. To reduce the effective diameter drift caused by cell deformation, the channel height of the deterministic lateral displacement sorting area can be set to be consistent with the micropillar height and slightly larger than the diameter of a typical circulating tumor cell. In the example, the channel height is 35 micrometers, which restricts the cells in an approximately two-dimensional manner when passing through the array, thereby improving the separation consistency corresponding to the critical separation diameter. Optionally, to prevent large monocytes from the leukocyte subset from entering the circulating tumor cell collection channel, a straight rectifier channel can be arranged before and after the deterministic lateral displacement sorting area. A slight shearing pre-shaping section can be added within the rectifier channel, with the shear rate controlled at 3000 to 5000 per second. This ensures that the transient deformation of deformable cells is more uniform, reducing fluctuations in enrichment purity caused by "accidental crossings." After passing through the aforementioned deterministic lateral displacement sorting area, a secondary enriched suspension of circulating tumor cells is obtained and introduced into the single-cell microcavity capture array area.

[0021] The goal of the single-cell microcavity capture array region is to capture circulating tumor cells (CTCs) in a secondary enrichment suspension at the single-cell level and stably confine each captured CTC within a corresponding cylindrical microcavity, facilitating subsequent phase-sensitive single-cell oxygen metabolism dynamic monitoring. The bottom surface of the single-cell microcavity capture array region features a matrix arrangement of cylindrical microcavities. In this example, the diameter of each cylindrical microcavity is 18 micrometers, the depth is 25 micrometers, and the array size is 100 rows by 100 columns to provide 10,000 capture sites. A micropore is located at the center of the bottom of each cylindrical microcavity and is connected to a negative pressure chamber. The micropore diameter is 4 micrometers in this example, and the micropore length is 20 micrometers. At the start of capture, a negative pressure chamber is connected to a negative pressure source to establish a stable negative pressure, which can be selected from 8 kPa to 15 kPa; in this example, 12 kPa is used. This allows circulating tumor cells to flow through the microcavity along with the secondary enrichment suspension, causing local fluid to be attracted by the micropores, forming downward streamlines. Cells, under the combined action of streamline traction and the geometric constraint of the microcavity inlet, enter the cylindrical microcavity and block the micropores. After blockage, the fluid pathway within the cylindrical microcavity is blocked, the volumetric flow rate at the micropores drops sharply, and the pressure difference between the inside and outside of the cylindrical microcavity decreases rapidly, resulting in the disappearance of negative pressure within the cylindrical microcavity. This prevents subsequent cells from entering the same cylindrical microcavity, thus preventing single-cell occupancy. This "blockage-to-self-termination" mechanism allows the capture process to operate without complex valve switching and avoids individual control of each capture site under high-throughput conditions. The more capture sites, the greater the benefit. To reduce the probability of micropores being blocked prematurely by cell debris or protein clots, a 30-micron coarse filter grid can be arranged before the entrance of the single-cell microcavity capture array area. At the same time, the working solution of the whole blood sample to be tested is pre-filtered through a cell filter with a 40-micron pore size before entering the microfluidic chip, ensuring that the main flow enters the array with "cells" as the main particles.

[0022] An oxygen-sensitive phosphorescent sensing film is pre-deposited on the inner wall of a cylindrical microcavity. This film is formed by uniformly dispersing platinum porphyrin luminescent molecules within a polystyrene matrix. A feasible preparation method involves dissolving polystyrene in dichloromethane to obtain a 3% (w / w) polystyrene solution, then adding platinum porphyrin luminescent molecules to a concentration of 0.2 mg / mL. After thorough mixing by shaking, the solution is spin-coated onto the surface of the single-cell microcavity capture array region. The spin-coating speed is 1500 rpm for 30 seconds, followed by standing at 25°C for 10 minutes to allow solvent evaporation and form a uniform film layer. To confine the sensing film to the inner wall of the cylindrical microcavity and avoid covering the micropores, a short-term hydrophobic shielding treatment can be applied to the micropore locations before spin-coating. For example, a small amount of fluorinated silane reagent can be drop-coated and baked at 60°C for 5 minutes, creating a low-wetting boundary at the micropores. This causes the spin-coated liquid to adhere more readily to the non-micropore areas on the sidewalls and bottom of the cylindrical microcavity. After deposition, the unattached residue was removed by rinsing with phosphate-buffered saline at a rate of 0.5 mL / min for 1 minute, and then the microfluidic chip was sealed for later use. The advantage of this process is that the oxygen-sensitive phosphorescence sensing membrane is located inside the cylindrical microcavity. Subsequent cells entering the cylindrical microcavity are in the same microenvironment as the oxygen-sensitive phosphorescence sensing membrane. This results in a shorter coupling path for changes in oxygen concentration to the phosphorescence response, leading to a faster response. Simultaneously, perturbations in the oxygen concentration of the external main channel have a smaller impact on the readings within the single-cell microcavity.

[0023] To achieve "localization" in single-cell microcavity localization and capture, a occupancy identification process is typically performed after capture, and a cylindrical microcavity index is established. An example approach is to scan the single-cell microcavity capture array area using bright-field imaging, acquiring images covering the entire array, and then numbering the cylindrical microcavities according to the array rows and columns, for example, starting from the top left corner and prioritizing row-based numbering. Occupancy identification can be based on grayscale difference and circular contour matching: empty cylindrical microcavities exhibit a stable bright ring structure at the cavity opening, while cylindrical microcavities occupied by cells show a distinct central absorbent area at the cavity opening. This difference can be used to separate the "set of occupied cylindrical microcavities" from the "set of unoccupied cylindrical microcavities" for recording. The direct benefit of this approach is that subsequent dynamic monitoring of oxygen metabolism in each circulating tumor cell can be performed by target tracking based on the "cylindrical microcavity number," avoiding misreading or duplicate readings in high-density arrays. It also facilitates consistent correlation between "capture time, capture location, and pre-capture flow conditions" and subsequent single-cell data. Optionally, occupancy identification can also be combined with the flow indication signal of the micropore at the bottom of the microcavity for cross-validation: a micro flow sensor is connected in series at the inlet of the negative pressure chamber. At the beginning of the capture, the total flow rate is relatively large. As the cylindrical microcavity is gradually blocked, the total flow rate decreases in a step-like manner. When the total flow rate drops to 10% to 20% of the initial flow rate, the array occupancy rate usually reaches a high level. At this time, the injection of the secondary enrichment suspension of circulating tumor cells is stopped and the buffer solution is switched to maintain the flow, which can reduce the subsequent accumulation of cells at the array inlet.

[0024] In terms of alternative implementation methods, if the viscosity of the whole blood sample to be tested is high or the proportion of red blood cells is high, the focusing band of the helical inertial focusing channel may drift slightly. In this case, the constant flow rate can be reduced from 1.5 ml / min to 1.2 ml / min, and the upper limit of the helical radius of the helical inertial focusing channel can be adjusted from 12 mm to 15 mm to expand the "focusing stability range" by reducing the Dean drag force. Correspondingly, in order to keep the overall processing time from increasing significantly, the array size of the single-cell microcavity capture array region can be increased from 100 rows by 100 columns to 120 rows by 120 columns, so that the number of capture sites increases with the throughput demand. If the sample contains a large number of leukocyte subsets with diameters similar to circulating tumor cells, the critical separation diameter of the deterministic lateral displacement sorting region can be adjusted to approximately 11 micrometers to achieve a stricter size threshold. Specifically, this involves reducing the microcolumn gap to 16 micrometers and decreasing the fixed lateral offset between adjacent microcolumns to 1.5 micrometers. This makes the "cumulative lateral displacement condition" required for cells to enter the circulating tumor cell collection channel more stringent, thereby improving the enrichment consistency of the secondary enrichment suspension of circulating tumor cells. If greater emphasis is placed on post-capture cell viability, the negative pressure in the negative pressure chamber can be reduced from 12 kPa to 8 kPa, and the capture time extended, for example, from 3 minutes to 6 minutes. This gentler aspiration process reduces the peak cell membrane stress, resulting in a decrease in the occupancy velocity.

[0025] Through the above process, the whole blood sample to be tested completes helical inertial focusing and deterministic lateral displacement sorting, and completes single-cell microcavity localization and capture of circulating tumor cells in the single-cell microcavity capture array area, providing a stable, numberable, and traceable single-cell carrier environment for subsequent phase-sensitive single-cell oxygen metabolism dynamic monitoring based on lock-in amplification.

[0026] Step 2: Obtain the oxygen molecule concentration value over time for each captured circulating tumor cell and obtain a set of metabolic function characteristic parameters; Step 3: Release circulating tumor cells and obtain gene expression profile data, protein expression profile data and metabolite profile data respectively, and combine them with metabolic function characteristic parameter groups to form a multi-omics joint feature vector.

[0027] After capturing circulating tumor cells in the single-cell microcavity capture array region, the external environment of the microfluidic chip is stabilized to repeatable conditions before dynamic monitoring of oxygen metabolism begins. A common practice is to fix the microfluidic chip on a temperature-controlled platform of the microscope stage, maintaining the platform temperature at 37 degrees Celsius, and continuously introducing culture buffer containing a bicarbonate buffer system at the inlet of the single-cell microcavity capture array region to stabilize the solution environment within the cylindrical microcavity within 10 minutes. The direct benefit of this is that the temperature-related drift in the oxygen concentration-time curve is compressed to a smaller range, and the cellular respiration rate is closer to steady state within the observation window, facilitating subsequent estimation of the basal oxygen consumption rate using the slope. To avoid background rise caused by oxygen molecules from the outside air permeating along the polymer material, a barrier membrane can be coated around the culture buffer flow path, or the microfluidic chip can be placed in a sealed cavity and the cavity atmosphere can be maintained using a gas mixture with a fixed oxygen volume fraction, such as maintaining calibration and measurement under air-level conditions, ensuring consistent boundary conditions between different batches.

[0028] Phase-locked amplification (LLDPA) phase-sensitive single-cell oxygen metabolism dynamic monitoring begins with modulation excitation. A direct digital frequency synthesis chip outputs a sinusoidal modulated voltage signal with a fixed frequency; an example frequency of 5000 Hz can be used. The choice of 5000 Hz is based on the fact that the phosphorescence lifetime of the oxygen-sensitive phosphorescent sensing film is typically on the order of microseconds. If the modulation frequency is too low, the ability to resolve lifetime changes due to phase delay decreases; if the modulation frequency is too high, the photoelectric detection link becomes more sensitive to high-frequency noise and experiences more significant signal amplitude attenuation. 5000 Hz represents a trade-off between phase resolution and electronic bandwidth constraints. The sinusoidal modulated voltage signal drives an ultraviolet (UV) light-emitting diode (LED) via a voltage-to-current converter. The center wavelength of the UV LED can be 405 nm, and the output power can be any value between 2 mW and 8 mW; an example of 5 mW is used. The modulated excitation light, collimated by a collimating lens, vertically illuminates the single-cell microcavity trapping array region, causing the oxygen-sensitive phosphorescent sensing film on the inner wall of the cylindrical microcavity to emit red-band phosphorescence. The center wavelength of the red-band can be 650 nm. Phosphorescence emission intensity varies with oxygen concentration, but more importantly, the phosphorescence emission phase exhibits a phase delay relative to the modulated excitation light. This phase delay is determined by the lifetime of the excited state of the luminescent molecules, which is in turn affected by oxygen quenching. Therefore, the phase delay provides a more stable readout path for oxygen concentration. Using phase readout instead of relying solely on intensity readout can reduce systematic errors caused by fluctuations in optical transmittance, uneven film thickness, and differences in illumination between microcavities. These factors are more likely to affect the intensity amplitude, while the phase delay is mainly determined by excited state dynamics and is insensitive to slowly changing amplitude perturbations.

[0029] Phosphorescent emission is collected by the objective lens and then filtered by a red bandpass filter to remove excitation light scattering. The center wavelength of the bandpass filter can be 650 nm, and the bandwidth can be 40 nm. The filtered phosphorescent signal enters a photomultiplier tube (PMT), with a gain on the order of 1,000,000. The output photocurrent signal is converted into a voltage signal by a transimpedance amplifier and then input to an analog-to-digital converter (ADC) for digital sampling. The ADC sampling rate needs to cover the modulation frequency and its main noise components; an example sampling rate of 200,000 times per second and a quantization bit depth of 16 bits are used to ensure sufficient numerical resolution for the phase estimation after phase-locked loop demodulation. To reduce interference from power frequency and ambient light, a light shield can be used at the PMT input, and an analog high-pass filter can be added at the front end of the transimpedance amplifier to suppress DC drift. The analog high-pass cutoff frequency can be any value between 50 Hz and 200 Hz; an example of 100 Hz is used.

[0030] The digitized phosphorescent signal is subjected to phase-locked demodulation (PLM) within a field-programmable gate array (FPGA). The core of PLM is shifting the signal component at the target frequency to DC and then suppressing unrelated frequency components using a low-pass filter. Specifically, the digitized phosphorescent signal is split into two paths: one path is multiplied point-by-point with a reference sine signal, and the other path is multiplied point-by-point with a cosine signal obtained by shifting the reference sine signal by 90 degrees. These two multiplication operations yield in-phase and quadrature mixing sequences, respectively. Subsequently, the in-phase and quadrature mixing sequences are passed through finite impulse response (FIR) low-pass filters. The passband of the low-pass filter can be between 0 Hz and 20 Hz, and the filter order can be between 512 and 2048; in this example, 1024 is used, thus outputting DC in-phase and DC quadrature component signals. The reason for setting the passband of the low-pass filter to within 20 Hz is that the change in oxygen molecule concentration caused by single-cell respiration is a slow variable, with the change characteristics usually on the order of seconds to tens of seconds. Phase-locked demodulation does not need to retain dynamic components higher than tens of Hz. Narrowing the passband can significantly improve the signal-to-noise ratio and reduce the phase noise caused by the small jitter of the excitation light power.

[0031] After obtaining the DC in-phase component signal and the DC quadrature component signal, the phase delay of phosphorescent emission relative to the modulated excitation light is calculated. The phase delay can be calculated as follows: Calculation, where Indicates phase delay, This represents the amplitude of the DC in-phase component signal. This represents the amplitude of the DC quadrature component signal. This represents the arctangent function. The advantage of using the arctangent is that... and Simultaneously, when affected by the overall scaling of optical power, its ratio changes little, and the phase delay estimate remains stable. To avoid The phase transition occurs when the phase approaches zero, which can be implemented using the four-quadrant arctangent operation and... , For saturation protection, use the arctangent operation in the four quadrants. and The sign of the phase determines the quadrant in which the phase is located, thus keeping the phase varying within a continuous interval.

[0032] The relationship between phase delay and excited state lifetime can be used It means that among them Indicates the excited state lifetime. Represents the tangent function. Indicates phase delay, Indicates the frequency of the modulated excitation light. This represents pi. The significance of writing it in this form is that, with a fixed modulation frequency, the phase delay monotonically maps to the excited-state lifetime, and lifetime changes can be directly derived from phase changes. The correspondence between the excited-state lifetime and the oxygen molecule concentration in the cylindrical microcavity can be established through pre-calibration. During calibration, a cell-free calibration buffer is injected into the single-cell microcavity capture array region, and gas mixtures with different oxygen molecule volume fractions are applied to the external atmosphere-controlled cavity. For example, five points are selected: 0%, 5%, 10%, 15%, and 21%. Data is collected after each point has stabilized for 5 minutes. and And calculate and This yields a curve showing the relationship between lifetime and oxygen concentration. To convert volume fraction to oxygen concentration in solution, Henry's Law can be used for a consistent conversion, or dissolved oxygen can be simultaneously measured in the calibration buffer using a commercially available dissolved oxygen electrode, and a lookup table mapping can be generated. If a quenching model is used, it can be expressed as follows: ,in This represents the excited-state lifetime under anaerobic conditions. Indicates the current excited state lifetime. Denotes the quenching constant. This represents the oxygen concentration within the cylindrical microcavity. By fixing the quenching constant and the film material under the same process conditions, the response differences between different cylindrical microcavities can be compressed to a smaller range. Finally, during operation, the phase delay is converted into an oxygen concentration value through a lookup table. To balance speed and stability, a combination of lookup table and linear interpolation is commonly used during operation: the calibrated discrete points are stored in the table, and when the current phase delay falls between two adjacent points, the oxygen concentration value is calculated proportionally and output.

[0033] After completing the conversion from phase delay to oxygen concentration value, oxygen concentration values ​​are continuously acquired for each captured circulating tumor cell, and the curve of oxygen concentration change over time is recorded. The acquisition interval can be between 0.5 seconds and 2 seconds, with 1 second used in this example. Within each acquisition point, lock-in demodulation accumulation is still used to improve the signal-to-noise ratio. For example, setting the lock-in demodulation accumulation window to 100 modulation cycles results in an accumulation time of 0.02 seconds for a 5000 Hz modulated signal. This ensures that each output point contains sufficient averaging to suppress noise while maintaining the ability to track second-level metabolic changes. To reduce the impact of fluctuations in the main oxygen concentration outside the microcavity on the microcavity readings, a constant flow rate of culture buffer is maintained during acquisition. The flow rate can be between 5 μL / min and 30 μL / min, with 10 μL / min used in this example. This allows for slow fluid turnover above the cylindrical microcavity without causing cell dislocation, while also ensuring more consistent external oxygen supply conditions across the entire array.

[0034] The baseline oxygen consumption rate is extracted from the curve of oxygen molecule concentration versus time. One feasible method is to perform a linear regression on a stability window, with the regression slope representing the oxygen consumption rate. This represents the basal oxygen consumption rate, where Indicates the basal oxygen consumption rate. This indicates the change in oxygen molecule concentration. This indicates the change over time; the negative sign is used to represent a decrease in oxygen concentration as a positive rate of oxygen consumption. To allow... and The values ​​are unaffected by the initial transient. After capture, a 60-second wait can be initiated, followed by a subsequent 180-second regression window. This approach is advantageous because the capture instant may introduce fluid disturbance within the microcavity, causing the oxygen concentration to initially exhibit a non-linear adjustment. Avoiding this interval allows the basal oxygen consumption rate to more closely approximate the steady-state respiration level of circulating tumor cells. Optionally, to address the slow recovery of the oxygen concentration curve due to minor leakage in some microcavities, the regression model can be extended to a robust regression with intercept and drift terms, or the curve can be first median-filtered before regression. The filtering window can be set to 5 points to suppress impulse noise.

[0035] After the basal oxygen consumption rate was measured, metabolic regulatory reagents were sequentially perfused into the single-cell microcavity capture array region through the reagent injection channel of the microfluidic chip to obtain a set of metabolic function characteristic parameters. The reagent perfusion followed an "inject-dwell-rinse" rhythm to ensure the drug effect within the microcavity reached a quasi-steady state before reading the oxygen consumption rate. First, an ATP synthase inhibitor solution was perfused; oligomycin solution, at a concentration of 2 μmol / L, was used. The injection time was 60 seconds, the dwell time was 180 seconds, and then the solution was rinsed with culture buffer for 120 seconds. After oligomycin inhibited ATP synthase, mitochondrial proton reflux was obstructed, electron transport flux decreased, and the oxygen consumption rate decreased accordingly. Therefore, the basal oxygen consumption rate before perfusion of the ATP synthase inhibitor solution was subtracted from the steady-state oxygen consumption rate after perfusion to obtain the ATP synthesis-related oxygen consumption. To avoid mistaking short-term flow disturbances caused by drug injection for metabolic changes, it is recommended to extract the steady-state oxygen consumption rate from the latter half of the residence phase, for example, by regressing from 60 seconds after the start of residence to 120 seconds after the end of residence to obtain the post-drug oxygen consumption rate.

[0036] Subsequently, the mitochondrial uncoupling agent solution was perfused. A cyano-p-trifluoromethoxyphenylhydrazine solution at a concentration of 1 μmol / L was suitable. The injection time was 60 seconds, and the residence time was 240 seconds, followed by rinsing with culture buffer for 120 seconds. After the uncoupling agent eliminated the mitochondrial membrane potential gradient, the electron transport chain operated at a higher flux under conditions without ATP synthesis coupling limitations, resulting in an increased oxygen consumption rate and a peak value. Therefore, the peak oxygen consumption rate measured after perfusion with the mitochondrial uncoupling agent solution was taken as the maximum respiratory capacity. Peak value extraction could be performed using a sliding window to determine the maximum slope, with a window length of 60 seconds to reduce peak fluctuations caused by single-point noise.

[0037] Finally, an electron transport chain complex inhibitor solution is infused. This solution can be a mixture of rotenone and antimycin A, with both rotenone and antimycin A concentrations set at 0.5 μmol / L. The infusion time is 60 seconds, and the residence time is 180 seconds. This mixed inhibitor blocks mitochondrial respiration, and residual oxygen consumption originates from non-mitochondrial pathways. Therefore, the residual oxygen consumption after infusion of the electron transport chain complex inhibitor solution is considered as non-mitochondrial oxygen consumption. Residual oxygen consumption can be directly calculated based on the steady-state oxygen consumption rate after drug administration, or the average rate of decrease in oxygen concentration during the inhibitor residence phase can be used as the residual oxygen consumption. At this point, a set of metabolic functional characteristic parameters is obtained for each circulating tumor cell. This set includes basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption. A one-to-one correspondence is established between these parameters and the corresponding cylindrical microcavity number of the circulating tumor cell to avoid sample misclassification during subsequent release and sorting.

[0038] Regarding alternative implementation methods, to further improve phase delay stability under strong background noise, the phase-locked demodulation accumulation window can be extended from 100 modulation cycles to 300 modulation cycles, while the low-pass filter passband can be narrowed from 20 Hz to 10 Hz. This comes at the cost of a decrease in output update speed, but is still sufficient for second-level metabolic changes. If phase delay bias is found in certain cylindrical microcavities due to membrane thickness variations, a partitioned lookup table can be established for different regions during the calibration phase. For example, the array can be divided into four quadrants for separate calibration, and the corresponding lookup table can be selected based on the microcavity location during runtime, improving the consistency of oxygen concentration values. If greater emphasis is placed on cell viability and subsequent omics quality, the residence time of metabolic regulatory reagents can be shortened, for example, oligomycin residence time can be shortened to 120 seconds and uncoupling agent residence time to 180 seconds, to reduce the total drug exposure time. Furthermore, a buffer rinse time of 180 seconds can be added after each drug treatment to reduce the impact of residues on subsequent omics.

[0039] After acquiring the metabolic function characteristic parameters, the process proceeds to releasing circulating tumor cells and collecting gene expression profiles, protein expression profiles, and metabolite profiles. During release, the perfusion of metabolic regulation reagents is stopped, and the single-cell microcavity capture array area is flushed with culture buffer at a flow rate of 20 μL / min for 300 seconds to reduce drug concentrations inside and outside the cylindrical microcavities, minimizing interference with subsequent sequencing and mass spectrometry. Subsequently, positive pressure gas is introduced into the negative pressure chamber to release the negative pressure adsorption at the micropores. The positive pressure can be 10 kPa to 20 kPa, with 15 kPa used in this example. Simultaneously, cell release buffer is perfused into the single-cell microcavity capture array area. The cell release buffer can be phosphate-buffered saline containing 1% bovine serum albumin, at a flow rate of 30 μL / min for 60 seconds. The positive pressure changes the fluid flow from "inhalation" to "expulsion" in the micropore direction. Circulating tumor cells within the cylindrical microcavities leave the cylindrical microcavities under the combined action of upward local flow and mainstream shear traction, entering the single-cell sorting and collection area. To maximize the release recovery rate, microscopic observation can be performed simultaneously during the release process to confirm that the cylindrical microcavities change from the "occupying state" to the "empty state". After the release is completed, a 60-second buffer rinse is performed to completely remove the released cells from the array area.

[0040] Circulating tumor cells were individually aspirated from the single-cell sorting and collection area using micromanipulation capillaries and transferred to individual reaction tubes. To stably correlate metabolic functional parameters with subsequent omics data, the individual reaction tubes were numbered and barcoded before transfer, with the numbering rules corresponding one-to-one with the cylindrical microcavity numbers; for example, the row and column coordinates of the cylindrical microcavities were converted into individual reaction tube numbers. The inner diameter of the micromanipulation capillaries could be 20 to 30 micrometers, and the aspiration volume could be 0.5 to 1.5 microliters, with 1 microliter used as an example. During aspiration, the morphology of the target cells was first confirmed in the single-cell sorting and collection area using bright-field imaging. Then, the capillary tip was aligned with the target cells, and aspiration was performed slowly under negative pressure, with the negative pressure variation controlled within 5 kPa to reduce cell shear damage. After each transfer, the capillary was rinsed three times with washing buffer to reduce the probability of cross-contamination, with a washing volume of 5 microliters each time.

[0041] Gene expression profiling data were acquired using a single-cell transcriptome sequencing workflow. One possible implementation involves adding a single circulating tumor cell from an independent reaction tube to lysis buffer and immediately protecting it with ribonucleic acid (RNA). The lysis buffer volume can be 2 μL. After lysis, reverse transcription primers and reverse transcriptase are added for reverse transcription at 42°C for 90 minutes. To improve the detection sensitivity at the single-cell level, the complementary DNA after reverse transcription is amplified. The number of amplification cycles can be 18 to 24, with 22 cycles used in the example. The amplified products are purified and used to construct sequencing libraries, with fragment lengths controlled between 300 and 600 base pairs. For high-throughput sequencing, the read length per cell can be 150 base pairs, and the sequencing data volume can be 1,000,000 to 3,000,000 reads, with 2,000,000 reads used in the example to cover more low-abundance transcripts. Subsequently, the sequencing reads were de-adapted, aligned, deduplicated, and each gene was counted to generate gene expression profile data. To ensure that the gene expression profile data could be aligned with the metabolic function characteristic parameter set, the data file name was consistent with the individual reaction tube number, and information such as the cylindrical microcavity number, acquisition start and end time, and exposure time of metabolic regulation reagents was retained in the metadata table.

[0042] Protein expression profiling data were acquired using a single-cell proteomics mass spectrometry (SMS) workflow. To keep single-cell protein loss within acceptable limits, low-adsorption tubing was used in individual reaction tubes, and the total reaction volume was kept below 10 μL. An example procedure involved first adding denaturing buffer to allow for thorough protein development at 60°C for 30 minutes, followed by the addition of a reducing agent and alkylating agent to treat disulfide bonds, and then enzymatic digestion with trypsin at 37°C for 12 hours. The digested peptides were desalted using a solid-phase extraction microcolumn and then identified and quantified by liquid chromatography-tandem mass spectrometry (LC-MS / MS). The LC gradient time was 30 to 60 minutes, with an example of 45 minutes. To improve single-cell proteomics coverage, a data-independent acquisition strategy could be used, or the same cell peptide could be injected twice and the identification results merged during data processing. The final output protein expression profiling data included protein identifiers and corresponding abundance values, linked to the individual reaction tube numbers.

[0043] Metabolite profiling data were acquired using a single-cell metabolomics mass spectrometry (MS / MS) detection workflow. Metabolites are time-sensitive; therefore, it is recommended to terminate metabolism and extract small molecules as soon as possible after release and transfer. One possible implementation is to add a pre-cooled methanol solution (approximately 20 μL) to a separate reaction tube, gently mix, and place on ice for 5 minutes to rapidly halt intracellular metabolic reactions. Water and acetonitrile are then added to form an extraction system. After centrifugation, the supernatant is collected for liquid chromatography-tandem mass spectrometry (LC-MS / MS) analysis. LC-MS can employ a hydrophilic column to separate polar metabolites, while a reversed-phase column separates hydrophobic metabolites, forming a two-injection coverage strategy. Mass spectrometry detection can be performed once in positive ion mode and once in negative ion mode to cover small molecules with different ionization characteristics. Data processing includes peak extraction, retention time calibration, isotope deconvolution, and matching with a standard library to obtain metabolite identifiers and corresponding abundance values, thus generating metabolite profiling data.

[0044] After simultaneously acquiring gene expression profiles, protein expression profiles, and metabolite profiles, metabolic function feature parameter sets are fused with the three types of omics data to form a multi-omics joint feature vector. Before fusion, label alignment is performed: gene expression profile data uses gene symbols or transcript identifiers as keys, protein expression profile data uses protein identifiers as keys, and metabolite profile data uses metabolite names and library identifiers as keys. A unified feature index table is established for each cell. To control vector dimensionality and reduce noise, features are selected from a subset of genes with the highest variability in gene expression profile data (e.g., 2000 genes with the highest variability); features are selected from proteins with a high intercellular reproducibility rate in protein expression profile data (e.g., 500 proteins with high detection rates); and features are selected from metabolite profile data with qualitative confidence levels reaching the library matching threshold (e.g., 200 metabolites with high matching confidence). Subsequently, standardization is performed within the same metabolic phenotype layer. Standardization can be achieved using... ,in This represents the standardized numerical value. This represents the raw numerical value of a cell on a certain molecular index. This represents the arithmetic mean of this molecular index across all cells within the metabolic phenotype layer. This represents the standard deviation of this molecular indicator for all cells within the metabolic phenotype layer. Performing standardization within the metabolic phenotype layer, rather than across all cells, helps reduce the impact of overall scale differences between different metabolic phenotype layers on clustering and distance calculations, while preserving finer-grained molecular differences within the same metabolic phenotype layer. After standardization, the three types of standardized data are concatenated end-to-end in the order of gene-protein-metabolite to obtain a multi-omics joint feature vector for each circulating tumor cell. Metabolic functional feature parameter sets are then appended to the vector metadata, enabling subsequent tumor feature atlas construction to utilize both molecular features and oxygen metabolism phenotype features.

[0045] Step 4: Construct circulating tumor cell subpopulations based on multi-omics joint feature vectors, infer the subpopulation evolution directed graph, extract landmark molecular features, and integrate and output tumor feature atlas.

[0046] After generating the multi-omics joint feature vector for each circulating tumor cell, the multi-omics joint feature vectors of all circulating tumor cells are first aggregated along the cell dimension to obtain the multi-omics joint feature matrix. Each row of the multi-omics joint feature matrix corresponds to one circulating tumor cell, and each column corresponds to one molecular indicator. The molecular indicators are arranged in the following order: molecular indicators corresponding to gene expression profile data, molecular indicators corresponding to protein expression profile data, and molecular indicators corresponding to metabolite profile data, concatenated end to end. To ensure that the subsequent principal component analysis is more balanced for features of different dimensions, the multi-omics joint feature matrix maintains a numerical system consistent with the normalization process within the metabolic phenotype layer before entering the principal component analysis. A consistent imputation strategy is adopted for molecular indicators with missing values, such as replacing missing values ​​with the normalized value of 0 within the same metabolic phenotype layer, so that the missing values ​​do not introduce additional bias and do not destroy the relative differences between different cells.

[0047] Principal component analysis (PCA) is used to map high-dimensional multi-omics joint feature vectors to a low-dimensional, computable distance-reduced feature space, thereby reducing the impact of noise and redundancy on similarity metrics. One feasible procedure is to first center the multi-omics joint feature matrix column-wise; the centered matrix is ​​denoted as... ,in This represents the centered multi-omics joint feature matrix. The centering process involves calculating the column mean for each column of molecular indicators and subtracting the column mean from each value in that column, ensuring each column is centered at 0. The covariance matrix is ​​then calculated. Writing the covariance matrix ,in Represents the covariance matrix. Representation matrix transpose, Indicates the number of circulating tumor cells. This represents the correction term for the degrees of freedom. (Regarding the covariance matrix) Eigenvalue decomposition yields a set of eigenvectors and a set of eigenvalues. The eigenvalues ​​are sorted from largest to smallest, and the eigenvector corresponding to each eigenvalue is considered a principal component direction. Larger eigenvalues ​​indicate that the principal component direction can explain more variance and better preserve intercellular differences. To determine the number of retained principal components, the variance contribution rate is calculated and summed. The variance contribution rate can be written as... ,in Indicates the first The variance contribution rate of each principal component Indicates the first The eigenvalues ​​corresponding to each principal component This represents the sum of the eigenvalues ​​of all principal components. Then, the cumulative variance contribution rate is calculated. ,in Indicates the preceding The cumulative variance contribution rate of each principal component Indicates the number of principal components retained. When When the preset explanatory power threshold is reached, the accumulation stops and the previous values ​​are retained. The principal components constitute the dimensionality-reduced feature space. The preset explanatory power threshold can be selected from 0.80 to 0.95; the example uses 0.90. This preserves most of the variations related to the subtype while compressing more noise dimensions, thus making distance calculations more stable. Taking an example dataset containing 120 circulating tumor cells, each with a multi-omics joint feature vector of dimension 2700, it is common practice to retain 40 to 60 principal components to achieve a cumulative variance contribution rate of 0.90. In this example, 50 principal components are retained to construct the dimensionality-reduced feature space. Projecting the multi-omics joint feature vector of each circulating tumor cell onto the dimensionality-reduced feature space yields the dimensionality-reduced feature vector for that cell, written as... ,in Represents the dimensionality-reduced feature vector. This represents the centered multi-omics joint feature vector of the cell. Indicates from the previous The projection matrix is ​​formed by concatenating the principal component eigenvectors column by column.

[0048] After obtaining the dimensionality-reduced feature vectors of all circulating tumor cells, a similarity adjacency graph of circulating tumor cells is constructed, forming a nearest-neighbor undirected graph. The construction of the nearest-neighbor undirected graph prioritizes local similarity: each cell is connected only to a small subset of its closest neighbors, avoiding noise diffusion from fully connected graphs and reducing subpopulation fragmentation caused by isolated points. Specifically, the Euclidean distance is calculated for the dimensionality-reduced feature vectors of any two circulating tumor cells. The Euclidean distance is written as... ,in Represents cells With cells Euclidean distance, Represents cells The dimensionality reduction eigenvectors in the th dimension The values ​​on each principal component dimension Represents cells The dimensionality reduction eigenvectors in the th dimension The values ​​on each principal component dimension This represents summing over all retained principal component dimensions. The Euclidean distances between a given cell and all other cells are sorted in ascending order. A preset number of nearest neighbors are selected as neighboring cells, and undirected edges are established. The preset number of nearest neighbors can be between 10 and 30; this example uses 15. The reason for choosing 15 is that, with a sample size in the hundreds, 15 ensures that the undirected graph of nearest neighbors generally maintains overall connectivity, while avoiding forcibly pulling obviously dissimilar, distant cells into the same local neighborhood. To reduce accidental connections caused by "one-way nearest neighbors," a symmetry strategy can be optionally adopted: as long as the cell... cells Select as a neighbor, or cell cells Selecting a neighbor as a neighbor means establishing an undirected edge between the two; this improves the robustness of the nearest neighbor undirected graph in the sparse case.

[0049] Constrained community detection based on metabolic continuity is performed on a nearest-neighbor undirected graph to obtain circulating tumor cell subpopulations. The core idea of ​​constrained community detection is to simultaneously satisfy two conditions: one is that the joint feature vectors of multi-omics are similar in the dimensionality-reduced feature space, and the other is that the metabolic function feature parameter set remains continuous among adjacent cells. This binds "molecular similarity" and "metabolic compatibility" into the same subpopulation standard, reducing false clustering caused solely by molecular noise. Specifically, a community assignment list and a queue of nodes to be expanded are established. The community assignment list records whether each circulating tumor cell has been assigned to a community and which community it has been assigned to; the queue of nodes to be expanded is used to expand the community boundaries in a first-in, first-out (FIFO) order. A node not assigned to any community is selected from the nearest-neighbor undirected graph as a seed node. A new community is created, and the seed node is assigned to this new community. Simultaneously, the seed node is added to the queue of nodes to be expanded. Then, the head node of the queue of nodes to be expanded is taken as the current expansion node, and all neighbor nodes of the current expansion node are traversed. For each neighboring node, if it has already been assigned to a community, skip that node; otherwise, perform a difference analysis on each of the metabolic function characteristic parameters between the neighboring node and the currently expanding node. The metabolic function characteristic parameter set includes four parameters: basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption. Calculate the absolute value of the difference between each of the four parameters, written as ______. ,in This represents the numerical value of a certain metabolic function characteristic parameter of a neighboring node. This represents the numerical value of the metabolic function characteristic parameter corresponding to the current extended node. This indicates taking the absolute value. Only when the absolute values ​​of the differences among all four parameters are less than the preset metabolic continuity threshold is a neighboring node assigned to the community of the currently expanding node and added to the queue of nodes to be expanded. The preset metabolic continuity threshold needs to balance "not being too strict, leading to community fragmentation" and "not being too lenient, leading to confounding different metabolic phenotypes." One feasible approach is to first calculate the median absolute deviation for each metabolic function characteristic parameter of all circulating tumor cells, and then use the median absolute deviation to generate the threshold. The median absolute deviation is written as... ,in Indicates the absolute deviation of the median. This represents the set of values ​​for a specific metabolic function parameter across all circulating tumor cells. This represents the median function. The threshold can be between 2 and 4 times the absolute deviation of the median; in this example, it is 3 times. The advantage of this is that the absolute deviation of the median is not sensitive to outliers, and the threshold better represents the "typical fluctuation range," preventing the entire community expansion condition from being skewed when individual cells exhibit extreme metabolic values. Community expansion continues until the queue of nodes to be expanded is empty, indicating that the community can no longer expand under the metabolic continuity constraint. Then, it checks if there are still nodes in the nearest undirected graph that have not been assigned to any community. If so, one is selected as the new seed node, and the above process is repeated until all nodes have been assigned to communities. Each community corresponds to a circulating tumor cell subpopulation, which is identified by consecutive numbers, for example, starting from 1.

[0050] After obtaining circulating tumor cell subsets, a directed graph of subset evolution is inferred. First, metabolic activity indices are calculated for each circulating tumor cell subset. These indices are expressed as the average basal oxygen consumption rate within the subset, and are written as... ,in Indicating circulating tumor cell subsets The metabolic activity index value, This represents the summation of the basal oxygen consumption rates of all member cells within this subgroup. This represents the number of member cells within the subpopulation. The reason for using the average basal oxygen consumption rate (BAOC) is that it directly reflects the stable respiratory level, is obtained before the intervention of metabolic regulators, is less affected by external disturbances, and is more consistent as a comparison indicator between subpopulations. Then, all undirected edges in the nearest neighbor undirected graph are traversed. If two cells connected by an undirected edge belong to different circulating tumor cell subpopulations, an adjacency relationship is recorded between the two subpopulations. The introduction of adjacency relationships ensures that evolutionary inference is established only between subpopulations that are close to each other in molecular space, avoiding overextension caused by directly connecting edges between completely dissimilar subpopulations based on the magnitude of metabolic activity index values. For any pair of adjacent circulating tumor cell subpopulations, the magnitude of the metabolic activity index values ​​of the two subpopulations is compared, and a directed edge is established from the circulating tumor cell subpopulation with the smaller metabolic activity index value to the circulating tumor cell subpopulation with the larger metabolic activity index value, thus forming a directed graph of subpopulation evolution. To reduce directional instability when the metabolic activity index values ​​of two subgroups are very close, a minimum difference threshold can be introduced. For example, if the difference between the two metabolic activity index values ​​is less than 0.05, a directed edge is not established. Alternatively, the number of undirected edges across subgroups can be used as a second criterion: directions with a larger number of undirected edges across subgroups are preferred to be retained to avoid direction reversal caused by small numerical fluctuations.

[0051] refer to Figure 1 , Figure 1This paper demonstrates the core principle of a circulating tumor cell (CTC) subpopulation boundary identification method based on a nearest-neighbor undirected graph. The graph presents the distribution patterns and interconnections of 120 CTCs in a two-dimensional feature space after dimensionality reduction via principal component analysis (PCA). The horizontal axis represents Principal Component 1 (PC1), and the vertical axis represents Principal Component 2 (PC2). These two principal components are the most important variation directions extracted from a multi-omics joint feature vector composed of gene expression profiles, protein expression profiles, and metabolite profiles. All cell points in the graph are assigned to four different CTC subpopulations based on their metabolic phenotypes, labeled red, blue, green, and orange, respectively. Each subpopulation exhibits a relatively concentrated regional distribution in the dimensionality-reduced feature space, but there is a certain spatial proximity between the subpopulations. This proximity provides a topological basis for subsequent subpopulation evolution inference.

[0052] The construction of the nearest-neighbor undirected graph follows a principle prioritizing local similarity. For any circulating tumor cell, the system first calculates the Euclidean distance between this cell and all other cells in the dimensionality-reduced feature space. The Euclidean distance comprehensively reflects the combined differences in the dimensionality of the multi-omics joint feature vectors after dimensionality reduction by principal component analysis while preserving the principal component dimensions. Subsequently, these distances are sorted in ascending order, and the 15 nearest neighbors are selected as the cell's neighbors. An undirected edge is then established between this cell and each of its neighbors. This connection strategy based on a fixed number of neighbors ensures that each cell can be connected to a sufficient number of similar cells, while avoiding forcibly pulling obviously dissimilar, distant cells into the same local neighborhood. Two types of connectors can be observed in the figure. One type is the thin gray line, which represents intra-subpopulation connectivity, meaning that the cells at both ends of the connector belong to the same circulating tumor cell subpopulation. This type of connectivity forms a dense network structure within each subpopulation, reflecting the high similarity of cells in the same subpopulation in terms of multi-omics characteristics. The other type is the purple dashed line, which represents inter-subpopulation connectivity, meaning that the cells at both ends of the connector belong to different circulating tumor cell subpopulations. This type of connectivity appears in the boundary region between different subpopulations, indicating that there is a gradual transition rather than abrupt separation between subpopulations in the feature space.

[0053] The nearest-neighbor undirected graph serves a dual function in this invention. Firstly, it acts as a topological constraint framework for community detection based on metabolic continuity. The community detection algorithm selects unassigned nodes from the nearest-neighbor undirected graph as seed nodes to create new communities, then expands along undirected edges to neighboring nodes. Neighboring nodes are only included in the same community if their metabolic function feature parameter sets satisfy the continuity constraint with the currently expanded node. This expansion mechanism ensures that each circulating tumor cell subpopulation maintains internal consistency in metabolic function while sharing similar multi-omics features, thus avoiding the masking of metabolic heterogeneity caused by clustering based solely on molecular features. Secondly, it provides adjacency criteria for establishing edges in the directed graph of subpopulation evolution. Evolutionary directions can only be established between two circulating tumor cell subpopulations connected by undirected edges across subpopulations in the nearest-neighbor undirected graph. This constraint ensures that evolutionary inference is only performed between subpopulations that are close to each other in molecular space, avoiding over-extrapolation caused by directly connecting edges between completely dissimilar subpopulations based on the magnitude of metabolic activity index values.

[0054] The spatial distribution of the four circulating tumor cell subpopulations in the figure exhibits typical metabolic phenotypic stratification. The red subpopulation, located in the lower left region, shows high values ​​for both basal oxygen consumption (BAC) and ATP synthesis-related oxygen consumption (ATP-related oxygen consumption), corresponding to the metabolic characteristics of a highly oxidative phosphorylation-active layer. The blue subpopulation, located in the lower right region, has a lower BAC and the quotient of ATP-related oxygen consumption to BAC is less than a preset threshold, belonging to the glycolysis-dependent layer. The green subpopulation, located in the upper left region, shows a large difference between its maximum respiratory capacity and BAC, indicating ample respiratory reserves, belonging to the metabolically adequate layer. The orange subpopulation, located in the upper right region, exhibits intermediate states in multiple metabolic parameters, representing a mixed metabolic phenotype. The connectivity density of the nearest neighbor undirected graph reveals numerous inter-subpopulation connections between the red and green subpopulations, and a certain number between the blue and orange subpopulations. However, relatively few inter-subpopulation connections exist between the red and blue subpopulations. This connectivity pattern suggests differences in the ease of transition between different metabolic phenotypes.

[0055] The construction parameters of the nearest neighbor undirected graph have a significant impact on subgroup boundary identification. Setting the preset number of nearest neighbors to 15 ensures overall connectivity of the graph when the sample size is in the hundreds, while preventing clearly dissimilar, distant cells from being forcibly pulled into the same local neighborhood. If the preset number of nearest neighbors is too small, such as 5, the graph is prone to having multiple connected components, making subgroup identification impossible globally. If the preset number of nearest neighbors is too large, such as 30, the graph introduces a large number of long-range connections, making cross-subgroup connections too dense and weakening the clarity of subgroup boundaries. The graph also shows that the number of neighbor connections for each cell is not completely equal. This is due to a symmetry strategy: an undirected edge is established between cell A and cell B whenever cell A selects cell B as a neighbor or cell B selects cell A as a neighbor. This strategy improves the robustness of the nearest neighbor undirected graph in sparse cases, ensuring that important local similarities are not missed due to the randomness of unidirectional nearest neighbor selection.

[0056] After forming the subpopulation evolution directed graph, evolutionary stage numbers are assigned to each circulating tumor cell subpopulation to give the tumor feature atlas a readable stage structure. First, the circulating tumor cell subpopulation with the lowest metabolic activity index value is identified in the subpopulation evolution directed graph as the starting subpopulation, and its evolutionary stage number is set to 1. Then, a stage assignment queue is established, and the starting subpopulation is added to the queue. The following process is executed iteratively: the first subpopulation in the stage assignment queue is taken as the current subpopulation; all directed edges in the subpopulation evolution directed graph originating from the current subpopulation are traversed; for each directed edge pointing to a target subpopulation, the target subpopulation is checked. If the target subpopulation has not yet been assigned an evolutionary stage number, its evolutionary stage number is set to the current subpopulation's evolutionary stage number plus 1, and the target subpopulation is added to the stage assignment queue; when the target subpopulation has already been assigned an evolutionary stage number, the smaller assigned number is retained, and optionally another candidate number is recorded for subsequent consistency checks. This queue expansion method ensures that the nearest neighbor stages starting from the initial subpopulation are assigned values ​​first, forming an overall stage progression from low to high metabolic activity. If there are multiple incoming edges in the subpopulation evolution directed graph that cause stage conflicts, a "shortest path stage" strategy can be optionally introduced: the final stage of the target subpopulation is taken as the shortest path length from the initial subpopulation to the target subpopulation plus 1, in order to reduce stage expansion caused by local branches.

[0057] After determining the circulating tumor cell subpopulations and evolutionary stages, characteristic molecular features of each subpopulation are extracted to provide interpretable molecular descriptions when outputting tumor feature atlases. For each molecular indicator in the multi-omics joint feature vector, the mean within and outside the subpopulation are calculated, and the expression ratio is also calculated. The mean within the subpopulation is written as... ,in This represents the mean within a subgroup. This represents the sum of the values ​​of all member cells within this subgroup on this molecular index. This indicates the number of member cells within a subpopulation. The mean outside a subpopulation is written as... ,in This represents the out-of-group mean. Indicates the number of cells outside the subpopulation. The expression ratio is written as... ,in This represents the expression ratio. The reason for using a ratio instead of a difference is that ratios are less sensitive to changes in the overall scale, especially maintaining relatively comparable upregulation or downregulation directions after splicing different omics. Subsequently, the positive and negative biomarker sets are determined based on the expression ratio and preset upregulation and downregulation thresholds. The preset upregulation threshold can be selected from 1.5 to 3.0, with an example of 2.0; the preset downregulation threshold can be selected from 0.3 to 0.7, with an example of 0.5. Genes, proteins, or metabolites corresponding to molecular indicators with expression ratios greater than the preset upregulation threshold are included in the positive biomarker set, and those corresponding to molecular indicators with expression ratios less than the preset downregulation threshold are included in the negative biomarker set. To avoid instability caused by excessively small out-of-group averages, a lower limit truncation can be introduced for the out-of-group averages, for example, if the out-of-group average is less than 0.01, it is included in the calculation as 0.01, thus preventing a single rare feature from being amplified into a falsely high ratio biomarker.

[0058] The integrated output of the tumor feature atlas simultaneously stores structural and content information, facilitating retrieval and subsequent analysis. A feasible output organization includes two main tables and one relational table. The first main table records, at the circulating tumor cell level: cell ID, multi-omics joint feature vector, dimensionality-reduced feature vector, circulating tumor cell subpopulation ID, evolutionary stage number, and metabolic function feature parameter set. The second main table records, at the circulating tumor cell subpopulation level: circulating tumor cell subpopulation ID, metabolic activity index value, positive biomarker set, negative biomarker set, and a list of member cells within the circulating tumor cell subpopulation. The relational table records the directed edge connections in the subpopulation evolution directed graph. Each record includes the starting circulating tumor cell subpopulation ID, the ending circulating tumor cell subpopulation ID, the number of undirected edges across subpopulations in the adjacency relationship between the two subpopulations, and the starting and ending metabolic activity index values. The output format can be either a comma-separated file or a structured text file. The comma-separated file facilitates reading by statistical tools, while the structured text file preserves the hierarchical structure of set and relational fields. To ensure comparability across batches, the output can include the generation method and corresponding values ​​for the number of retained principal components, preset nearest neighbor number, preset explanatory power threshold, preset upregulation threshold, preset downregulation threshold, and preset metabolic continuity threshold for principal component analysis, so that the same process can be reproduced in different experimental batches.

[0059] Regarding optional implementation methods, if a small sample size leads to multiple connected components in the nearest neighbor undirected graph, the preset nearest neighbor number can be increased from 15 to 25 to improve connectivity before performing constrained community detection based on metabolic continuity. If greater emphasis is placed on the constraint strength of metabolic function characteristic parameter groups, the preset metabolic continuity threshold can be adjusted from 3 times the absolute deviation of the median to 2 times, making community expansion more stringent and obtaining more "metabolic consistent" circulating tumor cell subpopulations. If more attention is paid to continuous evolution at the molecular level and the direction is not desired to be determined solely by the basal oxygen consumption rate, the metabolic activity index value can be replaced from the "average value of the basal oxygen consumption rate" to the "average value of the difference between the maximum respiratory capacity and the basal oxygen consumption rate," which is closer to the potential evolutionary drivers related to respiratory reserve. After the replacement, the same method can be used to establish directed edges between adjacent circulating tumor cell subpopulations, from those with smaller metabolic activity index values ​​to those with larger metabolic activity index values, to maintain process consistency.

[0060] refer to Figure 2 This heatmap presents the differences in four circulating tumor cell subpopulations across four key metabolic functional parameters in matrix form. The horizontal axis, from left to right, represents subpopulation 1, subpopulation 2, subpopulation 3, and subpopulation 4. The vertical axis, from top to bottom, represents basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption. Each matrix cell is coded using a red-yellow-green color scheme: red indicates a positive and large standardized value for that metabolic parameter, green indicates a negative and large absolute value, and yellow indicates a standardized value close to zero. The color bar scale ranges from -1.5 to +1.5, covering a typical range of standardized values. Each cell also includes the specific standardized value, facilitating quantitative comparisons of metabolic functional differences between subpopulations. The acquisition of these metabolic functional parameter sets relied on a phase-sensitive single-cell oxygen metabolism dynamic monitoring protocol. The basal oxygen consumption rate (BAC) was obtained by linear regression of the oxygen concentration over time within the post-capture stabilization window; the negative regression slope represents the BAC. This parameter reflects the steady-state respiration level of circulating tumor cells (CMCs) without external metabolic regulatory intervention. ATP synthesis-related oxygen consumption was calculated by the difference in oxygen consumption rates before and after perfusion with an ATP synthase inhibitor; this parameter quantifies the portion of oxygen consumed for ATP synthesis during mitochondrial oxidative phosphorylation. Maximum respiratory capacity was characterized by the peak oxygen consumption rate measured after perfusion with a mitochondrial uncoupling agent; this parameter reflects the maximum electron transport flux of CMCs under conditions without ATP synthesis coupling limitations. Non-mitochondrial oxygen consumption was determined by the residual oxygen consumption rate after perfusion with an electron transport chain complex inhibitor; this parameter represents oxygen consumption from non-mitochondrial pathways, such as peroxisomal reactions or cytoplasmic redox reactions. These four parameters characterize the energy metabolism status of CMCs from different dimensions, providing a quantitative basis for metabolic phenotypic stratification.

[0061] Standardization is performed within the same metabolic phenotype layer, rather than across all cells. This design helps reduce the impact of overall scale differences between different metabolic phenotype layers on clustering and distance calculations, while preserving finer-grained molecular differences within the same metabolic phenotype layer. The specific standardization formula is z = v minus μ divided by σ, where z represents the standardized value, v represents the original value of a cell for a certain molecular indicator, μ represents the arithmetic mean of all cells within that metabolic phenotype layer for that molecular indicator, and σ represents the standard deviation of all cells within that metabolic phenotype layer for that molecular indicator. The standardized values ​​presented in the heatmap are actually the result of summarizing and then standardizing across metabolic phenotype layers to compare the metabolic characteristics of different subpopulations at a uniform scale. The heatmap clearly shows significant differences in the metabolic function parameters of the four circulating tumor cell subpopulations. Subgroup 1 showed positive standardized values ​​of 0.8 and 1.2 for both basal oxygen consumption rate and ATP synthesis-related oxygen consumption, respectively, appearing in a red hue. Its standardized value for maximum respiratory capacity was 1.0, also in red, while its standardized value for non-mitochondrial oxygen consumption was 0.3, appearing in a yellow-green hue. This parameter combination indicates that subgroup 1 cells primarily rely on mitochondrial oxidative phosphorylation for energy production, with high ATP synthesis coupling efficiency in this process, consistent with the typical characteristics of highly oxidative phosphorylation-active layers. Subgroup 2 showed negative standardized values ​​of -0.9, -1.1, and -0.8 for basal oxygen consumption rate, ATP synthesis-related oxygen consumption, and maximum respiratory capacity, respectively, appearing in a green hue. Its standardized value for non-mitochondrial oxygen consumption was 0.2, appearing in a yellow hue. This parameter combination indicates that subgroup 2 cells have generally lower mitochondrial respiratory activity and rely more on the glycolysis pathway for energy supply, consistent with the metabolic characteristics of glycolysis-dependent layers.

[0062] Subpopulation 3 exhibits an extremely high standardized value of 1.5 for maximum respiratory capacity (deep red), while the standardized values ​​for basal oxygen consumption rate (BAC) and ATP synthesis-related oxygen consumption (ATP-related oxygen consumption) are 0.3 and 0.5 respectively (light red / yellow), and the standardized value for non-mitochondrial oxygen consumption is 0.1 (close to zero, yellow). The key characteristic of this parameter combination is the large difference between maximum respiratory capacity and BAC, indicating that this subpopulation has ample respiratory reserves and can rapidly increase mitochondrial respiratory flux when faced with sudden increases in energy demand or stress conditions, consistent with the functional positioning of the metabolic reserve layer. Subpopulation 4 shows intermediate metabolic function parameters: the standardized value for BAC is 0.1 (close to zero, yellow), the standardized value for ATP synthesis-related oxygen consumption is -0.2 (light green), the standardized value for maximum respiratory capacity is 0.2 (yellow), and the standardized value for non-mitochondrial oxygen consumption is 0.9 (red). This parameter combination indicates that subpopulation 4 cells exhibit a relatively balanced distribution between mitochondrial respiration and non-mitochondrial oxygen consumption, and the ATP synthesis coupling efficiency is slightly below average, representing a mixed metabolic phenotype.

[0063] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for constructing tumor feature atlases based on CTC enrichment and multi-omics analysis, characterized in that, The method includes: Step 1: The whole blood sample to be tested is sequentially subjected to helical inertial focusing and deterministic lateral displacement sorting, and the single-cell microcavity localization and capture of circulating tumor cells is completed in the single-cell microcavity capture array area; Step 2: Obtain the oxygen molecule concentration value over time for each captured circulating tumor cell and obtain a set of metabolic function characteristic parameters; Step two includes: using a direct digital frequency synthesis chip to generate a sinusoidal modulated voltage signal with a fixed frequency; the sinusoidal modulated voltage signal drives an ultraviolet light-emitting diode to generate modulated excitation light through a voltage-to-current conversion circuit; the modulated excitation light is collimated by a collimating lens and then irradiates the single-cell microcavity trapping array region to excite the oxygen-sensitive phosphorescence sensing membrane to generate red-band phosphorescence emission; the phosphorescence emission is collected by an objective lens and filtered by a red bandpass filter to remove excitation light scattering before being input into a photomultiplier tube; the photocurrent signal output by the photomultiplier tube is converted into a voltage signal by a transimpedance amplifier and then digitally sampled by an analog-to-digital converter; the digitally sampled phosphorescence signal is subjected to phase-locked demodulation processing in a field-programmable gate array (FPGA); the phase-locked demodulation processing includes generating a reference sine signal and a cosine signal obtained by shifting the reference sine signal by 90 degrees; the digitally sampled phosphorescence signal is multiplied point-by-point with the reference sine signal and the cosine signal respectively, and then passed through a finite impulse response (FIR) low-pass filter to obtain a DC in-phase component signal and a DC quadrature component signal respectively; Step 3: Release circulating tumor cells and obtain gene expression profile data, protein expression profile data and metabolite profile data respectively, and combine them with metabolic function characteristic parameter groups to form a multi-omics joint feature vector; Step 4: Construct circulating tumor cell subpopulations based on multi-omics joint feature vectors, infer the subpopulation evolution directed graph, extract landmark molecular features and integrate and output tumor feature atlas; A constrained community detection based on metabolic continuity was performed on a nearest-neighbor undirected graph to obtain circulating tumor cell subsets. The metabolic function characteristic parameter set included four parameters: basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption. The absolute value of the difference between each of the four parameters was calculated, and the absolute value of the difference was written as... ,in This represents the numerical value of a certain metabolic function characteristic parameter of a neighboring node. This represents the numerical value of the metabolic function characteristic parameter corresponding to the current extended node. This means taking the absolute value. Only when the absolute value of the difference between the four parameters is less than the preset metabolic continuity threshold will the neighbor node be assigned to the community to which the current expansion node belongs and the neighbor node be added to the queue of nodes to be expanded.

2. The method according to claim 1, characterized in that, Step one, helical inertial focusing, includes: injecting the whole blood sample to be tested into the helical inertial focusing channel. The cross-section of the helical inertial focusing channel is rectangular and the width remains constant along the helical direction. When the whole blood sample to be tested flows in the helical inertial focusing channel, it is simultaneously subjected to inertial lift and Dean drag force. Cells with a diameter greater than a preset size threshold form a focusing band at the equilibrium position of the two forces and move along the inner wall of the channel. Red blood cells and platelets migrate to the outside of the channel with the Dean secondary flow and are discharged through the waste liquid outlet, thereby obtaining a primary enriched suspension containing white blood cells and circulating tumor cells.

3. The method according to claim 2, characterized in that, The deterministic lateral displacement sorting in step one includes: introducing the primary enriched suspension into the deterministic lateral displacement sorting area, setting up a periodically arranged array of micropillars in the deterministic lateral displacement sorting area, wherein adjacent columns of micropillars in the micropillar array have a fixed lateral offset along the flow direction, when the cell diameter is greater than the critical separation diameter, the cell undergoes lateral displacement along the micropillar offset direction after contacting the micropillar and enters the circulating tumor cell collection channel, when the cell diameter is less than the critical separation diameter, the cell passes through the micropillar gap along the streamline direction and enters the leukocyte discharge channel, thereby obtaining the secondary enriched suspension of circulating tumor cells.

4. The method according to claim 3, characterized in that, Step one, single-cell microcavity localization and capture, includes: guiding a secondary enrichment suspension of circulating tumor cells into the single-cell microcavity capture array area; the bottom surface of the single-cell microcavity capture array area is provided with a matrix of cylindrical microcavities, the diameter of each cylindrical microcavity being used to accommodate a single circulating tumor cell; a micropore is opened at the center of the bottom of each cylindrical microcavity and is connected to a negative pressure chamber, so that the circulating tumor cell is drawn into the cylindrical microcavity under negative pressure and the micropore is sealed; after sealing, the negative pressure in the cylindrical microcavity disappears, thus forming a single-cell occupancy; an oxygen-sensitive phosphorescent sensing membrane is pre-deposited on the inner wall surface of each cylindrical microcavity, the oxygen-sensitive phosphorescent sensing membrane being formed by uniformly dispersing platinum porphyrin luminescent molecules in a polystyrene matrix.

5. The method according to claim 4, characterized in that, Step two also includes: determining the phase delay of phosphorescence emission relative to the modulated excitation light based on the DC in-phase component signal and the DC quadrature component signal, and converting the phase delay into the oxygen molecule concentration value in the cylindrical microcavity based on the pre-calibrated quenching correspondence; continuously collecting the oxygen molecule concentration value in each cylindrical microcavity and recording the change curve of the oxygen molecule concentration value over time, and determining the decrease in the oxygen molecule concentration value per unit time as the basal oxygen consumption rate of the corresponding circulating tumor cells.

6. The method according to claim 5, characterized in that, Step two also includes: after the basal oxygen consumption rate is measured, metabolic regulation reagents are sequentially perfused into the single-cell microcavity capture array region through the reagent injection channel of the microfluidic chip. The metabolic regulation reagents include ATP synthase inhibitor solution, mitochondrial uncoupling agent solution, and electron transport chain complex inhibitor solution. The basal oxygen consumption rate before perfusion of ATP synthase inhibitor solution minus the oxygen consumption rate after perfusion of ATP synthase inhibitor solution is taken as ATP synthesis-related oxygen consumption. The peak oxygen consumption rate after perfusion of mitochondrial uncoupling agent solution is taken as maximum respiratory capacity. The residual oxygen consumption after perfusion of electron transport chain complex inhibitor solution is taken as non-mitochondrial oxygen consumption, thereby forming a set of metabolic function characteristic parameters, which includes basal oxygen consumption rate, ATP synthesis-related oxygen consumption, maximum respiratory capacity, and non-mitochondrial oxygen consumption.

7. The method according to claim 6, characterized in that, Step three includes: introducing positive pressure gas into the negative pressure chamber to release the negative pressure adsorption at the micropores and perfusing the single-cell microcavity capture array area with cell release buffer, allowing circulating tumor cells to be released from the cylindrical microcavities and flow out to the single-cell sorting and collection area; in the single-cell sorting and collection area, micromanipulating capillaries are used to aspirate circulating tumor cells one by one and transfer them to independent reaction tubes; for each single circulating tumor cell in each independent reaction tube, single-cell transcriptome sequencing is performed sequentially to obtain gene expression profile data, single-cell proteome mass spectrometry is performed to obtain protein expression profile data, and single-cell metabolome mass spectrometry is performed to obtain metabolite profile data.

8. The method according to claim 7, characterized in that, Step three also includes: calculating the median basal oxygen consumption rate of all circulating tumor cells as the basal oxygen consumption rate cutoff value; calculating the median maximum respiratory capacity of all circulating tumor cells as the maximum respiratory capacity cutoff value; calculating the difference between the maximum respiratory capacity and the basal oxygen consumption rate of all circulating tumor cells and taking the median as the respiratory reserve cutoff value; and stratifying all circulating tumor cells by metabolic phenotype based on metabolic function characteristic parameter groups. Circulating tumor cells with a basal oxygen consumption rate greater than the basal oxygen consumption rate cutoff value and a maximum respiratory capacity greater than the maximum respiratory capacity cutoff value are classified into oxidative phosphorylation active cells. Circulating tumor cells (CTCs) with a basal oxygen consumption rate (BAC) less than the BAC threshold and a quotient of ATP synthesis-related oxygen consumption and BAC less than a preset threshold are classified into a glycolysis-dependent layer. CTCs with a maximum respiratory capacity and a BAC greater than the respiratory reserve threshold are classified into a metabolic reserve-sufficient layer. CTCs belong to multiple metabolic phenotype layers. Within each metabolic phenotype layer, gene expression profiles, protein expression profiles, and metabolite profiles are standardized and concatenated end-to-end in the gene-protein-metabolite order to form a multi-omics joint feature vector.

Images