Experimental method based on depletion-type nk cells as lung cancer treatment target

Abstract

The invention discloses an experimental method based on depletion-type nk cells as a lung cancer treatment target, and relates to the technical field of bioinformatics. Comprising the following steps: S1, data processing; s2, identification of cell subset characteristic genes and cell type annotation; s3, re-clustering the NK cells; s4, NK cell function annotation; s5, differential gene expression analysis; s6, analyzing the difference of activation and inhibition of receptor expression by the RGS1 + NK cells through single cells; s7, performing quasi-time sequence analysis; s8, carrying out cell culture and siRNA transfection; s9, detecting the expression condition of an NK92mi cell activated receptor after the RGS1 is knocked down; s10, carrying out a lactic dehydrogenase release experiment; s11, measuring the IFN-gamma content of the cell supernatant; and S12, virtual screening of drugs. Through analysis and experimental verification, the key function of RGS1 in NK cell depletion in lung cancer is disclosed, and a new and potential molecular target is provided for lung cancer treatment.

Description

Technical Field

[0001] The invention relates to the technical field of bioinformatics, and in particular to an experimental method based on exhaustion-type NK cells as a target for lung cancer treatment. Background Art

[0002] Traditional high-throughput transcriptome sequencing (bulk RNA-seq) mainly reflects the average expression level of the cell population. This technology often makes the expression signals of cell clones with growth advantages dominate, thus masking the differences in cell type-specific expression patterns. This limitation hinders the in-depth exploration of the deep molecular mechanisms of tumors and intratumor heterogeneity. In contrast, single-cell RNA sequencing (scRNA-seq) technology provides a comprehensive and unbiased perspective that can reveal the diversity of cells in tumor tissues. Through single-cell RNA sequencing, it is possible to reveal the cellular heterogeneity in the tumor microenvironment, including different subpopulations of tumor cells, immune cells, stromal cells, etc.; analyze the gene expression patterns of individual cells, identify key molecular markers and signaling pathways; understand the interaction between tumor cells and surrounding cells, including signal transmission and mutual influence between cells; identify key driving factors for tumor progression and metastasis, and provide potential targets for targeted therapy, thereby providing a molecular basis for the formulation of individualized treatment strategies, improving treatment effects and reducing side effects.

[0003] Lung cancer is the leading cause of cancer-related deaths worldwide. In the tumor microenvironment (TME), immune cells play a key role in the occurrence, development and metastasis of tumors. In the tumor state, immune cells have impaired immune function, also known as "immune cell exhaustion", which is usually closely related to the uncontrolled proliferation of malignant tumors. In the lung microenvironment, the state of immune tolerance and suppression provides convenient conditions for the growth and spread of tumors.

[0004] NK cells are located in large numbers in the lungs and play a vital role in lung tumor monitoring and anti-tumor immunity. However, the immunosuppressive microenvironment of the lungs poses a major challenge to the functional characteristics of NK cells, which can lead to their functional decline, exhaustion, and impaired anti-tumor activity. Therefore, in-depth research on immune cells in the lung cancer microenvironment, especially the function and exhaustion mechanism of NK cells, is of great significance for developing new treatment strategies and improving the survival rate of lung cancer patients.

[0005] After searching, the application scheme of Chinese patent application number CN201410205571.9 discloses a high-throughput transcriptome sequencing data quality control method based on multi-core CPU hardware. It includes: using multi-core CPU to process high-throughput transcriptome sequencing data in parallel to obtain data with low sequencing quality sequences removed; using multi-core CPU to predict and remove rRNA sequences in the data with low sequencing quality sequences removed, and to perform qualitative identification of contaminated sequences; and performing statistics and evaluation on sequence alignment results. The high-throughput transcriptome sequencing data quality control method in the above patent has the following shortcomings: the data segmentation and allocation process is relatively complicated, which may lead to uneven allocation of computing resources, thereby affecting the overall processing efficiency; using hidden Markov model to predict and remove rRNA sequences may result in prediction errors, resulting in some non-rRNA sequences being mistakenly predicted as rRNA and removed. Summary of the invention

[0006] The purpose of the present invention is to solve the shortcomings of the prior art and propose an experimental method based on exhaustion of NK cells as a target for lung cancer treatment.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] An experimental method based on exhaustion of NK cells as a target for lung cancer treatment, comprising:

[0009] S1: Data processing, selection of normal lung tissue and lung cancer samples, sequencing data quality control and cell filtering, standardization, dimensionality reduction, clustering, selection of highly variable genes, principal component analysis using the RunPCA function in the Seurat package, and clustering of the first 20 principal components obtained by dimensionality reduction using the FindNeighbors and FindClusters functions in the Seurat package;

[0010] S2: Identification of cell subpopulation characteristic genes and cell type annotation: The FindAllMarkers function was used to obtain the characteristic genes of each subpopulation, and then the SingleR package was used for cell type annotation. Combined with the SingleR package for cell type annotation, matching was performed based on the database celldex:HumanPrimaryCellAtlasData to obtain the predicted label of each cell;

