A single-cell sequencing data analysis method, device, medium, and program product

Abstract

The present invention belongs to the field of intelligent medicine, and specifically relates to a single-cell sequencing data analysis method, system, and device. The method includes obtaining single-cell sequencing data; calculating a multi-modal affinity matrix based on at least two types of modal information in the single-cell sequencing data; and performing non-negative matrix factorization on the multi-modal affinity matrix to obtain a latent variable matrix W. This application integrates multi-modal information to obtain a low-dimensional representation with consistent multi-modal information. When the latent variable matrix W is applied to scenarios such as cell subset analysis, it exhibits superior performance compared to other state-of-the-art methods.

Description

Technical Field

[0001] The present invention relates to the field of intelligent medical technology, and more particularly, to a single-cell sequencing data analysis method, device, medium, and program product. Background Art

[0002] Cancer cells undergo various genetic and epigenetic changes that shape their characteristics and lead to the formation of distinct subpopulations during tumor progression. Accurately deciphering key cell subpopulations is crucial for formulating effective treatment strategies. Compared with traditional bulk studies, single-cell sequencing technology has greatly enhanced our ability to unravel the complex cellular composition of a given tumor. For example, single-cell RNA sequencing (scRNA-seq) is a widely used method for analyzing clinically relevant subpopulations by leveraging single-cell gene expression. scRNA-seq can also classify malignant and non-malignant cell populations by analyzing expression-derived genetic mutations such as copy number variations (CNVs). However, while scRNA-seq is promising in examining cell heterogeneity, it often fails to elucidate the regulatory mechanisms that define the transcriptomic landscape.

[0003] Single-cell epigenomic assays, such as single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq), are capable of robustly analyzing cell subpopulations beyond gene expression and can reveal diverse gene regulatory information that controls gene expression. However, the application of scATAC-seq in studying tumor subpopulations is computationally challenging. Traditional methods mainly rely on CNVs to determine the genetic heterogeneity of tumor cells. Despite these advancements, CNV alone often fails to detect key tumor subpopulations. For example, melanoma subpopulations with different anti-PD1 responses can only be distinguished based on their mutation profiles rather than CNV events. In addition, the interplay between genetic and epigenetic changes enables subpopulations to evade treatment barriers and drive cancer progression, highlighting the need to integrate these features. Unfortunately, existing methods have failed to effectively utilize the full spectrum of epigenetic and genetic information from scATAC-seq data. Moreover, scATAC-seq data is typically affected by inherent high sparsity and technical noise, which poses a significant challenge in determining informative features for dissecting tumor subpopulations. Summary of the Invention

[0004] In view of the above problems, the present invention provides a single-cell sequencing data analysis method called multi-modal based single-cell sequencing data analysis (MAAS), which accurately identifies tumor subpopulations by leveraging and integrating informative multi-modal features such as CNV, SNV, and chromatin accessibility data.

[0005] The first aspect of the present application discloses a single-cell sequencing data analysis method, the method comprising:

[0006] S1: Obtain single-cell sequencing data;

[0007] S2: Calculate a multimodal affinity matrix based on at least two of the modal information in the single-cell sequencing data;

[0008] S3: Perform non-negative matrix factorization on the multimodal affinity matrix to obtain a latent variable matrix W.

[0009] Furthermore, when performing non-negative matrix factorization on the multimodal affinity matrix, integrate the multimodal affinity matrix to obtain the latent variable matrix W, where the integration is to solve for the latent variable matrix W in a loss function that includes the multimodal affinity matrix;

[0010] Optionally, the loss function for solving the latent variable matrix W is expressed as follows:

[0011]

[0012] where W represents the latent variable matrix, A (i) represents the multimodal affinity matrix, i ∈ (1, k) represents traversing k multimodal affinity matrices, H (i) represents the diagonal coefficient matrix when the affinity matrix A (i) is mapped to the latent variable matrix W, W T is the transpose of the latent variable matrix W, represents the Frobenius norm;

[0013] Optionally, solve the latent variable matrix W using the multiplicative update rule of stochastic gradient descent;

[0014] Optionally, the multiplicative update rule is expressed as follows:

[0015]

[0016] where, represents the learning rate for updating W, η represents the learning rate for updating H (i) W represents the latent variable matrix, A (i) represents the multimodal affinity matrix, i ∈ (1, k) represents traversing k multimodal affinity matrices, H (i) represents the diagonal coefficient matrix when the affinity matrix A (i) is mapped to the latent variable matrix W.

[0017] Furthermore, the single-cell sequencing data includes any one or more of the following: scATAC-seq, scDNA-seq, scRNA-seq, single-cell methylation, single-cell Hi-C, single-cell surface protein expression data;

[0018] Optionally, the types of the modal information include any two or more of the following: chromatin accessibility, copy number variation, single nucleotide variation, gene expression, DNA methylation, chromosome crosslinking, surface protein expression;

[0019] Optionally, the multi-modal affinity matrix includes any two or more of the following affinity matrices: chromatin accessibility affinity matrix, copy number variation affinity matrix, single nucleotide variation affinity matrix;

[0020] Optionally, the calculation method of the multi-modal affinity matrix includes:

[0021] Calculating the chromatin accessibility affinity matrix using the cosine distance for chromatin accessibility;

[0022] Calculating the copy number variation affinity matrix using the hamming distance for copy number variation;

[0023] Calculating the single nucleotide variation affinity matrix using the hamming distance for single nucleotide variation;

[0024] Calculating the gene expression affinity matrix using the cosine distance for gene expression;

[0025] Calculating the surface protein expression affinity matrix using the Cosine distance for surface protein expression;

[0026] Calculating the DNA methylation affinity matrix using the Cosine distance for DNA methylation;

[0027] Calculating the surface protein expression affinity matrix using the Cosine distance for surface protein expression;

[0028] Calculating the surface protein chromosome crosslinking affinity matrix using the Cosine distance for chromosome crosslinking;

[0029] Optionally, the affinity matrix is self-expressive, and the self-expression is represented as:

[0030] A (i) ~A (i) WH (i) W T

[0031] Wherein, A (i) represents the affinity matrices of different modalities, W represents the latent variable matrix, and H (i) represents the diagonal coefficient matrix when the affinity matrix A (i) is mapped to the latent variable matrix W, W T is the transpose of the latent variable matrix W, ~ represents approximately equal, i = 1, 2, …, n, and n is the number of modalities;

[0032] Optionally, the chromatin accessibility affinity matrix calculated using the cosine distance for chromatin accessibility is represented as follows:

[0033] A pq = 1 - cosine(x p , x q )

[0034] where x p represents the number of peaks of cell p, x q represents the number of peaks of cell q, cosine() represents calculating the cosine correlation, and A pq represents the element in the chromatin accessibility affinity matrix;

[0035] Optionally, the copy number variation affinity matrix calculated using the Hamming distance for copy number variation is represented as follows:

[0036] A pq = 1 - Hamming(x p , x q )

[0037] where x p represents the copy number variation of cell p, x q represents the copy number variation of cell q, Hamming() represents calculating the Hamming distance, and A pq represents the element in the copy number variation affinity matrix;

[0038] Optionally, after correcting the abnormal high peaks of chromatin accessibility, the chromatin accessibility affinity matrix is calculated. The correction of the abnormal high peaks of chromatin accessibility is represented as follows:

[0039]

[0040] where x pj represents the number of peaks of cell p before correction, represents the number of peaks of cell p before correction, j represents region j, r pε(j) represents the observed copy number of the region containing chromatin accessibility region j, and represent the regions of copy number gain and copy number reduction in cell p respectively, and ε represents the copy number region covering region j;

[0041] Optionally, noise reduction is performed on the copy number variation and / or single nucleotide variation to obtain a low-rank matrix. Based on the low-rank matrix of the copy number variation and / or single nucleotide variation, the copy number variation affinity matrix and / or single nucleotide variation affinity matrix is calculated using the Hamming distance;

[0042] Optionally, the low-rank matrix is obtained by robust principal component analysis;

[0043] Optionally, one or more of the following tools or software are used to calculate the copy number variation: InferCNV, epiAneufinder, Copy-scAT;

[0044] Optionally, one or more of the following tools or software are used to calculate the single nucleotide variation: HaplotypeCaller, Mutect2, Samtools, BCFtools, FreeBayes;

[0045] Optionally, one or more of the following tools or software are used to calculate the gene expression: Seurat, Scanpy, Scran, SingleCellExperiment;

[0046] Optionally, one or more of the following tools or software are used to calculate the DNA methylation: methylKit, MethyPipe, Methylscaper;

[0047] Optionally, one or more of the following tools or software are used to calculate the surface protein expression: Seurat, CITE-seq-Count, SingleCellExperiment;

[0048] Optionally, one or more of the following tools or software are used to calculate the chromosome crosslinking: scHiCExplorer, singleCellHiC, scHiCluster.

[0049] Furthermore, the method further includes S4: performing cluster analysis based on the latent variable matrix W to obtain cell subsets;

