A method of screening for biomarkers of b-cell lymphoma

Abstract

The application discloses a method for screening B-cell lymphoma biomarkers and belongs to the technical field of bioinformatics. The application aims to provide a more accurate and simple method for screening B-cell lymphoma biomarkers. The application provides a method for screening B-cell lymphoma biomarkers, which utilizes scRNA-seq data of Burkitt lymphoma, diffuse large B-cell lymphoma and follicular lymphoma, splenic marginal zone lymphoma and bulk RNA-seq data of Burkitt lymphoma and diffuse large B-cell lymphoma database to screen the relationship between T / NK and myeloid cells all cell subpopulation infiltration ratio and patient survival period. The application provides a theoretical basis for screening lymphoma disease survival analysis.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology, specifically relating to a method for screening biomarkers for B-cell lymphoma. Background Technology

[0002] The core idea behind screening biomarkers for B-cell lymphoma is to compare the molecular differences between tumors and normal tissues (or different subtypes / prognostic groups) using high-throughput technologies, identify molecules (such as gene mutations, differentially expressed genes / proteins, etc.) related to tumor occurrence, progression, or treatment response, and verify their clinical application value. Commonly used methods include genomics (whole exon / targeted sequencing), transcriptomics (RNA-seq, microarrays), proteomics (mass spectrometry, antibody microarrays), and single-cell omics (scRNA-seq), and the workflow includes sample collection, molecular detection, differential analysis, functional validation, and clinical translation.

[0003] Genomics screens for markers by identifying driver mutations (such as MYD88 L265P and BCL2 t(14;18)), but it is easy to confuse driver and passenger mutations, and the detection of copy number variations is limited by sample purity; transcriptomics focuses on differentially expressed genes, but there are problems such as the disconnect between mRNA and protein function and batch effect interference; proteomics directly detects functional molecules, but low-abundance proteins have insufficient sensitivity and ignore post-translational modifications; single-cell omics analyzes intratumoral heterogeneity, but the results are unstable due to data sparsity and the subjectivity of cell annotation.

[0004] The shortcomings of existing methods are present throughout the entire process: At the sample level, tumor heterogeneity (subtype / spatial differences), sample quality (FFPE degradation, low abundance in body fluids), and individual patient differences (treatment history, comorbidities) can easily interfere with the authenticity of signals; At the technical level, the inherent limitations of each method (such as functional blind spots in genomics and regulatory disconnects in the transcriptome) lead to false positives / false negatives; In data analysis, uncorrected multiple testing leads to a proliferation of false positives, differences in bioinformatics tools and algorithms reduce the consistency of results, and there is a lack of pathway integration of molecular mechanisms; In the validation and application stages, in vitro models cannot simulate the tumor microenvironment, and clinical validation often fails due to small sample sizes, lack of independent cohorts, or improper endpoint design (such as emphasizing molecular differences over clinical outcomes). At the same time, it is difficult to balance sensitivity and specificity of markers, and the feasibility of dynamic monitoring is low. Summary of the Invention

[0005] The purpose of this invention is to provide a more accurate and concise method for screening biomarkers for B-cell lymphoma.