[0011] S3: NK cell re-clustering, NK cells were extracted for separate clustering, the resolution was controlled at 0.1, and a total of 7 cell populations were obtained;

[0012] S4: NK cell function annotation;

[0013] S5: Differential gene expression analysis: Analysis of the high expression of RGS1 in exhausted NK cells and drawing of volcano plots;

[0014] S6: Single-cell analysis of differential expression of activation and inhibitory receptors in RGS1+ NK cells;

[0015] S7: Pseudo-time series analysis was used to infer the dynamic changes of cell states, and NK cells were analyzed to infer pseudo-time trajectories; group-specific marker genes were selected, cells were subjected to dimensionality reduction, cells were sorted in pseudo-time, and the dynamic changes of cell states were inferred;

[0016] S8: Cell culture and siRNA transfection;

[0017] S9: qRT-PCR detection of activation receptor expression in NK92mi cells after RGS1 knockdown;

[0018] S10: lactate dehydrogenase release assay;

[0019] S11: IFN-γ content in cell supernatant was measured by enzyme-linked immunosorbent assay;

[0020] S12: Virtual drug screening using RGS1 as a target.

[0021] Preferably: the NK cell function annotation of S4 is as follows:

[0022] Functional enrichment analysis was performed using Gene Set Enrichment Analysis (GSEA) to complete KEGG pathway enrichment analysis, and the c2.cp.kegg gene set of the Molecular Signatures Database was used for annotation; the functional differences of different cell subpopulations in specific biological processes were evaluated using the enrichment score (ES) indicator.

[0023] Preferably: in S5, the differential gene expression analysis comprises:

[0024] S51: Analysis of RGS1 overexpression. The overexpression of RGS1 in exhausted NK cells was analyzed based on p_val_adj and avg_log2FC. A volcano plot was drawn to show the significant differences and expression changes of genes.

[0025] Preferably, in S6, the single cell analysis of the differences in the expression of activation and inhibitory receptors of RGS1+NK cells is as follows:

[0026] S61: Differential expression of activating and inhibitory receptors, differential analysis of inhibitory and activating receptors: The FeaturePlot function in the Seurat software package was used to perform a detailed analysis of the differential expression of activating receptors NCR3, FCGR3A and inhibitory receptors TIGIT, PDCD1, and LAG3 in NK cell clusters. The point size was set to 1.0 to optimize the visualization effect, and label display was enabled to directly identify the expression differences of these key genes in the figure.

[0027] Preferably, in S7, the dynamic changes of the cell state inferred by pseudo-time series analysis include:

[0028] S71: NK cells were analyzed using Monocle 2 software to infer pseudo-time trajectories;

[0029] S72: Select group-specific marker genes using the detectGenes function;

[0030] S73: Use reduceDimension function to reduce the dimension of cells;

[0031] S74: Use the orderCells function to perform pseudo-time sorting of cells by constructing a tree structure based on cell expression patterns to infer the dynamic change trajectory of cell states;

[0032] S75: Visualizing the dynamics of expressions along a trajectory using the BEAM function.

[0033] Preferably: in said S8, cell culture and siRNA transfection include:

[0034] NK92mi cells were cultured in MEMα and 12.5% ​​fetal bovine serum, 12.5% ​​horse serum at 37°C and 5% CO2. si-RGS1: 5′-CUGUAAAGCAGAAGAGAUA(dT)(dT)-3′ and 3′-UAUCUCUUCUGCUUUACAG(dT)(dT)-5′ were transfected into NK92mi cells, and the knockdown efficiency was detected by qRT-PCR.

[0035] The beneficial effects of the present invention are:

[0036] 1. Through bioinformatics analysis and experimental verification, the present invention reveals the key role of RGS1 in NK cell exhaustion in lung cancer, providing a new and potential molecular target for the treatment of lung cancer.

[0037] 2. This invention uses integrated analysis of single-cell RNA sequencing data combined with in vitro experiments to provide a comprehensive method to understand the dynamic changes and functional status of NK cells in lung cancer.

[0038] 3. The present invention not only identifies new therapeutic targets, but also conducts virtual drug screening based on RGS1, which helps to quickly discover and develop new therapeutic drugs.

[0039] 4. The present invention is verified by in vitro experiments, and the present invention demonstrates that knocking down RGS1 can significantly enhance the killing effect of NK92mi cells on lung cancer cells, which indicates that the method of the present invention may provide more effective therapeutic effects in clinical practice.

[0040] 5. The strategy of the present invention is not only applicable to lung cancer, but may also be extended to other types of cancer, in which NK cell exhaustion is a common problem, thus having broader application prospects.

[0041] 6. The research results of the present invention are based on solid scientific data and experimental evidence, including gene expression analysis, functional validation experiments and drug screening data, providing a reliable basis for further clinical application. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 This is a schematic diagram showing a positive correlation between the KEGG_NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY pathway in the NK cell 0 subset;

[0043] Figure 2 This is a schematic diagram showing a negative correlation between the KEGG_NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY pathway in NK cell subsets 3 and 4;