[0050] Optionally, the clustering method includes one or more of the following: hierarchical clustering, K-means clustering, fuzzy C-means clustering;

[0051] Optionally, the optimal number of clusters for the cluster analysis is selected based on the silhouette index, Davies-Bouldin index, Dunn validity index, and Calinski-Harabasz index;

[0052] Optionally, the guiding score for the optimal number of clusters for the cluster analysis is represented as follows:

[0053]

[0054] Among them, S represents the guiding score for the optimal number of clusters, and the number of clusters corresponding to the largest S is selected as the optimal number of clusters. SI represents the silhouette index, DBI represents the Davies-Bouldin index, DVI represents the Dunn validity index, and CHI represents the Calinski-Harabasz index;

[0055] Optionally, the contribution of each modality to clustering is quantified based on the diagonal coefficient matrix;

[0056] Optionally, the quantification of the contribution of each modality to clustering based on the diagonal coefficient matrix is expressed as:

[0057] contribution of modality i=trace(H (i) )

[0058] Among them, trace(H (i) ) represents the sum of the diagonal values, and i represents the corresponding numbers of different modalities respectively.

[0059] Furthermore, the method further includes S5: performing functional enrichment analysis based on the cell clustering results to obtain the gene pathway differences between different cell clustering results.

[0060] Furthermore, the method further includes S4': reconstructing a phylogenetic tree based on the latent variable matrix W;

[0061] Optionally, the steps of reconstructing the phylogenetic tree include: calculating a consensus affinity matrix based on the latent variable matrix W, and reconstructing the phylogenetic tree based on the consensus affinity matrix.

[0062] The second aspect of the present application discloses a method for predicting the survival rate of renal tumors, and the method includes:

[0063] S201: Obtain the sequencing data of the sample to be tested;

[0064] S202: Extract target genes based on the sequencing data;

[0065] S203: Judge the survival rate of the sample to be tested based on the expression levels of the target genes;

[0066] Optionally, the target genes include any one or more of the following: LOC388242, SULT1A4, SLX1A-SULT1A3;

[0067] Optionally, the target genes further include any one or more of the following: BOLA2, SLX1A, LOC648987, FNBP1, LINC-PINT;

[0068] Optionally, the target gene further includes any one or more of the following: NFATC2, IKZF3, CRIP1, NSD3, FAM117A;

[0069] Optionally, the target gene further includes any one or more of the following: PER1, CCDC88C, RARA, KDM6B, LINC01816, CDC42SE2;

[0070] Optionally, the target gene further includes any one or more of the following: TEDC1, MIR3194, LRRC23, JUND, C2orf42, NFKB2, ZFPM1, MIR23A, PPP1R16B, GPR68, TNFRSF1B.

[0071] A third aspect of the present application discloses a computer device, which includes: a memory and a processor; the memory is used to store program instructions; the processor is used to call the program instructions, and when the program instructions are executed, it is used to execute the steps of the above method.

[0072] A fourth aspect of the present application discloses a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the above method are implemented.

[0073] A fifth aspect of the present application discloses a computer program product or system, including a computer program, and when the computer program is executed by a processor, the steps of the above method are implemented.

[0074] A single-cell sequencing data analysis system includes:

[0075] An acquisition unit 301: used to acquire single-cell sequencing data;

[0076] An affinity matrix calculation unit 302: used to calculate a multi-modal affinity matrix based on at least two of the modal information of the single-cell sequencing data;

[0077] A latent variable matrix solving unit 303: used to perform non-negative matrix factorization on the multi-modal affinity matrix to obtain a latent variable matrix W.

[0078] Advantages of the present application:

[0079] 1. The present application innovatively discloses a single-cell sequencing data analysis method. This method first incorporates multiple different-modal data, and then finds the latent variable matrix W of the multiple different-modal data through dimensionality reduction. The latent variable matrix W realizes finding a consistent low-dimensional representation for the multi-modal information;

[0080] 2. This method can integrate and reduce the dimension of multimodal information extracted from single-omics sequencing data to find the latent variable matrix W. For example, scATAC-seq can simultaneously extract chromatin accessibility, copy number variation, and single nucleotide variation; scRNA-seq can simultaneously extract gene expression, CNV, and SNV.

[0081] 3. This method can integrate and reduce the dimension of multi-omics sequencing data to find the latent variable matrix W, thus realizing the effective utilization of multi-omics data in data analysis.

[0082] 4. After obtaining the latent variable matrix W of multimodal information based on the data analysis method of this application, when applying the latent variable matrix W to cell subpopulation analysis, compared with other state-of-the-art methods, the method of this application shows the best performance; after obtaining the latent variable matrix W of multimodal information based on the data analysis method of this application, when applying the latent variable matrix W to construct a phylogenetic tree, it shows a good construction effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present invention. For those skilled in the art, without creative efforts, other drawings can be obtained according to these drawings.

[0084] Figure 1 It is a schematic flowchart of the method provided in the first aspect of the embodiment of the present invention;

[0085] Figure 2 It is a schematic diagram of a single-cell sequencing data analysis system provided by the embodiment of the present invention;

[0086] Figure 3 It is a schematic diagram of a computer device provided by the embodiment of the present invention;

[0087] Figure 4 It is a schematic diagram of the architecture of an exemplary computing device provided by the embodiment of the present invention;

[0088] Figure 5 It is a schematic diagram of a storage medium provided by the embodiment of the present invention;

[0089] Figure 6It is a MAAS workflow provided by an embodiment of the present invention. a. The MAAS input includes a cell-by-peak matrix, a cell-by-CNV matrix, and a cell-by-SNV matrix; the original peaks are corrected according to the copy values in the corresponding regions. Robust principal component analysis is performed on the SNV to reduce noise and obtain a low-rank matrix; b. The cell similarity of each omics layer is estimated based on Euclidean or Hamming distance and is included in multimodal integration through an improved matrix factorization strategy, and the latent space that captures genetic and epigenetic features can be inferred through iterative updates; c. Tumor subpopulations are identified according to the latent variables, and the consensus cell distance can also be obtained by calculating the Euclidean distance based on the cell-by-latent factory matrix, which is included in reconstructing the neighbor-joining phylogenetic tree;

[0090] Figure 7 Schematic diagram of the application of a MAAS workflow provided by an embodiment of the present invention to the benchmark analysis of tumor subpopulation identification: a. UMAP embedding of the genetic characteristics of three subpopulations as the ground truth. b. UMAP embedding of the MAAS latent factors of the three identified subpopulations. c. Consistency of the cell distribution of the three subpopulations between the ground truth and the MAAS results. d. UMAP embedding of the genetic characteristics of four subpopulations as the ground truth. e. UMAP embedding of the MAAS characteristics of the four identified subpopulations. f. Consistency of the cell distribution of the four subpopulations between the ground truth and the MAAS results. g. Ten measurements of the accuracy of tumor cell subpopulations identified by different methods using random seeds and ARI (left), NMI (middle), and V-measure (right). These points are colored according to the integration method. h. Accuracy of tumor cell subpopulations identified by different methods, measured by ARI (left), NMI (middle), and V-measure (right) for different numbers of cells. Accuracy of tumor cell subpopulations identified by different methods, measured by ARI (left), NMI (middle), and V-measure (right) for different numbers of subpopulations. For the box plot, the center line represents the median. The y-axis of panels g-I shows the first and third quartiles, and the upper and lower whiskers extend from the hinges to the maximum or minimum values;

[0091] Figure 8Schematic diagram of the analysis of a MAAS workflow provided by an embodiment of the present invention applied to glioma subgroup identification, inferring three glioma subgroups with highly different mutation spectra: a. UMAP embedding of three tumor cell subgroups determined by MAAS. b. Copy number profiles of tumor cell subgroups identified by MAAS and CNV. Each column represents each bin (100,000 bp), and each row represents a cell. Blue indicates copy number loss; white indicates copy number neutrality; red indicates copy number gain. c. Enriched cancer hallmark features with significant differences between subgroups: red indicates enriched upregulation, and blue indicates downregulation. d. The left part shows the phylogeny of tumor subgroups and normal cells based on CNV; e. The number of mutated genes in each subgroup. f. The proportion of cells related to chemosensitivity (green) and drug resistance (orange) in each glioma cell subgroup determined by MAAS; drug sensitivity is estimated by IC50, and the lower the IC50, the more sensitive the drug response; g. The proportion of cells related to AKC treatment response (pink) and non-response (blue) in each glioma cell subgroup determined by MAAS; h. Cells related to AKC treatment response and non-response in each glioma cell subgroup predicted by random forest; the statistical P value is determined using the chi-square test;