[0006] This invention provides a method for screening biomarkers for B-cell lymphoma, the steps of which are as follows: Step 1: Select data: scRNA-seq data of Burkitt lymphoma, diffuse large B-cell lymphoma and follicular lymphoma, and splenic marginal zone lymphoma, and bulk RNA-seq data of Burkitt lymphoma and diffuse large B-cell lymphoma databases; Step 2: Use the R package Seurat to filter cells based on quality control indicators, and use the LogNormalize method to normalize the data with a scaling factor of 10,000 to correct for changes in sequencing depth between units. Step 3: Using the data obtained in Step 2, select 5000 highly variable genes as the core features for downstream dimensionality reduction and clustering, perform principal component analysis, and then generate an Elbow Plot to determine the optimal number of principal components during clustering. Based on the optimal number of principal components, construct a K-nearest neighbor graph and optimize the neighborhood connection weights using Jaccard similarity. Step 4: The cell data obtained in Step 3 is processed by unsupervised cell clustering to separate cells into groups. Cells showing mixed lineage expression patterns are excluded using the marker gene method. Finally, single-cell data are obtained for subsequent merging analysis. Step 5: Add batch data, patient information, and disease type data as annotation columns, and exclude all low-quality cells; Step 6: Use the R package decontX to remove environmental RNA contamination, retaining all cells with contamination levels <0.2; Step 7: Normalize using NormalizeData in Seurat with a scaling factor of 10,000 as the parameter, detect consistent HVG across datasets using FindVariableFeatures, select the top 2,000 features for data integration, scale and center gene expression using ScaleData, and remove batch effects between data batches and patients using the RunHarmony function in the Harmony package. Step 8: Retain the first 30 principal components, use Seurat's FindNeighbors to calculate the Euclidean distance in the Harmony dimensionality reduction space to construct the KNN graph, use FindClusters to perform unsupervised clustering at a resolution of 0.5, further reduce the dimensionality through RunUMAP in Seurat, and visualize it through unified manifold approximation and projection; Step 9: Screen the cell population using marker gene methods and the SingleR tool, take the intersection data, and remove cells that show mixed lineage expression to ensure that confounding cells are excluded; Step 10: Use the R package DoubletFinder to calculate the number of mixed duplexes in each subpopulation, use paramSweep to scan parameters, find the optimal pK value by combining find.pK with the BCreal index, use the clustering results as annotations, and use modelHomotypic to estimate the proportion of homologous duplexes; calculate the expected number of heterologous duplexes based on 14.5% of the total number of cells, and adjust it based on the homologous proportion, and finally identify and delete duplexes in each cluster; Step 11: Extract cells of each type individually for a second round of dimensionality reduction clustering to further characterize the subsets of T / NK cells, myeloid cells, and fibroblasts: Step 12: Based on gene markers and differentially expressed genes, 5 NK cell subsets, 12 CD8 T cell subsets, 7 CD4 T cell subsets, 8 myeloid cell subsets, and 4 fibroblast subsets were identified. Step 13: Use cibersort to deconvolve all bulk RNA-seq datasets to quantify the abundance of cell types. It was found that myeloid cells and T / NK were most significantly correlated. Spearman correlation analysis was performed on the infiltration patterns of cell types, and it was found that there was a significant positive correlation between T / NK and myeloid cells in all 5 independent cohorts. Step 14: Based on the T / NK and myeloid cell data obtained in Step 12, differentially regulated genes at the RNA level were screened using the criteria of logFC ≥ 0.25 and min.pct ≥ 0.1. High-confidence marker genes were further screened using the criteria of ogFC > 0.5 and expression row sum > 1. A mean version of the Pseudobulk reference matrix was constructed based on the high-confidence genes. CIBERSORT deconvolution was run on the DLBCL and BLbulk datasets to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. Step 15: Perform stratified sampling on the T / NK and myeloid cell subpopulations obtained in Step 12, extract the subset counts matrix, and use the full subpopulation Pseudobulk matrix as the reference feature. Similarly, run CIBERSORT deconvolution on the 9 DLBCL and BL bulk datasets to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. Step 16: Take the intersection of the results of Step 14 and Step 15 to obtain the relationship between T / NK and the infiltration ratio of all cell subpopulations of myeloid cells and patient survival.

[0007] To further define the quality control criteria in step 2, the criteria are: only retain cell data in which 200-6000 genes were detected and the proportion of mitochondrial genes was less than 15%.

[0008] Further, the excluded indicators in step 5 are: nFeature_RNA < 200, nCount_RNA > 50,000, and percentage.mt > 15.

[0009] Further, the cell contamination rate was calculated using a binary mixture model via the decontX function. A series of contamination thresholds, ranging from 0.1 to 0.25, were tested, and the real-time visualization of cell retention at different thresholds was achieved using the Shiny interactive tool.

[0010] To further specify, in step 7, the function RunHarmony uses the parameters "theta = 1, lambda = 1.5, sigma = 0.1" to perform two iterations; in step 8, the parameters for visualization are "dims = 1:30, reduction ='harmony'".

[0011] Further specifying, in step 10, the paramSweep scan parameters are PCs=1:40 and sct=FALSE.

[0012] Further specifying, in step 11, for bone marrow cells, the first 29 PCs from 2,000 highly variable genes are used, with a resolution of 0.3; for T / NK cells, the first 40 PCs are extracted from 2,000 highly variable genes, with a resolution of 0.5.

[0013] Further, in step 14, the parameters for running CIBERSORT deconvolution on the DLBCL and BL bulk datasets are perm=0 and QN=TRUE to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample; in step 15, stratified sampling is performed by sampling 9% from large / medium subpopulations, 20% from small subpopulations, and retaining all cells from very small subpopulations, with at least 30 cells retained from each subpopulation.

[0014] This invention provides a method for screening biomarkers for B-cell lymphoma, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the method for screening biomarkers for B-cell lymphoma as described above.