[0044] Figure 3 Schematic diagram showing negative correlation of KEGG_OXIDATIVE_PHOSPHORYLATION pathway in NK cell 2 subsets;

[0045] Figure 4 The volcano plot shows that RGS1 is highly expressed in exhausted NK cells;

[0046] Figure 5 Schematic diagram showing that the expression level of activation receptors in RGS1+ depleted NK cells is significantly decreased;

[0047] Figure 6 Schematic diagram of the corresponding increase in inhibitory receptors in RGS1+ depleted NK cells;

[0048] Figure 7 Schematic diagram of the differentiation trajectory of NK cells analyzed by Monocle2;

[0049] Figure 8Schematic diagram of RGS1 mRNA expression level in NK92mi cells transfected with siRNA;

[0050] Fig. 9 Schematic diagram of increased expression of activated receptors after RGS1 knockdown;

[0051] Fig.10 This is a schematic diagram showing that the IFN-γ content in the cell supernatant significantly increased after RGS1 knockdown;

[0052] Fig.11 Schematic diagram showing that knockdown of RGS1 enhances the cytotoxic activity of NK92mi cells against target cells. DETAILED DESCRIPTION

[0053] The technical solution of the present invention is further described in detail below in conjunction with specific implementation methods.

[0054] Embodiment 1:

[0055] An experimental method based on exhaustion of NK cells as a target for lung cancer treatment, comprising:

[0056] S1: Data processing, selection of normal lung tissue and lung cancer samples, sequencing data quality control and cell filtering, standardization, dimensionality reduction, clustering, selection of highly variable genes, principal component analysis using the RunPCA function in the Seurat package, and clustering of the first 20 principal components obtained by dimensionality reduction using the FindNeighbors and FindClusters functions in the Seurat package;

[0057] S2: Identification of cell subpopulation characteristic genes and cell type annotation: The FindAllMarkers function was used to obtain the characteristic genes of each subpopulation, and then the SingleR package was used for cell type annotation. Through SingleR analysis, each cell was matched with the most similar reference cell type to obtain the predicted label of each cell;

[0058] S3: NK cell re-clustering, NK cells were extracted for separate clustering, the resolution was controlled at 0.1, and a total of 7 cell populations were obtained;

[0059] S4: NK cell function annotation;

[0060] S5: Differential gene expression analysis: Analysis of the high expression of RGS1 in exhausted NK cells and drawing of volcano plots;

[0061] S6: Single-cell analysis of the differences in the expression of activation and inhibitory receptors in RGS1+NK cells: The FeaturePlot function in the Seurat software package was used for analysis; inhibitory receptors TIGIT, PDCD1, and LAG3 were selected; activation receptors NCR3 and FCGR3A;

[0062] S7: Pseudo-time series analysis was used to infer the dynamic changes of cell states, and NK cells were analyzed to infer pseudo-time trajectories; group-specific marker genes were selected, cells were subjected to dimensionality reduction, cells were sorted in pseudo-time, and the dynamic changes of cell states were inferred;

[0063] S8: Cell culture and siRNA transfection;

[0064] S9: qRT-PCR detection of activation receptor expression in NK92mi cells after RGS1 knockdown;

[0065] S10: lactate dehydrogenase (LDH) release experiment;

[0066] S11: Enzyme-linked immunosorbent assay (ELISA) was used to measure the IFN-γ content in the cell supernatant;

[0067] S12: Virtual drug screening using RGS1 as a target.

[0068] Wherein, in said S1, cluster analysis includes:

[0069] Clustering method selection: Use the FindNeighbors and FindClusters functions in the Seurat package to cluster the first 20 principal components (dims=1:20) after dimensionality reduction, and control the resolution at 0.5 (resolution=0.5);

[0070] Graphical display,The clustering results are displayed through UMAP graphs to observe the spatial distribution of different cell subpopulations.

[0071] Among them, in S2, the cell type annotation includes:

[0072] S21: Identification of characteristic genes. The FindAllMarkers function was used to obtain the characteristic genes of each subgroup. The SingleR package was used to annotate cell types and match them based on the celldex:HumanPrimaryCellAtlasData database.

[0073] Among them, the NK cell function annotation of S4 is as follows:

[0074] Functional enrichment analysis was performed using Gene Set Enrichment Analysis (GSEA) to complete KEGG pathway enrichment analysis, and the c2.cp.kegg gene set of the Molecular Signatures Database was used for annotation; the functional differences of different cell subpopulations in specific biological processes were evaluated using the enrichment score (ES) indicator.

[0075] Wherein, in said S5, the differential gene expression analysis comprises:

[0076] S51: RGS1 high expression analysis, based on p_val_adj and avg_log2FC analysis of RGS1 high expression in exhausted NK cells, drawing volcano plots to show the significant differences and expression changes of genes;

[0077] In S6, single cell analysis of the differences in the expression of activation and inhibitory receptors of RGS1+NK cells was performed as follows:

[0078] S61: Differential expression of activating and inhibitory receptors, differential analysis of inhibitory and activating receptors: The FeaturePlot function in the Seurat software package was used to perform a detailed analysis of the differential expression of activating receptors NCR3, FCGR3A and inhibitory receptors TIGIT, PDCD1, and LAG3 in NK cell clusters. The point size was set to 1.0 to optimize the visualization effect, and label display was enabled to directly identify the expression differences of these key genes in the figure.

[0079] Among them, in S7, the dynamic changes of the cell state inferred by pseudo-time series analysis include:

[0080] S71: NK cells were analyzed using Monocle 2 software to infer pseudo-time trajectories;

[0081] S72: Select group-specific marker genes using the detectGenes function;

[0082] S73: Use reduceDimension function to reduce the dimension of cells;

[0083] S74: Use the orderCells function to perform pseudo-time sorting of cells by constructing a tree structure based on cell expression patterns to infer the dynamic change trajectory of cell states;

[0084] S75: Visualizing the dynamics of expressions along a trajectory using the BEAM function.

[0085] Wherein, in said S8, cell culture and siRNA transfection include:

[0086] NK92mi cells were cultured in MEMα and 12.5% ​​fetal bovine serum, 12.5% ​​horse serum at 37°C and 5% CO2. si-RGS1: 5′-CUGUAAAGCAGAAGAGAUA(dT)(dT)-3′ and 3′-UAUCUCUUCUGCUUUACAG(dT)(dT)-5′ were transfected into NK92mi cells, and the knockdown efficiency was detected by qRT-PCR.

[0087] Among them, in S9, qRT-PCR detection of NK92mi cell activation receptor expression after RGS1 knockdown includes:

[0088] Total RNA was extracted from cells using the Sparkeasy Tissue / Cell RNA Rapid Extraction Kit; PrimeScript with gDNAEraser was used TM The mRNA was reverse transcribed into cDNA using the RT kit; TBGreenExTaqII premix and cDNA were added to a 10 μL system for real-time quantitative PCR; the RGS1 primer sequences were: ForwardPrimer: TCTTCTCTGCTAACCCAAAGGA, ReversePrime: TGCTTTACAGGGCAAAAGATCAG; the primer sequences for activating receptor NKG2D were: ForwardPrimer: ACTCAGAGATGAGCAAATGCC, ReversePrime: CAGGTTGACTGGTAGTTAGTGC, the mRNA expression level was normalized to the GAPDH expression level, and the relative quantification was expressed by 2-△△CT; all experiments were repeated 3 times.

[0089] Wherein, in said S10, the lactate dehydrogenase release experiment comprises:

[0090] S101: Use the lactate dehydrogenase (LDH) cytotoxicity detection kit (Shanghai Biotech Co., Ltd.) to prepare the number of cells required for the experiment according to the effector-target ratio of 3:1;

[0091] S102: The density of target cell A549 cells was adjusted to 1×10^5 cells / mL, 100 μL was taken and added to a 96-well culture plate, and effector cells NK92mi cells were added according to different effector-target ratios, and the volume was fixed to 200 μL; 100 μL of target cells and culture medium were added to the target cell natural release well, and 100 μL of LDH release reagent provided by the kit diluted 10 times with PBS and target cells were added to the target cell maximum release well, and placed in the incubator for 4 hours;

[0092] S103: Then centrifuge the 96-well culture plate at 250g for 5 minutes, take 120μl of the supernatant from each well, add it to the corresponding well of a new 96-well plate, add 60μl of LDH working solution, incubate at room temperature (about 25°C) in the dark for 30 minutes, and measure the absorbance value at 490nm on an enzyme reader;

[0093] S104: The effector cytotoxicity (%) was calculated according to the following formula = [reaction pore D (490 nm) - natural release pore D (490 nm) / maximum release pore D (490 nm) - natural release pore D (490 nm)] × 100%; all experiments were repeated 3 times.

[0094] Wherein, in said S11, the ELISA method for determining the IFN-γ content in the cell supernatant comprises:

[0095] S111: Prepare the number of cells required for the experiment according to the effector-target ratio of 3:1; adjust the density of target cell A549 cells to 1×10^5 cells / mL, take 100μL and add it to a 96-well culture plate, add effector cell NK92mi cells according to different effector-target ratios, make the volume to 200μL, take the cell culture supernatant and centrifuge at 1000xg for 15 minutes; use the Human Interferonγ(IFN-γ) ELISA Kit (CSB-E04577h, Huamei Bio) to determine the IFN-γ content in the cell supernatant;

[0096] S112: Move various reagents to room temperature (18-25°C) to equilibrate for at least 30 minutes, prepare reagents according to the above method, and set aside;

[0097] S113: Sample addition: set up standard wells and test sample wells respectively; add 100ul of standard or test sample to each well, shake gently to mix, cover with plate sticker, and incubate at 37℃ for 2 hours;

[0098] S114: Discard the liquid, spin dry, and do not wash; add 100 μL of biotin-labeled antibody working solution to each well, cover with a new plate sticker, and incubate at 37°C for 1 hour;

[0099] S115: Discard the liquid in the wells, spin dry, wash the plate 3 times; soak for 2 minutes each time, 200 μL / well, spin dry;