[0092] Figure 9 MAAS identification results of ccRCC provided by an embodiment of the present invention: a. UMAP embedding of two tumor cell subgroups; b. The left part shows the phylogeny of tumor subgroups and normal cells, with driver SNVs covered in the branches, and the right panel shows the number of shared mutations and subgroup-specific mutations; c. The metastatic potential of each tumor cell subgroup is characterized by MET gene activity, metastasis signature score, and a protective factor calculated by enrichment score through metabolic-related pathways. The center line represents the median, and the boundaries show the first and third quartiles. The upper and lower whiskers extend from the hinges to the maximum or minimum values, and the distance from the hinges does not exceed 1.5 times the interquartile range. The P value is determined by the two-tailed Wilcoxon rank-sum test; d. Pseudotime-ordered ccRCC cell analysis; e. Cells are colored according to the order of differentiation, with yellow indicating less differentiation potential and blue indicating greater differentiation potential. f. Kaplan-Meier survival curves show the clinical relevance of the dedifferentiation score in two independent datasets, TCGA and E-MTAB-1980; the statistical P value is determined using the two-tailed log-rank test. G. Forest plot shows the hazard ratio and 95% confidence interval of the dedifferentiation score and clinical information according to the multivariate Cox model. Red dots represent the hazard ratio, and the horizontal bars extend from the lower limit to the upper limit of the 95% confidence interval of the hazard ratio estimate; the statistical P value is calculated by the two-tailed Wald test;

[0093] Figure 10The gene expression differences between two subpopulations identified by MAAS of ccRCC provided by an embodiment of the present invention; red represents genes specific to C1, and purple represents genes specific to C2;

[0094] Figure 11 Schematic diagram of three B-cell lymphoma subpopulations identified by MAAS provided by an embodiment of the present invention, which have highly different SNV and chromatin accessibility states. a. Principal component analysis map of three tumor subpopulations identified by MAAS; b–d: Principal component analysis maps of tumor subpopulations identified by EpiAneufinder through (b) gene expression, (c) chromatin accessibility, and (d) CNV; e. Copy number profile of tumor subpopulations identified by MAAS. Chromosomes are colored black and gray. f. Mutation frequencies repeatedly mutated in at least five cells in MAAS subpopulations, with the top bar showing the average frequency of each subpopulation; g. Decomposition of the mutation spectrum into COSMIC signatures, including three substances with known etiologies: tobacco carcinogens (SBS4), defective DNA mismatch repair (SBS15 and SBS20), and aflatoxin (SBS24); h. Marker peaks of each subpopulation, determined by comparing each subpopulation with the background of all other subpopulations for significance; i: Ranking plot showing enriched motifs in differential peaks of C1 (left), C2 (middle), and C3 (right); j: Difference (residual) between chromatin accessibility and gene expression of EBF1, with peak and expression counts being min-max normalized.

[0095] Figure 12 Schematic diagram of a method for predicting the survival rate of renal tumors provided by an embodiment of the present invention;

[0096] Figure 13 Schematic diagram of a system for predicting the survival rate of renal tumors provided by an embodiment of the present invention. Detailed implementation manners

[0097] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention.

[0098] In some of the processes described in the specification, claims, and above-mentioned drawings of the present invention, a plurality of operations appear in a specific order. However, it should be clearly understood that these operations may not be executed in the order in which they appear herein or may be executed in parallel. The serial numbers of the operations, such as S101, S102, etc., are only used to distinguish different operations, and the serial numbers themselves do not represent any execution order. In addition, these processes may include more or fewer operations, and these operations may be executed in sequence or in parallel. It should be noted that the descriptions such as "first", "second", etc. in this article are used to distinguish different messages, devices, modules, etc., do not represent a sequence, and do not limit that "first" and "second" are different types.

[0099] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative efforts belong to the scope of protection of the present invention.

[0100] Figure 1 is a schematic flowchart of a single-cell ATAC-seq data analysis method provided by an embodiment of the present invention. Specifically, the method includes the following steps:

[0101] S101: Obtain single-cell sequencing data;

[0102] In one embodiment, the scATAC-seq data of primary glioma is obtained under the accession number GSE139136 in the NCBI database (https: / / www.ncbi.nlm.nih.gov / ).

[0103] In one embodiment, the scATAC-seq data of ccRCC is obtained under the accession number PRJNA768891 in the NCBI database. The gene expression of ccRCC patients is obtained from the E-MTAB-1980 cohort of the European Bioinformatics Institute (https: / / www.ebi.ac.uk / ).

[0104] In one embodiment, from

[0105] the 10×Genomics (https: / / www.10xgenomics.com / resources / datasets / ) database to obtain the scATAC-seq data of primary B-cell lymphoma.

[0106] In one embodiment, the scATAC-seq dataset of the SNU601 cell line was obtained under the NCBI database accession number PRJNA674903, and the single-cell whole-genome sequencing data of the SNU601 cell line was obtained under the accession number PRJNA498809.

[0107] In one embodiment, paired tumor and normal samples of patients SU006 and SU008 were obtained under the NCBI database accession number PRJNA533341. Gene expression and clinical characteristics of patients in the Cancer Genome Atlas (TCGA) cohort were obtained from the GDC portal (https: / / portal.gdc.cancer.gov / ).

[0108] S102: Calculate a multi-modal affinity matrix based on at least two of the modal information in the single-cell sequencing data;

[0109] In some embodiments, the single-cell sequencing data includes any one or more of the following: scATAC-seq, scDNA-seq, scRNA-seq, single-cell methylation, single-cell Hi-C, single-cell surface protein expression data.

[0110] In some embodiments, the modal information includes any two or more of the following: chromatin accessibility, copy number variation, single nucleotide variation, gene expression, DNA methylation, chromosome crosslinking, surface protein expression, where chromatin accessibility is represented by Peaks, single nucleotide variation is represented by SNV, and copy number variation is represented by CNV.

[0111] In some embodiments, the multi-modal affinity matrix includes any two or more of the following affinity matrices: chromatin accessibility affinity matrix, copy number variation affinity matrix, single nucleotide variation affinity matrix, gene expression affinity matrix, DNA methylation affinity matrix, chromosome crosslinking affinity matrix, surface protein expression affinity matrix.

[0112] In some embodiments, the chromatin accessibility affinity matrix is calculated using the cosine distance for chromatin accessibility.

[0113] In some embodiments, the copy number variation affinity matrix is calculated using the Hamming distance for copy number variation.

[0114] In some embodiments, the single nucleotide variation affinity matrix is calculated using the Hamming distance for single nucleotide variation.

[0115] In some embodiments, the gene expression affinity matrix is calculated using the cosine distance for gene expression.

[0116] In some embodiments, a surface protein expression affinity matrix is calculated using Cosine distance for surface protein expression.

[0117] In some embodiments, a DNA methylation affinity matrix is calculated using Cosine distance for DNA methylation.

[0118] In some embodiments, a chromosome crosslinking affinity matrix is calculated using Cosine distance for chromosome crosslinking.

[0119] In some embodiments, a surface protein expression affinity matrix is calculated using Cosine distance for surface protein expression.

[0120] 1. scATAC-seq Data Processing

[0121] In one embodiment, we used SRA Toolkit (v2.10.9) to obtain FASTQ files of raw sequencing data from NCBI Sequence Read Archive, and then used 10×Genomics Cell Ranger ATAC (v2.1.0, https: / / support.10xgenomics.com / single-cell-atac / software) to align it with the GRCh38 reference genome under default parameter settings. To control data quality and obtain Peaks (cell-by-peak matrix), the fragment file was used as input to the ArchR package (v1.0.2). High-quality cells were retained based on enrichment of transcription start sites (>3) and the number of unique fragments (>3000). Doublets were identified and removed using addDoubletScores and filterDoublets, with a default resolution of 1.5. As part of the construction of the arrow file, a tiling matrix consisting of the tiling of reads in 5000 base pairs was constructed. This tiling matrix was used as input, and iterative latent semantic index dimensionality reduction was calculated using the addIterativeLSI function (set with two iterations, 20000 variable features, and 30 dimensions). After dimensionality reduction, we performed clustering using the addClusters function with a default resolution of 0.2. Pseudo-bulk group coverages based on cluster assignments were generated using the addGroupCoverages function, and peak calling was performed using the addReproduciblePeakSet function and the MACS2 peak caller (extendSummits = 2500). Finally, Peaks (cell-by-peak matrix) were constructed using the addPeakMatrix function. Additionally, a gene-by-gene score matrix for functional enrichment analysis was obtained using the addGeneScoreMatrix.

[0122] 2. Calculate peak variability

[0123] To identify genes with the greatest variance in chromatin accessibility that may affect expression, we calculated the total number of peaks (Peaks) within 2 kb regions upstream and downstream of each gene promoter in the cells. We then calculated the absolute difference in peak counts between subsets. The promoter regions (hg38) were obtained from GENCODE (https: / / www.gencodegenes.org / human / release_38.html).

[0124] 3. Calculate residuals

[0125] We considered all accessible chromatin regions within a 500 kb window around the transcription start site of EBF1. The residual was defined as the square of the difference between the normalized chromatin accessibility and the expression of the EBF1 gene.

[0126] 4. Motif enrichment between different peaks

[0127] First, addMotifAnnotation was used to determine the presence of motifs in the peak set of Peaks. Then, we used the peakAnnoEnrichment function to examine these differentially accessible peaks for enrichment of different motifs.

[0128] 5. Mutation calling