[0015] The present invention provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for screening B-cell lymphoma biomarkers as described above.

[0016] Beneficial effects: To determine the clinical relevance of cell type infiltration, the relationship between the infiltration ratios of all T / NK and myeloid cell subpopulations and patient survival (OS) was examined in all bulk datasets. Ultimately, it was found that the infiltration ratios calculated by both methods appeared in multiple bulk datasets. The cell subpopulations that were significantly associated with survival analysis results were NK_Active_CCL3+ and CD4_Tem_Th1like_MKI67. The enrichment of both subpopulations was significantly associated with better patient survival prognosis. Detailed Implementation

[0017] Example 1. 1. Data preprocessing: (1) Data sources and preprocessing: Data were obtained from multiple public databases and previously published articles, including scRNA-seq and bulk RNA-seq data for various B-cell lymphomas, including Burkitt lymphoma (BL), diffuse large B-cell lymphoma (DLBCL), follicular lymphoma (FL), and splenic marginal zone lymphoma (SMZL). Standardized gene expression data served as the basis for subsequent analyses.

[0018] (2) Data foundation analysis: A. First, basic processing was performed on all single-disease scRNA-seq datasets. After downloading the data, the R package Seurat was used to filter cells based on quality control indicators, retaining only cells with 200-6000 detected genes (nFeature_RNA) and a mitochondrial gene percentage (percent.mt%) below 15%. Then, the data was normalized using the LogNormalize method, with a scaling factor of 10,000 to correct for variations in sequencing depth between units (parameters). Next, 5000 hypervariable genes (HVGs) were selected as the core features for downstream dimensionality reduction and clustering, followed by principal component analysis (PCA). Elbow plots were then generated to determine the optimal number of principal components (PCs) for clustering. A K-nearest neighbor (KNN) graph was constructed using the retained PCs (optimal PCs), and neighborhood connection weights were optimized using Jaccard similarity. Subsequently, unsupervised cell clustering (a well-known method) was performed to group cells, and marker gene methods were used to exclude cells exhibiting mixed lineage expression patterns, ultimately yielding single-cell data for subsequent merging analysis. B. Add batch data, patient information, and disease type data as annotation columns. Exclude all low-quality cells (nFeature_RNA < 200, nCount_RNA > 50,000, and percentage.mt > 15); then use the R package decontX to remove environmental RNA contamination: this package uses a binary mixture model to calculate the cell contamination rate through the decontX function, testing a series of contamination thresholds (0.1-0.25), and using the Shiny interactive tool to achieve real-time visualization of cell retention at different thresholds, finally selecting to retain all cells with a contamination level < 0.2; then continue to use NormalizeData in Seurat to normalize the data with a scaling factor of 10,000. Consistent HVG across datasets was detected using FindVariableFeatures. The top 2,000 features were selected for data ensemble. Gene expression was scaled and centered using ScaleData. Batch effect removal was performed between data batches and patients using the RunHarmony function from the Harmony package, with two iterations performed using parameters "theta = 1, lambda = 1.5, sigma = 0.1". The top 30 principal components were retained to support downstream analysis steps. Euclidean distances in the Harmony dimensionality reduction space were calculated using Seurat's FindNeighbors to construct a KNN graph. Unsupervised clustering was performed using FindClusters at a resolution of 0.5. Dimensionality reduction was further achieved using RunUMAP from Seurat and visualized using the Unified Manifold Approximation and Projection (UMAP) with parameters "dims = 1:30, reduction = 'harmony'". Classic literature markers and SingleR content tools were used to cross-reference and take intersections to remove cells exhibiting mixed lineage expression, ensuring the exclusion of confounding cells. The main cell types identified include B cells (CD79A, CD79B, MS4A1), T / NK cells (CD3D, CD3E, NKG7), bone marrow cells (LYZ, S100A8, S100A9, CD14, FCN1, CLEC4C, LILRA4, PACSIN1), oligodendrocytes (MOG), erythrocytes (HBA1, HBA2, HBB), and fibroblasts (COL1A1, COL3A1, PDPN). C. Subsequently, the R package DoubletFinder was used to calculate the heterologous pairs (individual cells) in each subpopulation. The parameters were scanned using paramSweep (PCs=1:40, sct=FALSE). The optimal pK value was found by combining find.pK with the BCreal index. The proportion of homologous pairs was estimated using modelHomotypic with the clustering results (RNA_snn_res.0.3) as annotation. The expected number of heterologous pairs was calculated at 14.5% of the total number of cells and adjusted in conjunction with the homologous proportion. Finally, pairs in each cluster were identified and removed. A second round of dimensionality reduction clustering was performed separately for each cell type to further characterize the subpopulations of T / NK cells, myeloid cells, and fibroblasts: For myeloid cells, the top 29 PCs from 2,000 highly variable genes (resolution=0.3) were used; for T / NK cells, the top 40 PCs from 2,000 highly variable genes (resolution=0.5) were extracted. Finally, based on classical markers and differentially expressed genes (DEGs) in the literature, we identified 5 NK cell subsets, 12 CD8 T cell subsets, 7 CD4 T cell subsets, 8 myeloid cell subsets, and 4 fibroblast subsets. D. Ultimately, more than 170,000 cells were obtained, and 39 cell subpopulation types were identified.