[0100] S116: Add 100ul of horseradish peroxidase-labeled avidin working solution to each well, cover with a new plate sticker, and incubate at 37°C for 1 hour;

[0101] S117: Discard the liquid in the wells, spin dry, wash the plate 5 times; soak for 2 minutes each time, 200ul per well, spin dry;

[0102] S118: Add 90ul of substrate solution to each well in sequence and color at 37℃ in the dark for 15-30 minutes; add 50ul of stop solution to each well in sequence to terminate the reaction;

[0103] S119: Within 5 minutes after the reaction is terminated, the optical density (OD value) of each well is measured in sequence at a wavelength of 450 nm using a microplate reader.

[0104] Wherein, in S12, taking RGS1 as a target, performing drug virtual screening includes:

[0105] S121: Preparation of RGSS1 protein and compound library: The protein crystal structure was first pre-processed to add hydrogen to all atoms lacking hydrogen in the structure; the protonation state of the amino acid residues was adjusted to optimize the hydrogen bond distribution of the protein; water molecules were removed according to the hydrogen bond pattern, and the energy was minimized using the OPLS3e force field; compounds with clear or potential immunopotentiating effects were obtained from the natural product library in the MCE immunopotentiator library database, and geometry was optimized using the LigPrep module; Epik was used to generate the ionized state and tautomers of the ligand at pH = 7.0 ± 2.0; for stereoisomers, a maximum of 32 isomers were obtained for each natural product;

[0106] S122: RGS1 binding pocket prediction: Use DoGSiteScorer to predict the binding site of RGS1; apply Gaussian difference filter to predict chemical drugability; predict pocket drugability score values ​​between 0 and 1, the higher the score, the more druggable the pocket;

[0107] S123: Screening of compound library: Based on the identified binding sites, natural compounds in the MCE database were subjected to structure-based virtual screening using the Glide module; the prepared database was pre-filtered using the LigFilter module to remove molecules that did not meet the drug-likeness characteristics; the obtained compounds were preliminarily screened using high-throughput virtual screening (HTVS) and ranked using the Glide score (GScore); standard precision (SP) and extra precision (XP) were used to screen possible best candidates; after HTVS docking, the poses of 20% of the best compounds were output for SP docking, and then 20% of the compounds were screened for XP docking.

[0108] In the present invention, the NK cell subsets are annotated as follows:

[0109] Through Gene Set Enrichment Analysis (GSEA), the functional exploration of NK cell subsets was carried out in depth. The analysis results showed that the KEGG_NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY pathway was positively correlated in NK cell subset 0, but negatively correlated in subsets 3 and 4. Specifically, in the positively correlated subgroups, the enrichment score (ES) of this pathway was significantly higher, indicating that these cell subsets may have stronger cytotoxic functions, which is consistent with the enhanced cytotoxic activity mediated by NK cells. In contrast, in the negatively correlated subgroups, the ES of this pathway was lower, suggesting that the cytotoxic function of these cell subsets may be inhibited and their efficiency is reduced when performing their natural killing function. This difference in pathway activity may reflect the different functional states of NK cells in the tumor microenvironment, where the positively correlated subgroups may play a more active role in the anti-tumor immune response, while the negatively correlated subgroups may be associated with immunosuppression or tolerance and present an exhausted phenotype. Therefore, cell cluster 0 was defined as the “Active nk” subpopulation, and cell clusters 3 and 4 were defined as the “exhausted nk” subpopulations.

[0110] It was also observed that in NK cell subset 2, the KEGG_OXIDATIVE_PHOSPHORYLATION pathway showed a significant negative correlation with the gene set. In the resting state, NK cells mainly rely on oxidative phosphorylation to produce energy, while in the activated state, glycolysis and oxidative phosphorylation levels are significantly enhanced to support their effector functions. Oxidative phosphorylation not only provides energy for NK cells, but also participates in regulating the activation, survival and function of NK cells. Blocking the oxidative phosphorylation pathway can significantly inhibit the activation of NK cells and their ability to secrete cytokines and kill target cells. Therefore, cell cluster 2 is defined as the "exhausted nk_lowoxphos" subset.

[0111] Through Gene Set Enrichment Analysis (GSEA), the functional status of NK cell subsets was deeply analyzed. The analysis results revealed that the pathways related to cell proliferation in subgroup 6: KEGG_CELL_CYCLE, KEGG_OOCYTE_MEIOSIS, and KEGG_DNA_REPLICATION showed significant positive correlation. In particular, cell cycle genes, including PCNA and MKI67, were observed to have significantly increased expression levels in the subgroup. PCNA, as a key protein involved in DNA replication and repair, and MKI67, as a marker of cell proliferation, their high expression indicate that these NK cell subsets may be in a highly proliferative state. This finding prompted the definition of subgroup 6 as a "proliferating nk" subgroup. At the same time, it was found that through DEG combined with pseudo-time series analysis, NK cell subgroup 1 was an excessive NK cell cluster, so subgroup 1 was defined as a "transitonal nk" subgroup. Therefore, based on gene set enrichment analysis and differential gene expression in different NK cell subsets, NK cell subsets were re-annotated as “Active nk”, “exhausted nk”, “exhausted nk_low oxphos”, “proliferating nk”, and “transitonal nk”.