[0129] We used GATK Mutect2 and HaplotypeCaller to verify the reliability of SNV calling.

[0130] We also benchmarked two CNV callers (epiAneufinder and Copy-scAT) for scATAC-seq data, and finally we used epiAneufinder for CNV calling for MAAS analysis.

[0131] S103: Calculate a multimodal affinity matrix based on the multimodal information

[0132] 6. Low-rank approximation of the mutation matrix

[0133] To separate noise and missing values from the original data and accommodate sparsity, we used robust principal component analysis (robust PCA) to obtain a low-rank matrix of the mutation data. Specifically, our goal was to decompose the original matrix X (i) into a low-rank matrix and a noise matrix

[0134]

[0135] where λ and P Ω (X (i) ) is a linear operator:

[0136]

[0137] In the above formula, (m, n) represents the matrix size, and N(x = 0) represents the number of missing values. Under the convex condition, the constrained problem can be written as:

[0138]

[0139] This problem is solved by the inexact augmented Lagrangian multiplier algorithm.

[0140] 7. Estimation of cell affinity in each feature layer

[0141] We first corrected the Peaks (i.e., chromatin accessibility) map based on the prior knowledge that an increased copy number leads to an abnormally high peak density and vice versa:

[0142]

[0143] where x pj represents the number of peaks of cell p before correction, represents the number of peaks of cell p before correction, j represents region j, and r pε(j) represents the observed copy number of the chromatin accessibility region j, and represent the regions of copy number gain and copy number reduction in cell p respectively, and ε represents the copy number region;

[0144] Then, we calculated the affinity (default: cosine) between cells p and q. The Hamming distance was used to estimate cell similarity based on CNV or SNV:

[0145] A pq = 1 - Hamming(x p , x p )

[0146] where x p represents the CNV or SNV of cell p, and x q represents the CNV or SNV of cell q,

[0147] Hamming() represents calculating the Hamming distance, and A pq represents the element in the CNV or SNV affinity matrix;

[0148] S103: Perform non - negative matrix factorization on the multi - modal affinity matrix to obtain the latent variable matrix W.

[0149] In one embodiment, scATAC - seq data is obtained, and the following three - modal information is extracted: Peaks, SNV, CNV; then the cosine distance is used to obtain the Peaks affinity matrix, the Hamming distance is used to obtain the SNV and CNV affinity matrices, non - negative matrix factorization is performed on the affinity matrices of these three modalities to obtain the latent variable matrix W, and the latent variable matrix W is used for clustering of cell sub - populations.

[0150] In one embodiment, scRNA - seq data of liver cancer tumors is obtained, and the following three - modal information is extracted: gene expression data, SNV, CNV; then the cosine distance is used to obtain the gene expression affinity matrix, the Hamming distance is used to obtain the SNV affinity matrix and the CNV affinity matrix, non - negative matrix factorization is performed on the affinity matrices of these three modalities to obtain the latent variable matrix W, and the latent variable matrix W is used for the clustering analysis of liver cancer cell sub - populations.

[0151] In one embodiment, scRNA - seq data and scATAC - seq data of gastric cancer are obtained, and the following multi - modal information is extracted: gene expression, Peaks, SNV, CNV; then the cosine distance is used to obtain the gene expression affinity matrix and the Peaks affinity matrix, the Hamming distance is used to obtain the SNV affinity matrix and the CNV affinity matrix, non - negative matrix factorization is performed on the affinity matrices of these four modalities to obtain the latent variable matrix W, and the latent variable matrix W is used for the clustering analysis of gastric cancer cell sub - populations.

[0152] In one embodiment, scDNA - seq data is obtained, and the following multi - modal information is extracted: CNV, SNV; then the Hamming distance is used to obtain the CNV affinity matrix and the SNV affinity matrix, non - negative matrix factorization is performed on the affinity matrices of these two modalities to obtain the latent variable matrix W, and the latent variable matrix W is a consistent low - dimensional representation of the multi - modal affinity matrix.

[0153] In one embodiment, single - cell methylation, scRNA - seq data, and scDNA - seq data are obtained, and the following multi - modal information is extracted: DNA methylation, gene expression matrix, SNV, CNV; then the cosine distance is used to obtain the DNA methylation affinity matrix and the gene expression affinity matrix, the hamming distance is used to obtain the SNV affinity matrix and the CNV affinity matrix, non - negative matrix factorization is performed on the affinity matrices of these four modalities to obtain the latent variable matrix W.

[0154] In one embodiment, scATAC-seq, scRNA-seq, single-cell methylation, single-cell Hi-C, and single-cell surface protein expression data are obtained; the extracted modal information includes the following seven types: Peaks, CNV, SNV, gene expression, DNA methylation, surface protein expression, and chromosome crosslinking; after calculating the corresponding affinity matrices in sequence, non-negative matrix factorization is performed on the affinity matrices of the seven modalities to obtain the latent variable matrix W.

[0155] 8. MAAS Structure

[0156] We use improved non-negative matrix factorization to jointly integrate multiple modalities to achieve dimensionality reduction of the affinity matrix A (i) The model aims to find a consistent low-dimensional space W to simultaneously represent different layers (each modality represents a different feature layer), and uses a diagonal matrix H (i) to represent the coefficients of the latent factors to be projected into the space:

[0157] A (i) ~A (i) WH (i) W T

[0158] We note that A (i) As a self-expression term, it indicates that our multi-modal can learn and maintain the local structure of subspace clustering. Given the input terms, our model minimizes the loss function as follows:

[0159]

[0160] Using multiplicative update rules through stochastic gradient descent, as follows:

[0161]

[0162] where

[0163]

[0164] Therefore, we can obtain

[0165]

[0166]

[0167] Based on the derivative, the regular learning rate is expressed as:

[0168]

[0169] A block coordinate descent scheme is implemented, where we only optimize based on one rule and keep other rules unchanged. Finally, we achieved the decomposition using a manual solution of the equation

[0170]

[0171] When the condition is satisfied the gradient descent terminates.

[0172] 9. Contribution of each modality

[0173] We quantify the contribution of each modality by calculating the trace (sum of diagonal values).

[0174] contribution of modality i=trace(H (i) )=‖H (i) ‖1

[0175] The larger the trace, the higher the weight of the clustering assignment.

[0176] 10. Cell clustering

[0177] In our study, we applied hierarchical clustering, K-means, and fuzzy C-means clustering to classify tumor subpopulations using a custom resolution. At the same time, we designed a comprehensive score S based on existing metrics to guide the optimal number of clusters:

[0178]

[0179] This metric evaluates the distance between internal objects in a subpopulation and distant objects in different subpopulations. To simplify the presentation, we omit the superscripts. The detailed introduction of each term is as follows:

[0180] 1. Silhouette Index (SI). For sample i, the formula for SI is

[0181]

[0182] a(i) represents the average distance between sample i and any other sample within the subpopulation, and b(i) represents the minimum average distance between sample i and other subpopulations.

[0183] 2. Davies-Bouldin Index (DBI). For subpopulations i and j, the DBI formula is:

[0184]

[0185] where represents the average distance between the samples in i and the clustering center of; δ(C i ,C j ) represents the distance between the clustering centers of i and j. The smaller the DBI, the more satisfactory the partition.

[0186] 3. Dunn Validity Index (DVI). For subgroups i and j, DVI is defined as:

[0187]

[0188] represents the maximum distance between any two samples within subgroup u.

[0189] 4. Calinski-Harabasz Index (CHI). This index uses covariance to evaluate clustering performance as follows:

[0190]

[0191] B k is the between-subgroup covariance matrix, and F k is the within-subgroup covariance matrix. Additionally, n and k represent the sample size and the number of clusters. The partition with the maximum S will be selected for comparison with our model. Uniform Manifold Approximation and Projection (UMAP) embedding is generated based on the best principal components determined using the R package findPC.

[0192] 11. Phylogenetic Reconstruction

[0193] First, calculate the consensus affinity matrix based on the latent variable matrix W. Alternatively, it can be calculated based on single-modal features. The phylogenetic tree is reconstructed using the neighbor-joining algorithm implemented in the R package ape.

[0194] 12. Simulation Settings

[0195] We randomly generated 5000 accessible chromatin regions, 800 copy number regions, and 300 SNVs. We first created a cell correlation matrix ranging from 0.3 to 1. Using the R package MASS, the accessible chromatin and copy number regions follow a multivariate Gaussian distribution with a mean of 4. The SNVs for each cell are generated according to a binomial distribution with a probability of 0.3, and the correlation matrix is composed of the input covariance matrix of the mvnorm function. To generate subgroups with unique accessible chromatin profiles and no CNVs, we increased the count of 10 regions in the cell subgroup with log2 fold change (FC) values in the range of 0.65 to 1.5 compared to other cells. Thus, this subgroup can be considered an additional subgroup that cannot be distinguished by SNVs or CNVs. Precision is defined as the true positive of the real cells identified by MAAS among all cells within the same subgroup.

[0196] 13. Benchmark Analysis