[0019] 2. Cibersort immune infiltration analysis deconvolution: A. Cibersort was used to deconvolve all bulk RNA-seq datasets to quantify the abundance of nine major cell types. The study found that myeloid cells and T / NK cells showed the most significant correlation. Spearman correlation analysis of the infiltration patterns of the nine major cell types revealed a significant positive correlation between T / NK cells and myeloid cells in all five independent cohorts. Further analysis was then conducted on the subgroups of T / NK cells and myeloid cells. B. Two methods were used to construct CIBERSORT reference matrices for T / NK and myeloid cell subpopulations, respectively, to complete deconvolution analysis of bulk samples: The first method first screened differentially regulated genes at the RNA level (logFC≥0.25, min.pct≥0.1), and further screened high-confidence marker genes (logFC>0.5, expression row sum>1). Based on the high-confidence genes, a mean version of the Pseudobulk reference matrix was constructed, and this matrix was used as a feature reference to run CIBERSORT deconvolution (perm=0, QN=TRUE) on the DLBCL and BL bulk datasets to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample; The second method involves stratified sampling of T / NK and myeloid cell subpopulations (9% of large / medium subpopulations, 20% of small subpopulations, and all cells of the very small subpopulation are retained, with at least 30 cells retained in each subpopulation). The resulting subset counts matrix is ​​used as a reference feature, and CIBERSORT deconvolution (perm=0, QN=TRUE) is run on 9 DLBCL and BL bulk datasets (such as GSE10846 and TCGA) to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. The result is the intersection. 3. Survival Analysis: To determine the clinical relevance of cell type infiltration, the relationship between the infiltration ratios of all T / NK and myeloid cell subpopulations and patient survival (OS) was examined in all bulk datasets. Ultimately, it was found that the infiltration ratios calculated by both methods appeared in multiple bulk datasets. The cell subpopulations that were significantly associated with survival analysis results were NK_Active_CCL3+ and CD4_Tem_Th1like_MKI67, and the enrichment of both subpopulations was significantly associated with better patient survival prognosis.

Claims