[0112] The present invention also includes a virtual drug screening method based on RGS1 as a target:

[0113] ① Preparation of RGSS1 protein and compound library

[0114] The crystal structure of RGS1 (PDB ID: 2BV1) was downloaded from the RCSB PDB database and imported into the ProteinPreparation module for protein preparation. The protein crystal structure was first preprocessed to add hydrogen to all atoms that lacked hydrogen in the structure. The protonation state of the amino acid residues was then adjusted to optimize the hydrogen bond distribution of the protein. Water molecules were removed based on their hydrogen bonding patterns, and energy was minimized using the OPLS3e force field. Compounds with clear or potential immunoenhancing effects were obtained from the natural product library in the MCE immunoenhancer library database and geometry optimized using the LigPrep module. Epik was used to generate ionized states and tautomers of the ligand at pH = 7.0 ± 2.0. For stereoisomers, up to 32 isomers were obtained for each natural product.

[0115] ②Prediction of RGS1 binding pocket

[0116] To predict the drugability of the docking pocket, the binding site prediction of RGS1 was performed using DoGSiteScorer (DGSS). It applies a Difference of Gaussian (DoG) filter to predict chemical drugability. The predicted pocket drugability score value is between 0 and 1, and the higher the score, the more druggable the pocket.

[0117] ③Screening of compound library

[0118] Based on the identified binding sites, natural compounds in the MCE database were subjected to structure-based virtual screening using the Glide module. The prepared database was pre-filtered using the LigFilter module to remove molecules that did not meet the drug-likeness properties. The obtained compounds were first preliminarily screened using high-throughput virtual screening (HTVS) and ranked using the Glide score (GScore). Standard precision (SP) and extra precision (XP) were then used to screen possible best candidates. After HTVS docking, the poses of 20% of the best compounds were output for SP docking, and then 20% of the compounds were screened for XP docking.

[0119] In summary, drug virtual screening is based on predicted binding sites. The top 20% of compounds screened are filtered through HTVS and docked using SP. Then the top 20% of compounds are screened for XP docking. Finally, 5 drugs are retained based on the final XP docking results. As follows:

[0120]

[0121] Experimental results analysis:

[0122] 1. Single-cell transcriptome data processing and NK cell subset annotation

[0123] By integrating 20 single-cell transcriptome data and undergoing strict quality control, 130,984 high-quality cells were obtained. After removing batch effects through data integration and performing dimensionality reduction analysis using PCA, 10 major cell subsets were obtained through graphical clustering using UMAP and annotation using the SingleR package. NK cells were extracted for further clustering, resulting in 7 NK cell subsets. Through Gene Set Enrichment Analysis (GSEA), we conducted an in-depth functional exploration of NK cell subsets. Based on gene set enrichment analysis and differential gene expression in different NK cell subsets, NK cell subsets were re-annotated as "Active nk", "exhausted nk", "exhausted nk_low oxphos", "proliferating nk", and "transitonal nk".

[0124] 2. RGS1+ NK cells exhibit an exhausted phenotype

[0125] Pseudo-time series analysis showed that the expression level of activating receptors in RGS1+NK cell clusters decreased significantly when exhausted NK cells were at the end of differentiation trajectory, while the expression level of inhibitory receptors increased accordingly. This imbalance in the expression levels of activating and inhibitory receptors may have an important impact on the cytotoxicity and immune surveillance functions of NK cells.

[0126] 3. Knockdown of RGS1 expression in NK cells helps reverse NK cell exhaustion. The expression level of RGS1 mRNA in NK92mi cells after si transfection was detected by qRT-PCR. It was found that the RGS1 mRNA level in the si-RGS1 group was downregulated by about 50% compared with si-NC. The effect of RGS1 on the activation receptors of NK92mi cells was detected by qRT-PCR. The results showed that the activation receptors in the si-RGS1 group were significantly higher than those in the control group. The results of ELISA showed that the IFN-γ content in the cell supernatant of the si-RGS1 group was significantly increased compared with the control group. After knockdown of RGS1 in NK92mi cells, they were co-cultured with A549 cells, and the effector-target ratio of 3:1 was set to detect the release of LDH after 4 hours, and the effector cytotoxicity was calculated. The results showed that knockdown of RGS1 enhanced the killing activity of NK92mi cells against target cells.

[0127] 4. Muramyl dipeptide and other five drugs can be used as potential immune activators. Based on the predicted binding sites, the top 20% of the compounds screened out were filtered through HTVS and docked using SP. Then the top 20% of the compounds were screened for XP docking. Finally, based on the XP docking results, the five drugs with the strongest binding affinity were selected: Muramyl dipeptide, STING agonist-3, C12-TLRa, diABZI STING agonist-1, RS 09 (TFA) as potential immune activators for reversing NK cell exhaustion.