[0197] We used Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN) to estimate the performance of the UMAP method. Specifically, we combined three unimodal matrices into a meta-matrix for UMAP analysis. We did not use other distance-based clustering methods based on UMAP embedding coordinates because the distances between cells could not be directly interpreted. We also compared the accuracy of MAAS with that of five multi-model-based integrated clustering tools, including intNMF, PintNMF, SNF, LRACluster, and MCIA, using default parameters. The number of cells and subpopulations was initially set to 400 and 3, respectively. Each method was run 10 times with different random seeds. We also compared the robustness of each method by varying the number of cells (800, 1500, 2000, 2500, and 3000) and subpopulations (from 4 to 8). We used three metrics, including Normalized Mutual Information (NMI), ARI, and V-measure, to evaluate unsupervised clustering performance. These metrics evaluate the difference between the predicted clustering and the reference labels by quantifying scores. The ranges of AMI and V-measure are from 0 to 1, and the range of ARI is from -1 to 1. A higher score indicates a closer match between the clustering assignment obtained considering the inferred partition C = {C : , C1, …, C ＝} and the reference classes U = {U : , U, obtaining a closer match between the clustering assignment and the known clustering assignment for 1, …, U >}. The details of each metric are as follows.

[0198] 1. AMI. This metric only corrects for the consistency effect due to chance in the mutual information (MI) between subpopulations.

[0199]

[0200] where and represent the probability that an object comes from class C i or U j , and then we can obtain the representation of AMI:

[0201]

[0202] E[·] and H(·) represent expectation and Shannon entropy, respectively.

[0203] 2. ARI. Similarly, this metric normalizes the Rand index (RI) to ensure a value closer to 0 for a random partition.

[0204]

[0205] where a represents the number of sample pairs in the same subpopulation in C and U; b represents the number of inconsistent sample pairs.

[0206] ARI is expressed as

[0207]

[0208] 3. V - measure. This metric calculates the harmonic mean of homogeneity and completeness.

[0209]

[0210] h represents homogeneity and is defined as

[0211]

[0212] where N i→j represents the number of samples partitioned into U i in the cluster C j Similarly, c is defined as

[0213]

[0214] 14. Functional enrichment analysis

[0215] We used the R package ArchR to predict gene scores and motifs for activity. Additionally, we applied GSVA to calculate enrichment scores for pathway signatures of single cells based on curated gene sets (http: / / www.gsea - msigdb.org / gsea / msigdb / ). Differentially enriched pathways were defined by the R package limma with an adjusted P - value < 0.05 and a default threshold of |log2FC| > 0.075.

[0216] 15. Chemosensitivity and immunotherapy response analysis

[0217] We used the R package pRRophetic to estimate IC50 based on the gene expression profiles of the TCGA - GBM cohort and used the median as the cut - off to divide patients into a low - IC50 group and a high - IC50 group. The low - IC50 group indicates a subpopulation of patients sensitive to the drug. Then, we applied Scissor to predict cells significantly associated with sensitive and insensitive chemotherapeutic responses to each drug. For AKC treatment analysis, we collected NanoString nCounter panel datasets of 14 recurrent GBM patients, including 5 responders and 7 non - responders. We used the randomForest function of the R package randomForest and default parameters based on existing biomarkers to train a random forest classifier. We applied this model to predict the AKC treatment response of each cell.

[0218] 16. Survival analysis

[0219] We conducted multiple analyses to study the prognostic relevance of tumor subpopulation-specific genes. The R package survminer (https: / / cran.r-project.org / web / packages / survminer / index.html) was used. The survival curves of two patients were evaluated for each group using the Kaplan-Meier method. Statistical significance

[0220] was calculated using the two-tailed log-rank test. We performed multivariable Cox analysis and Wald test using the survival R package.

[0221] We applied MAAS to glioma tumors and identified a resistant subpopulation that was previously unknown and associated with clinical outcomes. Additionally, MAAS was able to discover new B cell lymphoma subpopulations by integrating scRNA-Seq data. Moreover, MAAS could even identify normal cell subpopulations in tumor tissues. In summary, MAAS is a reliable and powerful tool for identifying tumor subpopulations from single-cell sequencing data such as scATAC-Seq data, helping to uncover new disease mechanisms and improve tumor diagnosis and treatment strategies.

[0222] MAAS achieved excellent accuracy in predicting tumor subpopulations

[0223] To characterize cell heterogeneity in tumors, we developed an algorithm called MAAS to accurately identify tumor subpopulations by integrating informative genetic and epigenetic features obtained from scATAC-seq data ( Figure 6 method a). To ensure that our results were not biased by CNV, which usually affects the quantification of chromatin accessibility, MAAS implemented a weighted correction strategy to adjust for this confounding effect. Additionally, since SNVs derived from scATAC-seq data can be very sparse and noisy, MAAS used improved robust principal component analysis and inexact augmented Lagrangian multiplier algorithms to accurately detect SNV information ( Figure 6 a). Then, we used the cosine distance of chromatin accessibility and the Hamming distance of CNV and SNV to estimate cell similarity and integrated them using multimodal non-negative matrix factorization ( Figure 6 b). To further effectively cluster cells, MAAS decomposed cell similarity into a latent variable matrix W and a diagonal coefficient matrix H and used the correlation matrix A as self-expression for better clustering. To obtain optimal results, MAAS adopted a multiplicative update rule to minimize the value of the loss function. Finally, tumor cells were classified into subpopulations using K-means clustering, and a phylogenetic tree across subpopulations was constructed ( Figure 6 c).

[0224] To systematically evaluate the performance of MAAS, we conducted simulation analyses to assess whether MAAS could accurately deconvolute the labeled tumor subpopulations. We first generated three simulated cell subpopulations as the ground truth dataset, where Subpopulations 1 and 2 had different genetic characteristics from Subpopulation 3, and Subpopulations 1 and 2 had different chromatin accessibility characteristics ( Figure 7 method of a). MAAS could successfully separate Subpopulations 1 and 2 that could not be distinguished by unimodal methods ( Figure 7 a, b). Specifically, MAAS identified 83.3%, 100%, and 98.28% of the three subpopulations ( Figure 7 c), and in addition, MAAS accurately recovered 98.7% of the differences in accessible chromatin regions between Subpopulations 1 and 2. In addition, MAAS showed comparable performance with four simulated cell subpopulations, accurately identifying 97.8% and 97% of the cells in Subpopulations 1 and 2 respectively, and correctly distinguishing all cells in Subpopulations 3 and 4 ( Figure 7 d - f).

[0225] To further evaluate the robustness of the MAAS method, we compared it with other multi - omics integration and clustering tools, including Uniform Manifold Approximation and Projection (UMAP), intNMF, PintNMF, SNF, LRACluster, and MCIA methods. We first randomized the tumor cell subpopulations and evaluated the performance using three metrics: Adjusted Rand Index (ARI), Normalized Mutual Information (NMI), and V - measure (method). MAAS showed significantly better performance than UMAP multimodal clustering and other integration methods, with median ARI, NMI, and V - measure values of 0.992, 0.981, and 0.981 respectively ( Figure 7 g). In addition, we varied the number of cells in the tumor cell subpopulations in each simulation and found that MAAS obtained the highest ARI, NMI, and V - measure scores ( Figure 7 h). We also examined the effect of the number of subpopulations on subpopulation performance and found that MAAS showed the best performance among various subpopulations ( Figure 7 i). We also evaluated the clustering performance at different data sparsities (from 10% to 90%) and found that MAAS consistently showed excellent performance, with a classification rate exceeding 50% even at 90% data sparsity.

[0226] In addition, we benchmarked MAAS against CNV estimated from single-cell whole-genome sequencing data of gastric cancer (SNU601 cell line). MAAS not only accurately identified the subpopulations detected by CNV, including amplifications on chromosomes 1 and 3 and deletions on chromosomes 4 and 18. Additionally, it identified a new subpopulation, C2b, which had 184 differentially accessible regions. Compared with the traditional unimodal predicted subpopulation C2a. To investigate the potential biological functions of the new subpopulation, we performed gene ontology analysis on genes located near these differentially chromatin-accessible regions. We found that biological processes related to mitophagy, maintenance of RNA localization, and dephosphorylation were significantly enriched. These processes have previously been implicated in gastric cancer suppression, while dephosphorylation is involved in tumor cell proliferation. Overall, MAAS has an advantage in predicting tumor cell subpopulations compared to state-of-the-art methods.

[0227] MAAS delineates a new glioma subpopulation with drug resistance

[0228] Glioblastoma is the most common and aggressive primary brain malignancy in adults. To characterize the heterogeneity of glioblastoma, we applied MAAS to a scATAC-seq dataset of adult glioblastoma containing 703 tumor cells. MAAS reproduced the subpopulations identified by traditional unimodality ( Figure 8 a), as well as new subpopulations (Clusters) that could not be distinguished by traditional methods before ( Figure 8 C2 in b).