1. A method for screening biomarkers for B-cell lymphoma, characterized in that, The steps of the method are as follows: Step 1: Select data: scRNA-seq data of Burkitt lymphoma, diffuse large B-cell lymphoma and follicular lymphoma, and splenic marginal zone lymphoma, and bulk RNA-seq data of Burkitt lymphoma and diffuse large B-cell lymphoma databases; Step 2: Use the R package Seurat to filter cells based on quality control indicators, and use the LogNormalize method to normalize the data with a scaling factor of 10,000 to correct for changes in sequencing depth between units. Step 3: Using the data obtained in Step 2, select 5000 highly variable genes as the core features for downstream dimensionality reduction and clustering, perform principal component analysis, and then generate an Elbow Plot to determine the optimal number of principal components during clustering. Based on the optimal number of principal components, construct a K-nearest neighbor graph and optimize the neighborhood connection weights using Jaccard similarity. Step 4: The cell data obtained in Step 3 is processed by unsupervised cell clustering to separate cells into groups. Cells showing mixed lineage expression patterns are excluded using the marker gene method. Finally, single-cell data are obtained for subsequent merging analysis. Step 5: Add batch data, patient information, and disease type data as annotation columns, and exclude all low-quality cells; Step 6: Use the R package decontX to remove environmental RNA contamination, retaining all cells with contamination levels <0.2; Step 7: Normalize using NormalizeData in Seurat with a scaling factor of 10,000 as the parameter, detect consistent HVG across datasets using FindVariableFeatures, select the top 2,000 features for data integration, scale and center gene expression using ScaleData, and remove batch effects between data batches and patients using the RunHarmony function in the Harmony package. Step 8: Retain the first 30 principal components, use Seurat's FindNeighbors to calculate the Euclidean distance in the Harmony dimensionality reduction space to construct the KNN graph, use FindClusters to perform unsupervised clustering at a resolution of 0.5, further reduce the dimensionality through RunUMAP in Seurat, and visualize it through unified manifold approximation and projection; Step 9: Screen the cell population using marker gene methods and the SingleR tool, take the intersection data, and remove cells that show mixed lineage expression to ensure that confounding cells are excluded; Step 10: Use the R package DoubletFinder to calculate the number of mixed pairs in each subpopulation, use paramSweep to scan parameters, find the optimal pK value by combining find.pK with the BCreal index, use the clustering results as annotations, and use modelHomotypic to estimate the proportion of homologous pairs; calculate the expected number of heterologous pairs based on 14.5% of the total number of cells, and adjust it according to the homologous proportion, and finally identify and delete pairs in each cluster; Step 11: Extract cells of each type individually for a second round of dimensionality reduction clustering to further characterize the subsets of T / NK cells, myeloid cells, and fibroblasts: Step 12: Based on gene markers and differentially expressed genes, 5 NK cell subsets, 12 CD8 T cell subsets, 7 CD4 T cell subsets, 8 myeloid cell subsets, and 4 fibroblast subsets were identified. Step 13: Use cibersort to deconvolve all bulk RNA-seq datasets to quantify the abundance of cell types. It was found that myeloid cells and T / NK were most significantly correlated. Spearman correlation analysis was performed on the infiltration patterns of cell types, and it was found that there was a significant positive correlation between T / NK and myeloid cells in all 5 independent cohorts. Step 14: Based on the T / NK and myeloid cell data obtained in Step 12, differentially regulated genes at the RNA level were screened using the criteria of logFC ≥ 0.25 and min.pct ≥ 0.

1. High-confidence marker genes were further screened using the criteria of ogFC > 0.5 and expression line sum > 1. A mean version of the Pseudobulk reference matrix was constructed based on the high-confidence genes. CIBERSORT deconvolution was run on the DLBCL and BLbulk datasets to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. Step 15: Perform stratified sampling on the T / NK and myeloid cell subpopulations obtained in Step 12, extract the subset counts matrix, and use the full subpopulation Pseudobulk matrix as the reference feature. Similarly, run CIBERSORT deconvolution on the 9 DLBCL and BL bulk datasets to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. Step 16: Take the intersection of the results of Step 14 and Step 15 to obtain the relationship between T / NK and the infiltration ratio of all cell subpopulations of myeloid cells and patient survival.

2. The method according to claim 1, characterized in that, The quality control indicator in step 2 is: only retain cell data with 200-6000 detected genes and a mitochondrial gene ratio of less than 15%.

3. The method according to claim 1, characterized in that, The indicators excluded in step 5 are: nFeature_RNA < 200, nCount_RNA > 50,000, and percentage.mt > 15.

4. The method according to claim 1, characterized in that, The specific steps of step 6 are as follows: The cell contamination rate is calculated using a binary mixture model through the decontX function. A series of contamination thresholds are tested, ranging from 0.1 to 0.25, and the real-time visualization of cell retention under different thresholds is achieved using the Shiny interactive tool.

5. The method according to claim 1, characterized in that, In step 7, the function RunHarmony uses the parameters theta = 1, lambda = 1.5, sigma = 0.1 to perform two iterations; in step 8, the parameters for visualization are dims = 1:30, reduction = 'harmony'.

6. The method according to claim 1, characterized in that, In step 10, the paramSweep scan parameters are PCs=1:40 and sct=FALSE.

7. The method according to claim 1, characterized in that, In step 11, for bone marrow cells, the first 29 PCs from 2,000 highly variable genes were used at a resolution of 0.3; for T / NK cells, the first 40 PCs from 2,000 highly variable genes were extracted at a resolution of 0.

5.

8. The method according to claim 1, characterized in that, In step 14, CIBERSORT deconvolution is run on the DLBCL and BL bulk datasets with parameters perm=0 and QN=TRUE to obtain the proportion of T / NK and myeloid cell subpopulations in each bulk sample. In step 15, stratified sampling is performed by sampling 9% from large / medium subpopulations, 20% from small subpopulations, and retaining all cells from very small subpopulations, with at least 30 cells retained from each subpopulation.

9. A method for screening biomarkers for B-cell lymphoma, characterized in that, The device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements a method for screening B-cell lymphoma biomarkers as described in any one of claims 1-8.

10. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method for screening B-cell lymphoma biomarkers as described in any one of claims 1-8.