[0128] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. An experimental method based on exhaustion of NK cells as a target for lung cancer treatment, characterized in that: include: S1: Data processing, selection of normal lung tissue and lung cancer samples, sequencing data quality control and cell filtering, standardization, dimensionality reduction, clustering, selection of highly variable genes, principal component analysis using the RunPCA function in the Seurat package, and clustering of the first 20 principal components obtained by dimensionality reduction using the FindNeighbors and FindClusters functions in the Seurat package; S2: Identification of cell subpopulation characteristic genes and cell type annotation: The FindAllMarkers function was used to obtain the characteristic genes of each subpopulation, and then the SingleR package was used for cell type annotation. Combined with the SingleR package for cell type annotation, matching was performed based on the database celldex:HumanPrimaryCellAtlasData to obtain the predicted label of each cell; S3: NK cell re-clustering, NK cells were extracted for separate clustering, the resolution was controlled at 0.1, and a total of 7 cell populations were obtained; S4: NK cell function annotation; S5: Differential gene expression analysis: Analysis of the high expression of RGS1 in exhausted NK cells and drawing of volcano plots; S6: Single-cell analysis of differential expression of activation and inhibitory receptors in RGS1+ NK cells; S7: Pseudo-time series analysis was used to infer the dynamic changes of cell states, and NK cells were analyzed to infer pseudo-time trajectories; group-specific marker genes were selected, cells were subjected to dimensionality reduction, cells were sorted in pseudo-time, and the dynamic changes of cell states were inferred; S8: Cell culture and siRNA transfection; S9: qRT-PCR detection of activation receptor expression in NK92mi cells after RGS1 knockdown; S10: lactate dehydrogenase release assay; S11: IFN-γ content in cell supernatant was measured by enzyme-linked immunosorbent assay; S12: Virtual drug screening using RGS1 as a target.

2. According to claim 1, an experimental method based on exhausted NK cells as a target for lung cancer treatment, characterized in that: The S4 NK cell function annotations are as follows: Functional enrichment analysis was performed using Gene Set Enrichment Analysis to complete KEGG pathway enrichment analysis, and the c2.cp.kegg gene set of the Molecular Signatures Database was used for annotation; the functional differences of different cell subpopulations in specific biological processes were evaluated using the enrichment score indicator.

3. The experimental method according to claim 2, wherein the exhausted NK cells are used as a target for lung cancer treatment, In S5, the differential gene expression analysis includes: S51: Analysis of RGS1 overexpression. The overexpression of RGS1 in exhausted NK cells was analyzed based on p_val_adj and avg_log2FC. A volcano plot was drawn to show the significant differences and expression changes of genes.

4. The experimental method according to claim 3 based on exhaustion of NK cells as a target for lung cancer treatment, characterized in that: In S6, single cell analysis of the differences in the expression of activation and inhibitory receptors of RGS1+NK cells was performed as follows: S61: Differential expression of activating and inhibitory receptors, differential analysis of inhibitory and activating receptors: The FeaturePlot function in the Seurat software package was used to perform a detailed analysis of the differential expression of activating receptors NCR3, FCGR3A and inhibitory receptors TIGIT, PDCD1, and LAG3 in NK cell clusters. The point size was set to 1.0 to optimize the visualization effect, and label display was enabled to directly identify the expression differences of these key genes in the figure.

5. The experimental method according to claim 4, wherein the exhausted NK cells are used as a target for lung cancer treatment, In S7, the dynamic changes of cell states inferred by pseudo-time series analysis include: S71: NK cells were analyzed using Monocle 2 software to infer pseudo-time trajectories; S72: Select group-specific marker genes using the detectGenes function; S73: Use reduceDimension function to reduce the dimension of cells; S74: Use the orderCells function to perform pseudo-time sorting of cells by constructing a tree structure based on cell expression patterns to infer the dynamic change trajectory of cell states; S75: Visualizing the dynamics of expressions along a trajectory using the BEAM function.

6. The experimental method according to claim 5 based on exhaustion of NK cells as a target for lung cancer treatment, characterized in that: In S8, cell culture and siRNA transfection include: NK92mi cells were cultured in MEMα and 12.5% ​​fetal bovine serum, 12.5% ​​horse serum at 37°C and 5% CO2. si-RGS1: 5′-CUGUAAAGCAGAAGAGAUA(dT)(dT)-3′ and 3′-UAUCUCUUCUGCUUUACAG(dT)(dT)-5′ were transfected into NK92mi cells, and the knockdown efficiency was detected by qRT-PCR.

7. The experimental method according to claim 6, wherein the exhausted NK cells are used as a target for lung cancer treatment, In S9, qRT-PCR detection of NK92mi cell activation receptor expression after RGS1 knockdown includes: Total RNA was extracted from cells using the Sparkeasy Tissue / Cell RNA Rapid Extraction Kit; PrimeScript with gDNAEraser was used TM The mRNA was reverse transcribed into cDNA using the RT kit; TBGreenExTaqII premix and cDNA were added to a 10 μL system for real-time quantitative PCR; the RGS1 primer sequences were: ForwardPrimer: TCTTCTCTGCTAACCCAAAGGA, ReversePrime: TGCTTTACAGGGCAAAAGATCAG; the primer sequences for activating receptor NKG2D were: ForwardPrimer: ACTCAGAGATGAGCAAATGCC, ReversePrime: CAGGTTGACTGGTAGTTAGTGC, the mRNA expression level was normalized to the GAPDH expression level, and the relative quantification was expressed by 2-△△CT; all experiments were repeated 3 times.