[0229] To further characterize these new subpopulations, we performed functional enrichment analysis. We found that several well-known cancer hallmark pathways were significantly enriched in Subpopulation 2, such as myogenesis (moderated t-test, adjusted P-value = 9.56×10 -161 ) and epithelial-mesenchymal transition (moderated t-test, adjusted P-value = 7.41×10 -133 ). In contrast, compared with Subpopulation 2, several cancer hallmark pathways were highly enriched in Subpopulation 3, such as MYC target V1 (moderated t-test, adjusted P-value = 5.18×10 -79 ) and DNA repair (moderated t-test, adjusted P-value = 7.42×10 -73 )

[0230] To describe tumor progression, we constructed a phylogenetic tree to reveal tumor progression (Fig. 8d). We found that 64 mutant genes were present in all three subpopulations (not present in normal tissue cells), including EGFR and PDGFRA. Additionally, we also identified several subpopulation-specific SNVs, including SNVs related to AUTS2 in Subpopulation 1 and SNVs related to KCNJ12 and MACF1 in Subpopulation 2 (Figure 8 e). These genes have been reported as recurrently mutated genes.

[0231] To further define the clinical characteristics of the subgroups identified by MAAS, we first related the gene activities of each subgroup to the half-maximal inhibitory concentration (IC50) of first-line glioma chemotherapeutic drugs (Methods). Our findings showed that subgroup 3 was more sensitive to methotrexate and cisplatin ( Figure 8 f) Chi-square test, P values = 1.23 × 10 -22 and 0.002). Then we analyzed the responses of different cell subgroups to cell-targeted immunotherapy. Our results showed that subgroup 3 was most likely to respond to activated autologous killer cell (AKC) treatment compared to subgroup 2 ( Figure 8 g; Chi-square test, P value = 9.43 × 10−6). To predict the AKC treatment responses of each cell, we trained an accurate machine learning model. Our model showed that subgroup 2 (1.8%) had a higher response rate compared to subgroup 3 (6.2%) ( Figure 8 h; Chi-square test, P value = 0.004). In summary, MAAS identified a new glioma subgroup whose drug resistance was masked by traditional single-modal methods, highlighting the potential of the MAAS method as a powerful tool for accurately classifying glioma subgroups. MAAS anatomizes B cell lymphoma subgroups by multi-omics integration

[0232] To demonstrate the utility of the MAAS method for detecting subgroups by integrating with other single-cell multi-omics data, we applied MAAS to a 10× multi-omics dataset of B cell lymphoma 31 where scRNA-seq and scATAC-seq data were measured for each cell simultaneously. Our analysis showed that MAAS accurately predicted three tumor subgroups from 2,077 tumor cells. Compared to traditional single-modal methods, MAAS better separated these cells into distinct subgroups ( Figure 11 a-d). Additionally, we compared MAAS with inferCNV and CopyKAT, two widely used methods that use gene expression-derived copy numbers to identify tumor subgroups, and found that MAAS consistently outperformed them in cleanly separating cell subgroups.

[0233] To define the characteristics of the tumor cell subgroups predicted by MAAS, we first investigated the copy number profiles of the three subgroups and found that they all exhibited low CNV heterogeneity ( Figure 11 e). However, we observed that the subgroups identified by MAAS had distinct SNV profiles and mutation frequencies ( Figure 11f. The mutation frequency of subgroup 3 was the highest at 21.0%, while that of subgroup 1 was the lowest at 18.8%, followed by group 2 (19.6%) ( Figure 10 ; Kruskal-Wallis test, P value = 2.21×10 -63). We also performed mutation signature analysis (methods) and found that the distribution of these mutations differed among subgroups and was associated with single-base substitution (SBS) signatures related to tobacco carcinogens (SBS4), defective DNA mismatch repair (SBS15 and SBS20), aflatoxin (SBS24), SBS12, and SBS30 ( Figure 11 g).

[0234] Subsequently, we identified 2,146 distinct accessible region clusters in these three regions ( Figure 11 h). To further explore the underlying molecular mechanisms, we studied the potentially bound transcription factors (TFs) in differentially accessible chromatin regions (methods) and inferred TF activity by estimating gains or losses in chromatin accessibility. Our analysis showed that several TFs exhibited significant variability among subgroups ( Figure 11 i). For example, KLF activity increased in subgroup 2, while SOX and FOXP1 showed increased activity in subgroup 3. We also found that although EBF1 is a key factor in B cell lineage specification and had the same transcriptional level among subgroups, the activity level of EBF1-accessible chromatin regions increased in subgroup 2 ( Figure 11 j; methods), indicating that changes in EBF1 chromatin accessibility occurred before EBF1 gene expression. Overall, the MAAS method was able to integrate single-cell multi-omics data to characterize B cell lymphoma subgroups.

[0235] MAAS can be extended to identify normal cell subgroups

[0236] To further explore whether the MAAS method can be extended to detect subgroups of normal cells, we studied astrocytes in the glioma dataset and CD8 + T cells, circulating T cells, and monocytes in the B cell lymphoma dataset, with at least 300 cells of each cell type. We used two metrics, intra-subgroup distance and inter-subgroup distance, to evaluate the performance of subgroup identification. The smaller the intra-subgroup distance and the larger the inter-subgroup distance, the better the performance. Our results showed that compared with unimodal methods, MAAS achieved higher inter-subgroup distance values and smaller intra-subgroup distance values. The excellent performance of MAAS in dissecting normal cell heterogeneity may be attributed to its detection of subgroup-specific mutations. Overall, these findings suggest that MAAS can also be applied to normal subgroups, thus expanding the compatibility and scalability of the method.

[0237] To our knowledge, MAAS is the first computational method for multimodal integration of scATAC-seq data, which can identify key tumor subpopulations different from those determined by traditional unimodal methods (such as Copy-scAT and epiAneufinder). We found that MAAS provides higher accuracy in identifying tumor subpopulations compared to other available methods. By integrating multimodal data, we identified new biologically and clinically relevant tumor subpopulations. The MAAS method is fundamentally different from previous subpopulation prediction methods. First, the MAAS method maximizes the number of informative features extracted from scATAC-seq data, rather than relying on a single-modal feature that ignores other important biological aspects of tumor subpopulations. In addition, the self-expressive multimodal matrix factorization strategy enhances the multimodal signal, enabling more robust classification of tumor subpopulations. Moreover, MAAS is an interpretable multimodal integration method that can quantify the contribution of each modality to cell subpopulation assignment. For example, the three glioma subpopulations predicted by MAAS are mainly driven by CNV and chromatin accessibility, with weight increases of 22.62% and 20.76% respectively compared to SNV compared to the control group. The MAAS method is also feasible because it can simultaneously examine gene mutations and epigenetic variations without the need for additional single-cell assays.

[0238] In summary, the MAAS method highlights the ability of multimodal integration to dissect tumor heterogeneity using single-cell epigenomic data. We expect MAAS to be able to widely apply to widely available single-cell sequencing data of tumors and other diseases and help reveal key tumor subpopulations for cell-targeted therapies.

[0239] Figure 12 It is a schematic diagram of a method for predicting glioma drug resistance provided by an embodiment of the present invention.

[0240] The second aspect of the present application discloses a method for predicting the survival rate of renal tumors, the method comprising:

[0241] S201: Obtain sequencing data of a sample to be tested;

[0242] S202: Extract target genes based on the sequencing data;

[0243] S203: Judge the survival rate of the sample to be tested based on the expression level of the target genes;

[0244] Optionally, the target genes include any one or more of the following: LOC388242, SULT1A4, SLX1A-SULT1A3;

[0245] Optionally, the target gene further includes any one or more of the following: BOLA2, SLX1A, LOC648987, FNBP1, LINC-PINT;

[0246] Optionally, the target gene further includes any one or more of the following: NFATC2, IKZF3, CRIP1, NSD3, FAM117A;

[0247] Optionally, the target gene further includes any one or more of the following: PER1, CCDC88C, RARA, KDM6B, LINC01816, CDC42SE2;

[0248] Optionally, the target gene further includes any one or more of the following: TEDC1, MIR3194, LRRC23, JUND, C2orf42, NFKB2, ZFPM1, MIR23A, PPP1R16B, GPR68, TNFRSF1B.

[0249] We applied the described single-cell sequencing data analysis method to renal cancer and identified a progressive tumor subset associated with poor prognosis.

[0250] MAAS identified a new renal tumor subset associated with poor prognosis

[0251] To evaluate the utility of the MAAS method in identifying tumor subpopulations with low tumor purity, we analyzed the scATAC-seq dataset of a 75-year-old male diagnosed with clear cell renal cell carcinoma (ccRCC), which contained 390 tumor cells and a tumor purity of 3.96%. Compared with traditional methods, the MAAS method better showed the differences between tumor cell subpopulations ( Figure 9 a). Functional enrichment analysis of the MAAS-predicted subpopulations showed that the MAAS method not only recovered the hallmark pathways identified by traditional chromatin accessibility-based methods but also identified subsets of new cancer pathways that had been previously overlooked. These pathways included the subset 1 enrichment pathway of MYC target V1 (adjusted t-test, adjusted P-value = 1.10×10 -20 ) and TNFA signaling through NFKB (adjusted t-test, adjusted P-value = 1.39×10 -20 ) -14 ) and the subset 2 enrichment pathway epithelial-mesenchymal transition (adjusted t-test, adjusted P-value = 1.87×10 -11 ) and angiogenesis (adjusted t-test, P-value = 2.28×10 -4 ).

[0252] MAAS also generated a phylogenetic tree that depicts the evolutionary process of tumor cell subsets. We found that subset 1 had more mutations than subset 2 ( Figure 9 b), and identified 18 mutated genes specific to subset 1, including BACH2, which suppresses immune responses through epigenetic silencing in cancer, and 13 mutated genes specific to subset 2. Additionally, we used MET gene activity to evaluate the metastatic characteristics of each subset and established unique enriched transcriptional metastatic signatures and metabolic pathways in primary ccRCC tumors. Subset 1 exhibited significantly higher MET gene activity and metastatic scores than subset 2. On the other hand, subset 2 showed a higher primary tumor signature score ( Figure 9 c; Wilcoxon rank-sum test, P values = 0.0067, 0.04, and 7.6 × 10−16, respectively), indicating that subset 1 has a more significant metastatic potential and subset 2 may be a primitive subset.

[0253] To further validate the dynamic changes between subsets predicted by MAAS, we used Monocle to construct trajectories ( Figure 9 d), showing a stepwise transition from subset 2 to subset 1. Using CytoTRACE, we determined the degree of differentiation associated with the metastatic phenotype. Subset 1 had a higher CytoTRACE score and largely dominated the end of the trajectory ( Figure 9 e). We also identified genes significantly associated with the CytoTRACE score and calculated malignancy by averaging the top-ranked gene expressions ( Figure 10 , with red being genes specific to subset 1 and purple being genes specific to subset 2). This indicates that subset 1 has highly expressed genes, suggesting a high degree of dedifferentiation. Importantly, we found that ccRCC patients with increased expression of subset 1-specific genes in the tumor showed significantly worse survival rates ( Figure 9 f, g; log-rank test, P values = 8.44 × 10 -6 and 0.0015). Collectively, our MAAS analysis identified a new progressive subset of ccRCC that is most strongly correlated with low survival rates, which may contribute to a better understanding of the underlying pathogenesis of renal cancer.

[0254] Figure 3 is a schematic diagram of a computer device provided by an embodiment of the present invention, as Figure 3 shown, the device may include: one or more processors, and one or more memories; wherein, computer-readable code is stored in the memory, and when the computer-readable code is run by the one or more processors, the above-described method can be executed.

[0255] The processor in this embodiment may be an integrated circuit chip with signal processing capabilities. The above-mentioned processor may be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components. It can implement or execute the various methods, operations, and logic block diagrams disclosed in the embodiments of the present disclosure. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc., and it may be of the X86 architecture or the ARM architecture.

[0256] Generally speaking, the various exemplary embodiments of the present disclosure may be implemented in hardware or dedicated circuits, software, firmware, logic, or any combination thereof. Some aspects may be implemented in hardware, while other aspects may be implemented in firmware or software that can be executed by a controller, a microprocessor, or other computing devices. When aspects of the embodiments of the present disclosure are illustrated or described as block diagrams, flowcharts, or using some other graphical representation, it will be understood that the blocks, devices, systems, technologies, or methods described herein may be implemented as non-limiting examples in hardware, software, firmware, dedicated circuits or logic, general hardware or a controller or other computing devices, or some combination thereof.

[0257] For example, the method or device according to the embodiments of the present disclosure may also be implemented by means of Figure 4 the architecture of the computing device 3000 shown. As Figure 4 shown, the computing device 3000 may include a bus 3010, one or more CPUs 3020, a read-only memory (ROM) 3030, a random access memory (RAM) 3040, a communication port 3050 connected to a network, input / output components 3060, a hard disk 3070, etc. The storage device in the computing device 3000, such as the ROM 3030 or the hard disk 3070, may store various data or files used for the processing and / or communication of the method provided by the present disclosure and the program instructions executed by the CPU. The computing device 3000 may also include a user interface 3080. Of course, Figure 4 the architecture shown is only exemplary, and when implementing different devices, one or more components shown in the Figure 4 computing device may be omitted according to actual needs.

[0258] The embodiment of the present invention also provides a computer-readable storage medium, such as Figure 5As shown in the figure, it is a schematic diagram of a storage medium provided by an embodiment of the present invention. Computer-readable instructions 4010 are stored on the computer storage medium 4020. When the computer-readable instructions 4010 are run by a processor, the method according to the embodiments of the present disclosure described with reference to the above figures can be executed. The computer-readable storage medium in the embodiments of the present disclosure may be a volatile memory or a non-volatile memory, or may include both volatile and non-volatile memories. The non-volatile memory may be a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or a flash memory. The volatile memory may be a random access memory (RAM), which is used as an external cache. By way of example but not limitation, many forms of RAM are available, such as static random access memory (SRAM), dynamic random access memory (DRAM), synchronous dynamic random access memory (SDRAM), double data rate synchronous dynamic random access memory (DDR SDRAM), enhanced synchronous dynamic random access memory (ESDRAM), synchronous link dynamic random access memory (SLDRAM), and direct memory bus random access memory (DR RAM). It should be noted that the memories for the methods described herein are intended to include, but are not limited to, these and any other suitable types of memories. It should be noted that the memories for the methods described herein are intended to include, but are not limited to, these and any other suitable types of memories.

[0259] The embodiments of the present disclosure also provide a computer program product or system, and when the computer program is executed by a processor, the steps of the above method are implemented.

[0260] A single-cell sequencing data analysis system, as Figure 2 shown, includes:

[0261] An acquisition unit 301: used to acquire single-cell sequencing data;

[0262] An affinity matrix calculation unit 302: used to calculate a multi-modal affinity matrix based on at least two of the modal information of the single-cell sequencing data;

[0263] A latent variable matrix solving unit 303: used to perform non-negative matrix factorization on the multi-modal affinity matrix to obtain a latent variable matrix W.

[0264] A system for predicting the survival rate of renal tumors, as Figure 13 shown, includes:

[0265] A data acquisition unit 401: acquires sequencing data of a sample to be tested;

[0266] An extraction unit 402: extracts target genes based on the sequencing data;

[0267] Prediction unit 403: Determine the survival rate of the sample to be tested based on the expression level of the target gene.

[0268] It should be noted that the flowcharts and block diagrams in the accompanying drawings illustrate the possible architectures, functions, and operations of systems, methods, and computer program products according to various embodiments of the present disclosure. In this regard, each block in the flowchart or block diagram may represent a module, a program segment, or a part of code, and the module, program segment, or part of code contains one or more executable instructions for implementing the specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in a different order than that marked in the accompanying drawings. For example, two consecutive blocks shown may actually be executed substantially in parallel, and they may sometimes be executed in the reverse order, depending on the functions involved. It should also be noted that each block in the block diagram and / or flowchart, as well as the combination of blocks in the block diagram and / or flowchart, may be implemented by a dedicated hardware-based system for performing the specified functions or operations, or may be implemented by a combination of dedicated hardware and computer instructions.

[0269] Generally speaking, various example embodiments of the present disclosure may be implemented in hardware or a dedicated circuit, software, firmware, logic, or any combination thereof. Some aspects may be implemented in hardware, while other aspects may be implemented in firmware or software that can be executed by a controller, a microprocessor, or other computing devices. When aspects of the embodiments of the present disclosure are illustrated or described as block diagrams, flowcharts, or using some other graphical representation, it will be understood that the blocks, devices, systems, technologies, or methods described herein may be implemented as non-limiting examples in hardware, software, firmware, a dedicated circuit or logic, general hardware or a controller or other computing devices, or some combination thereof.

[0270] Those skilled in the art can clearly understand that for the convenience and brevity of description, the specific working processes of the systems, devices, and units described above may refer to the corresponding processes in the foregoing method embodiments and will not be elaborated herein.

[0271] In the several embodiments provided in this application, it should be understood that the disclosed systems, devices, and methods may be implemented in other ways. For example, the device embodiments described above are merely illustrative. For example, the division of the units is only a logical function division, and there may be other division methods in actual implementation. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Another point is that the couplings or direct couplings or communication connections shown or discussed with each other may be indirect couplings or communication connections through some interfaces, devices, or units, and may be in electrical, mechanical, or other forms.

[0272] The unit described as a separation component may or may not be physically separated. The component shown as a unit may or may not be a physical unit, that is, it may be located in one place, or it may be distributed over multiple network units. Some or all of the units can be selected according to actual needs to achieve the purpose of the solution of this embodiment.

[0273] In addition, in each embodiment of the present invention, each functional unit may be integrated in a processing unit, or each unit may exist physically alone, or two or more units may be integrated in one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of a software functional unit.