8. The experimental method according to claim 7, wherein the exhausted NK cells are used as a target for lung cancer treatment, In S10, the lactate dehydrogenase release experiment includes: S101: Use the lactate dehydrogenase cytotoxicity assay kit to prepare the number of cells required for the experiment based on a 3:1 effector-target ratio; S102: The density of target cell A549 cells was adjusted to 1×10^5 cells / mL, 100 μL was taken and added to a 96-well culture plate, and effector cells NK92mi cells were added according to different effector-target ratios, and the volume was fixed to 200 μL; 100 μL of target cells and culture medium were added to the target cell natural release well, and 100 μL of LDH release reagent provided by the kit diluted 10 times with PBS and target cells were added to the target cell maximum release well, and placed in the incubator for 4 hours; S103: Then centrifuge the 96-well culture plate at 250g for 5 minutes, take 120μl of the supernatant from each well, add it to the corresponding well of a new 96-well plate, add 60μl of LDH working solution, incubate at room temperature in the dark for 30 minutes, and measure the absorbance value at 490nm on an enzyme reader; S104: Calculate as follows: Effector cytotoxicity (%) = [reaction well D - natural release well D / maximum release well D - natural release well D] × 100%; all experiments were repeated 3 times.

9. The experimental method based on exhaustion of NK cells as a target for lung cancer treatment according to claim 8, characterized in that: In S11, the ELISA method for determining the IFN-γ content in the cell supernatant comprises: S111: Prepare the number of cells required for the experiment according to the effector-target ratio of 3:1; adjust the density of target cell A549 cells to 1×10^5 cells / mL, take 100μL and add it to a 96-well culture plate, add effector cell NK92mi cells according to different effector-target ratios, make the volume to 200μL, take the cell culture supernatant and centrifuge at 1000xg for 15 minutes; use the Human InterferonγELISA Kit to determine the IFN-γ content in the cell supernatant; S112: Move various reagents to room temperature, equilibrate for at least 30 minutes, prepare reagents and set aside; S113: Sample addition: set up standard wells and test sample wells respectively; add 100ul of standard or test sample to each well, shake gently to mix, cover with plate sticker, and incubate at 37℃ for 2 hours; S114: Discard the liquid, spin dry, and do not wash; add 100 μL of biotin-labeled antibody working solution to each well, cover with a new plate sticker, and incubate at 37°C for 1 hour; S115: Discard the liquid in the wells, spin dry, wash the plate 3 times; soak for 2 minutes each time, 200 μL / well, spin dry; S116: Add 100ul of horseradish peroxidase-labeled avidin working solution to each well, cover with a new plate sticker, and incubate at 37°C for 1 hour; S117: Discard the liquid in the wells, spin dry, wash the plate 5 times; soak for 2 minutes each time, 200ul per well, spin dry; S118: Add 90ul of substrate solution to each well in sequence and color at 37℃ in the dark for 15-30 minutes; add 50ul of stop solution to each well in sequence to terminate the reaction; S119: Within 5 minutes after the reaction is terminated, the optical density of each well is measured in sequence at a wavelength of 450 nm using an enzyme reader.

10. An experimental method based on exhaustion of NK cells as a target for lung cancer treatment according to claim 9, characterized in that: In S12, using RGS1 as a target, virtual drug screening includes: S121: Preparation of RGSS1 protein and compound library: The protein crystal structure was first pre-processed to add hydrogen to all atoms lacking hydrogen in the structure; the protonation state of the amino acid residues was adjusted to optimize the hydrogen bond distribution of the protein; water molecules were removed according to the hydrogen bond pattern, and the energy was minimized using the OPLS3e force field; compounds with clear or potential immunopotentiating effects were obtained from the natural product library in the MCE immunopotentiator library database, and geometry was optimized using the LigPrep module; Epik was used to generate the ionized state and tautomers of the ligand at pH = 7.0 ± 2.0; for stereoisomers, a maximum of 32 isomers were obtained for each natural product; S122: RGS1 binding pocket prediction: Use DoGSiteScorer to predict the binding site of RGS1; apply Gaussian difference filter to predict chemical drugability; predict pocket drugability score values ​​between 0 and 1, the higher the score, the more druggable the pocket; S123: Screening of compound library: Based on the identified binding sites, natural compounds in the MCE database were subjected to structure-based virtual screening using the Glide module; the prepared database was pre-filtered using the LigFilter module to remove molecules that did not meet the drug-likeness properties; the obtained compounds were preliminarily screened using high-throughput virtual screening and ranked using the Glide score; standard precision and ultra-precision were used to screen possible best candidates; after HTVS docking, the poses of 20% of the best compounds were output for SP docking, and then 20% of the compounds were screened for XP docking.

Images