[0274] Those of ordinary skill in the art can understand that all or part of the steps in the various methods of the above embodiments can be completed by instructing relevant hardware through a program, and the program can be stored in a computer-readable storage medium. The storage medium may include: read-only memory (ROM, Read Only Memory), random access memory (RAM, Random Access Memory), magnetic disk or optical disk, etc.

[0275] Those of ordinary skill in the art can understand that all or part of the steps in implementing the methods of the above embodiments can be completed by instructing relevant hardware through a program, and the said program can be stored in a computer-readable storage medium. The storage medium mentioned above may be a read-only memory, magnetic disk or optical disk, etc.

[0276] The above has introduced in detail a computer device provided by the present invention. For those of ordinary skill in the art, according to the idea of the embodiments of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A method for analyzing single-cell sequencing data, characterized in that, The method includes: S1: Obtain single-cell sequencing data; S2: Calculate a multi-modal affinity matrix based on at least two of the modal information in the single-cell sequencing data; S3: Perform non-negative matrix factorization on the multi-modal affinity matrix to obtain a latent variable matrix W, and the non-negative matrix factorization is expressed as: A (i) ~A (i) WH (i) W T Among them, A (i) represents the multi-modal affinity matrix, W represents the latent variable matrix, and H (i) represents the affinity matrix A (i) when mapped to the diagonal coefficient matrix of the latent variable matrix W, W T is the transpose of the latent variable matrix W, ∼ represents approximately equal, i = 1, 2, …, n, and n is the number of modalities. When performing non-negative matrix factorization on the multi-modal affinity matrix, integrate the multi-modal affinity matrix to obtain the latent variable matrix W, and the integration is to solve the loss function of the latent variable matrix W. The loss function for solving the latent variable matrix W is expressed as follows: Among them, W represents the latent variable matrix, and A (i) represents the multi-modal affinity matrix. i ∈ (1, k) represents traversing k multi-modal affinity matrices. H (i) represents the affinity matrix A (i) when mapped to the diagonal coefficient matrix of the latent variable matrix W, and W T is the transpose of the latent variable matrix W, represents the Frobenius norm.

2. The single-cell sequencing data analysis method according to claim 1, wherein Solve the latent variable matrix W using the multiplicative update rule of stochastic gradient descent.

3. The single-cell sequencing data analysis method according to claim 2, wherein The multiplicative update rule is expressed as follows: Among them, represents the learning rate for updating W, η represents the learning rate for updating H (i) W represents the latent variable matrix, A (i) represents the multi-modal affinity matrix, i ∈ (1, k) represents traversing k multi-modal affinity matrices, H (i) represents the affinity matrix A (i) is the diagonal coefficient matrix when mapping the affinity matrix A to the latent variable matrix W.

4. The single-cell sequencing data analysis method according to claim 1, characterized in that The single-cell sequencing data includes any one or more of the following: scATAC-seq, scDNA-seq, scRNA-seq, single-cell methylation, single-cell Hi-C, single-cell surface protein expression data.

5. The single-cell sequencing data analysis method according to claim 1, wherein The modal information includes any two or more of the following: chromatin accessibility, copy number variation, single nucleotide variation, gene expression, DNA methylation, chromosome crosslinking, surface protein expression.

6. The single-cell sequencing data analysis method according to claim 1, wherein The multi-modal affinity matrix includes any two or more of the following affinity matrices: chromatin accessibility affinity matrix, copy number variation affinity matrix, single nucleotide variation affinity matrix, gene expression affinity matrix, DNA methylation affinity matrix, chromosome crosslinking affinity matrix, surface protein expression affinity matrix.

7. The single-cell sequencing data analysis method according to claim 6, wherein The calculation method of the multi-modal affinity matrix includes: Calculate the chromatin accessibility affinity matrix using the cosine distance for chromatin accessibility; Calculate the copy number variation affinity matrix using the Hamming distance for copy number variation; Calculate the single nucleotide variation affinity matrix using the Hamming distance for single nucleotide variation; Calculate the gene expression affinity matrix using the cosine distance for gene expression; Calculate the surface protein expression affinity matrix using the Cosine distance for surface protein expression; Calculate the DNA methylation affinity matrix using the Cosine distance for DNA methylation; Calculate the chromosome crosslinking affinity matrix using the Cosine distance for chromosome crosslinking; Calculate the surface protein expression affinity matrix using the Cosine distance for surface protein expression.

8. The single-cell sequencing data analysis method according to claim 7, wherein The calculation of the chromatin accessibility affinity matrix using the cosine distance for chromatin accessibility is expressed as follows: A pq = 1 - cosine(x p , x q ) where x p represents the number of peaks of cell p, and x q represents the number of peaks of cell q. cosine() represents calculating the cosine correlation, and A pq represents an element in the chromatin accessibility affinity matrix.

9. The single-cell sequencing data analysis method according to claim 7, characterized in that The calculation of the copy number variation affinity matrix using the Hamming distance for copy number variation is expressed as follows: A pq = 1 - Hamming(x p , x q ) Among them, x p represents the copy number variation of cell p, and x q represents the copy number variation of cell q. Hamming() represents calculating the Hamming distance, and A pq represents the element in the copy number variation affinity matrix.

10. The single-cell sequencing data analysis method according to claim 8, wherein After correcting the abnormal peak values of chromatin accessibility, calculate the chromatin accessibility affinity matrix, and the correction of the abnormal peak values of chromatin accessibility is expressed as follows: where x pj represents the peak number of cell p before correction, represents the peak number of cell p after correction, j represents region j, r pε(j) represents the observed copy number of chromatin accessibility region j, and represent the regions of copy number gain and the regions of copy number loss in cell p, respectively, and ε represents the copy number region covering region j.

11. The single-cell sequencing data analysis method according to claim 5, characterized in that Denoise the copy number variation and / or single nucleotide variation to obtain a low-rank matrix, and calculate the copy number variation affinity matrix and / or single nucleotide variation affinity matrix using the Hamming distance based on the low-rank matrix of the copy number variation and / or single nucleotide variation.

12. The single-cell sequencing data analysis method according to claim 11, wherein The low-rank matrix is obtained by robust principal component analysis.

13. The single-cell sequencing data analysis method according to claim 5, wherein The copy number variation is calculated using one or more of the following tools or software: InferCNV, epiAneufinder, Copy-scAT.

14. The single-cell sequencing data analysis method according to claim 5, characterized in that The single nucleotide variation is calculated using one or more of the following tools or software: HaplotypeCaller, Mutect2, Samtools, BCFtools, FreeBayes.

15. The single-cell sequencing data analysis method according to claim 5, characterized in that The gene expression is calculated using one or more of the following tools or software: Seurat, Scanpy, Scran, SingleCellExperiment.

16. The single-cell sequencing data analysis method according to claim 5, wherein The DNA methylation is calculated using one or more of the following tools or software: methylKit, MethyPipe, Methylscaper.

17. The single-cell sequencing data analysis method according to claim 5, wherein The surface protein expression is calculated using one or more of the following tools or software: Seurat, CITE-seq-Count, SingleCellExperiment.

18. The single-cell sequencing data analysis method according to claim 5, wherein The chromosome crosslinking is calculated using one or more of the following tools or software: scHiCExplorer, singleCellHiC, scHiCluster.

19. The single-cell sequencing data analysis method according to claim 1, characterized in that The method further includes S4: performing cluster analysis based on the latent variable matrix W to obtain cell subsets.

20. The single-cell sequencing data analysis method according to claim 19, wherein, The clustering methods include one or more of the following: hierarchical clustering, K-means clustering, fuzzy C-means clustering.

21. The single-cell sequencing data analysis method according to claim 19, wherein, The optimal number of clusters for the cluster analysis is selected based on the silhouette index, Davies-Bouldin index, Dunn validity index, and Calinski-Harabasz index.

22. The single-cell sequencing data analysis method according to claim 21, wherein The guiding score for the optimal number of clusters in the cluster analysis is expressed as follows: Where S represents the guiding score for the optimal number of clusters, and the cluster number corresponding to the largest S is selected as the optimal number of clusters. SI represents the silhouette index, DBI represents the Davies-Bouldin index, DVI represents the Dunn validity index, and CHI represents the Calinski-Harabasz index.

23. The single-cell sequencing data analysis method according to claim 19, characterized in that, The method further includes S5: performing functional enrichment analysis based on the cell clustering results to obtain gene pathway differences between different cell clustering results.

24. The single-cell sequencing data analysis method according to claim 1, wherein The method further includes S4': reconstructing a phylogenetic tree based on the latent variable matrix W.

25. The single-cell sequencing data analysis method according to claim 24, wherein The steps for reconstructing the phylogenetic tree include: calculating a consensus affinity matrix based on the latent variable matrix W, and reconstructing the phylogenetic tree based on the consensus affinity matrix.

26. A computer device, characterized in that, The device includes: a memory and a processor; the memory is used to store a computer program; the processor executes the computer program to implement the steps of the method according to any one of claims 1-25.

27. A computer-readable storage medium, characterized in that, A computer program is stored thereon, and when the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1-25.

28. A computer program product, comprising a computer program, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1-25.

Images