Methods and systems for quantifying cellular activity from high-throughput sequencing data

Abstract

The present invention relates to a computer-implemented method and processing system for quantifying cellular activity from high-throughput sequencing data. The method includes the steps of generating a multimodal knowledge graph by combining gene regulatory networks (GRNs) with gene annotations derived from domain knowledge, where the vertices of the multimodal knowledge graph are genes and the gene annotations enrich the relationships between the genes, creating multiple gene modules (GMs) by clustering gene embeddings of the GRNs, embedding samples of the sequencing data into the multimodal knowledge graph, and generating, for each sample of the sequencing data, an activation vector represented as the distance between the embedding and the centroid of each of the multiple GMs.

Description

[Technical Field]

[0001] The present invention relates to a computer-implemented method and processing system for quantifying cellular activity from high-throughput sequencing data. [Background technology]

[0002] Recently, high-throughput sequencing methods (e.g., single-cell sequencing, sc-seq) have been developed that allow users to perform detailed analyses of the entire gene expression landscape (transcriptome) at single-cell resolution. Consequently, such methods have become useful techniques for investigating tumor heterogeneity (i.e., variability in expressed genes within different tumor cells of a patient) in order to predict treatment outcomes for patients before initiating therapy. [Prior art documents] [Non-patent literature]

[0003] [Non-Patent Document 1] D. van Dijk et al.: Recovering gene interactions from single-cell data using data diffusion. Cell 174(3), pp. 716–72927 (2018) [Non-patent document 2] M. Huang et al.: SAVER: gene expression recovery for single-cell RNA sequencing. Nature Methods 15, pp. 539-542 (2018) [Non-patent document 3] WV Li et al.: An accurate and robust imputation method scImpute for single-cell RNA-seq data. Nature Communications 9(997) (2018) [Non-patent document 4] G. Eraslan et al.: Single-cell RNA-seq denoising using a deep count autoencoder. Nature Communications 10(390) (2019) [Non-patent document 5] J. Rao et al.: Imputing Single-cell RNA-seq data by combining Graph Convolution and Autoencoder Neural Networks (2020). bioRxiv 2020.02.05.935296 [Non-patent document 6] X. Li et al: "Deep learning enables accurate clustering with batch effect removal in single-cell RNA-seq analysis", Nature Communications 11(2338) (2020) [Non-Patent Document 7] Sharma et al.: Longitudinal single-cell RNA sequencing of patient-derived primary cells reveals drug-induced infidelity in stem cell hierarchy. Nature Communications 9(4931) (2018) [Non-patent document 8] J. Fan, J. et al.: Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nature Methods 13, pp. 241-244 (2016) [Non-Patent Document 9] J. Ma et al.: Using deep learning to model the hierarchical structure and function of a cell. Nature Methods 15, pp. 290-298 (2018) [Non-Patent Document 10] FG Frost et al.: Pan-cancer RNA-seq data stratifies tumors by some hallmarks of cancer. Journal of Cellular and Molecular Medicine 24, pp. 418-430 (2020) [Non-Patent Document 11] B. Shickel et al.: Deep EHR: A survey of recent advances in deep learning techniques for electronic health record (EHR) analysis. IEEE Journal of Biomedical and Health Informatics 22(5), pp. 1589–1605 (2018) [Non-Patent Document 12] D. Toro-Dominguez et al.: Differential treatments based on drug-induced gene expression signatures and longitudinal systemic lupus erythematosus stratification. Scientific Reports 9(15502) (2019) [Non-Patent Document 13] Frohlich et al.: Premenopausal breast cancer: potential clinical utility of a multi-omics based machine learning approach for patient stratification. EPMA Journal 9, pp. 175-186 (2018) [Non-Patent Document 14] P. Jeyananthan et al.: Classification and regression analysis of lung tumors from multi-level gene expression data. Proceedings of the International Joint Conference on Neural Networks (2019) [Non-Patent Document 15] Y. Jang et al.: CaPSSA: visual evaluation of cancer biomarker genes for patient stratification and survival analysis using mutation and expression data. Bioinformatics 35(24), pp. 5341-5343 (2019) [Non-Patent Document 16] FJ Campos-Laborie et al.: DECO: decompose heterogeneous population cohorts for patient stratification and discovery of sample biomarkers using omic data profiling. Bioinformatics 35(19), pp. 3651-3662 (2019) [Non-Patent Document 17] L. Yu: RNA-Seq reproducibility assessment of the sequencing quality control project. Cancer informatics 19, pp. 1-5 (2020) [Non-Patent Document 18] K.-T. Kim et al.: Single-cell mRNA sequencing identifies subclonal heterogeneity in anti-cancer drug responses of lung adenocarcinoma cells. Genome Biology 16(127) (2015) [Non-Patent Document 19] Garcia-Duran, A., Niepert, M.: Learning graph representations with embedding propagation. Advances in Neural Information Processing Systems 30 (2017) [Non-Patent Document 20] Kotnis, B., Nastase, V.: Analysis of the impact of negative sampling on link prediction in knowledge graphs. KnoProceedings of the 1st Workshop on Knowledge Base Construction, Reasoning and Mining (2018) [Non-Patent Document 21] "Integrating single-cell transcriptomic data across different conditions, technologies, and species", Nature Biotechnology 36(5), pp. 411-420 (2018) [Non-Patent Document 22] Hug, N.: Surprise: A python library for recommender systems. Journal of Open Source Software 5, 2174 pages (2020) [Non-Patent Document 23] Rajkumar, SV: Multiple myeloma: 2012 update on diagnosis, risk-stratification, and management. American Journal of Hematology 87, pp. 77-88 (2012) [Non-Patent Document 24] Generic GO Slim Subset. http: / / current.geneontology.org / ontology / subsets / goslim_generic.obo [Non-Patent Document 25] McInnes, L., Healy, J., Melville, J.: UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (2018). arXiv:1802.03426 [stat.ML] [Non-Patent Document 26] Series GSE45719. https: / / www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE45719 [Non-Patent Document 27] MT Ribeiro et al.: “why should i trust you?”: Explaining the predictions of any classifier. In: Proceedings of the 22nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 2016 [Non-patent document 28] Tanaka, K.: The proteasome: Overview of structure and functions. Proceedings of the Japan Academy, Series B 85, ​​pp. 12–36, 2009 Summary of the Invention [Problem to be solved by the invention]

[0004] While high-throughput (and, to a large extent, sc-seq) data can provide significant insights into tumors, such data are affected by different forms of noise, such as batch effects and dropout events. For example, batch effects can arise when different expression values ​​of genes are observed within samples taken from the same patient. Dropouts represent the presence of ambiguous missing values ​​in the data, where it is unclear whether they result from non-expression of the gene or from non-detection of its expression value. However, the most significant problem is that when the data are used "as is" to perform predictions regarding the response of samples (e.g., patient cells, in the case of sc-seq) to specific treatments, the results are quite poor, making it impossible to make effective therapy recommendations.

[0005] The present invention therefore aims to improve and further develop methods and processing systems of the kind initially described for quantifying cellular activity from high-throughput sequencing data in such a way that a higher level of information for ML-based prediction is provided. [Means for solving the problem]

[0006] According to an embodiment of the present invention, the above object is achieved by a computer-implemented method for quantifying cellular activity from high-throughput sequencing data, the method comprising the steps of generating a multimodal knowledge graph by combining gene regulatory networks, GRNs, with gene annotations derived from domain knowledge, wherein the vertices of the multimodal knowledge graph are genes, and the gene annotations enrich the relationships between the genes; and clustering the gene embeddings of the GRNs to generate a multimodal knowledge graph by clustering the gene embeddings. multiple The steps of creating a genetic module, GM, of the sequencing data samples and embedding the sequencing data samples in a multimodal knowledge graph are as follows: multipleand generating an activation vector represented as a distance between the center of gravity of each of the GMs.

[0007] The above object is further achieved by providing a processing system for quantifying cellular activity from high-throughput sequencing data, the system comprising one or more processors, the one or more processors generating a multimodal knowledge graph by combining gene regulatory networks, GRNs, with gene annotations derived from domain knowledge, wherein the vertices of the multimodal knowledge graph are genes, and the gene annotations enrich relationships between the genes, and by clustering gene embeddings of the GRNs. multiple We create a gene module,GM,of the sequencing data samples, and embed the samples in a multimodal knowledge graph.,For each sample of the sequencing data, each sample is embedded, multiple generating an activation vector represented as a distance between the centroid of each of the GMs;

[0008] According to an embodiment of the present invention, it has been recognized that the above-mentioned objectives can be achieved by a method and system that uses an algorithm to create a low-level representation of noisy cellular expression obtained from high-throughput sequencing (e.g., single RNA cell sequencing). Specifically, an embodiment of the present invention relates to a method for representing each cell in the data as the distance between its embedding and the centroid of a gene cluster (gene module), obtained, for example, by appropriately clustering multimodal information derived from domain knowledge. The method may also include specific functions for imputing missing values ​​from raw high-throughput data. One important application resulting from this representation is patient stratification. By using activation vectors as input for machine learning algorithms, it is possible to predict the response of cells (and therefore patients) to specific therapies by using the activity levels of gene modules as highly discriminatory features.

[0009] According to an embodiment of the present invention, a key aspect of the method consists of constructing an activation vector (AV). The construction process can start with a gene regulatory network (GRN), a knowledge graph whose vertices are genes and associated with additional features that can be obtained from domain knowledge. Gene modules (GMs) can be obtained by clustering the embeddings of the vertices of the GRN. By linearly combining high-throughput sample data with gene embeddings, a single sample can be represented in an embedded form. Finally, since the GM clusters and samples are in the same embedding space, it is possible to calculate the distance of each sample (or cell, in the case of sc-seq) from each centroid of the different GMs, thus obtaining an AV representation.

[0010] According to one embodiment of a processing system for quantifying cellular activity from high-throughput sequencing data, the present invention proposes a multi-stage computational pipeline for patient stratification based on scRNA-seq data. The pipeline first combines multiple modalities of domain knowledge, including gene regulatory networks and gene annotations, and uses a graph convolutional neural network (GCNN) to learn low-dimensional representations of genes. Genes are further clustered into functionally related gene modules. Second, a collaborative filtering approach is used to remove artifacts and other noise sources from the scRNA-seq data. Next, the cleaned scRNA-seq data is combined with the low-dimensional gene module representations to create low-dimensional representations for each cell. A supervised machine learning model is trained to predict whether each cell comes from a sample that will respond to a particular treatment. Finally, these predictions are combined to stratify samples into responders and non-responders.

[0011] According to embodiments of the present invention, the sequencing data may include single RNA cell sequencing data, whereby each sample of the sequencing data is a single cell.

[0012] According to embodiments of the present invention, embedding samples of sequencing data into a multimodal knowledge graph may include linearly combining the sequencing data with gene embeddings.

[0013] According to an embodiment of the present invention, the method may further include using a graph convolutional neural network (GCNN) to create an embedding of each gene based on the GRN with gene annotations. The embeddings may then be clustered and provided as a GM.

[0014] According to an embodiment of the present invention, the method may further comprise applying a collaborative filtering algorithm to remove dropout values ​​and other noise sources from the high-throughput sequencing data before generating the activation vector.

[0015] According to an embodiment of the present invention, the method may further comprise using the activation vector as input for training a machine learning algorithm to predict the response of individual cells to the applied drug.

[0016] According to embodiments of the present invention, the method may further comprise using a trained machine learning algorithm to predict the label "responder" or "non-responder" for each cell of the sequencing data. In this context, a voting process may be performed that uses confidence scores to establish the overall response of the patient to a particular drug, where the proportion of cells predicted to respond to the drug represents the confidence that the patient will respond to the drug.

[0017] According to an embodiment of the present invention, the method may further comprise the step of providing as output one or more of the most significant activation vectors used in the prediction as an explanation for the obtained results.

[0018] According to embodiments of the present invention, the processing system can be integrated into a platform to be used, for example, within a hospital or insurance company. Alternatively, the processing system can be used within bioinformatics pipelines and products. In either case, before it can be operational, the processing system requires high-throughput training data with known responses to each drug of interest. Therefore, it cannot be effective in purely prospective studies where the outcome is unknown. Furthermore, according to embodiments, external domain knowledge is required for the generation of valid gene modules.

[0019] There are several ways in which the teaching of the present invention can be advantageously designed and further developed. To this end, reference is made on the one hand to the dependent claims and on the other hand to the following description of preferred embodiments of the invention, illustrated by the figures, by way of example. In connection with the description of preferred embodiments of the invention with the aid of the figures, generally preferred embodiments and further developments of the present teaching are described. [Brief explanation of the drawings]

[0020] [Figure 1] 1 is a flowchart for obtaining activation vector representations for patent stratification according to one embodiment of the present invention. [Figure 2] FIG. 1 is a schematic diagram illustrating the identification of gene clusters from extrinsic information embedding according to one embodiment of the present invention. [Figure 3] FIG. 1 is a schematic diagram showing the imputation of dropout expression values ​​according to one embodiment of the present invention. [Figure 4] FIG. 1 is a schematic diagram illustrating representations of individual cells clustered through hierarchical clustering according to one embodiment of the present invention. [Figure 5] FIG. 1 is a schematic diagram illustrating the architecture of a system for performing patient stratification from single-cell data according to an embodiment of the present invention. [Figure 6] FIG. 1 shows prediction results (true positive rate vs. false positive rate) obtained by a method according to an embodiment of the present invention compared to a prior art method. [Figure 7] FIG. 1 illustrates gene module clustering according to one embodiment of the present invention. [Figure 8] FIG. 1 illustrates imputation benchmarking by comparing imputed values ​​of two replicate cells according to one embodiment of the present invention. [Figure 9] FIG. 1 illustrates hierarchical clustering of activation vectors according to one embodiment of the present invention. [Figure 10] FIG. 10 is a diagram showing activation vector distribution by a gene module according to one embodiment of the present invention. [Figure 11]FIG. 1 shows gene expression distributions for top genes according to one embodiment of the present invention. [Figure 12] FIG. 1 illustrates a pipeline for a patient stratification process according to one embodiment of the present invention. [Figure 13] FIG. 1 shows a box plot indicating gene module importance for patient stratification according to one embodiment of the present invention. [Figure 14a] 1 shows a table illustrating gene ontologies contained within different clusters according to one embodiment of the present invention. [Figure 14b] 1 shows a table illustrating gene ontologies contained within different clusters according to one embodiment of the present invention. [Figure 14c] 1 shows a table illustrating gene ontologies contained within different clusters according to one embodiment of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0021] Tumor heterogeneity poses a serious challenge in cancer treatment. Specifically, high intratumor heterogeneity can generate subpopulations of cells with very different gene signatures within the tumor microenvironment (TME). As a result, even patients with seemingly similar conditions can respond very differently to the same treatment. Therefore, it is crucial to analyze the specific gene signatures that characterize tumor heterogeneity itself to understand which patients will respond to a particular treatment. In the context of this disclosure, this grouping of patients into responders and non-responders is referred to as patient stratification.

[0022] High-throughput, single-cell sequencing methods, such as single-cell RNA sequencing (scRNA-seq), enable detailed analysis of entire gene profiles at single-cell resolution. Therefore, scRNA-seq has become a useful technique for investigating the heterogeneity within the TME. However, analyzing scRNA-seq data is challenging due to various artifacts, such as batch effects, dropout events, technical noise, and statistical challenges such as the curse of dimensionality. In particular, dropout effects and technical noise are very frequent, both of which result in the observed expression of a particular transcript being zero. This makes it difficult to distinguish whether an observation is due to non-expression of the transcript or non-detection.

[0023] Various imputation methods have been developed to account for dropout and other artifacts. Conceptually, these algorithms borrow information across similar transcripts and cells to fill in missing values. Meanwhile, MAGIC (see, for reference, D. van Dijk et al.: Recovering gene interactions from single-cell data using data diffusion. Cell 174(3), pp. 716-72927 (2018)) and SAVER (see, for reference, M. Huang et al.: SAVER: gene expression recovery for single-cell RNA sequencing. Nature Methods 15, pp. 539-542 (2018)) use Markov transition matrices and gene-gene relationships, respectively, to remap all values ​​(not just missing ones). Another approach, called scImpute (see WV Li et al.: An accurate and robust imputation method scImpute for single-cell RNA-seq data. Nature Communications 9(997) (2018) for reference), relies on mixture models to predict the most likely dropouts and correct only those. Other techniques rely on deep learning (DL) algorithms to perform their tasks.DCA (for reference, see G. Eraslan et al.: Single-cell RNA-seq denoising using a deep count autoencoder. Nature Communications 10(390) (2019)) uses an autoencoder to capture nonlinear relationships between genes, resulting in more accurate imputation, while GraphSCI (for reference, see J. Rao et al.: Imputing Single-cell RNA-seq data by combining Graph Convolution and Autoencoder Neural Networks (2020). bioRxiv 2020.02.05.935296) uses gene-gene networks and the expression values ​​extracted from them to perform imputation by using graph convolutional networks and autoencoders. While such imputation techniques are important for preparing data for further analysis, they do not target patient stratification issues.

[0024] Another notable application of deep learning (DL) is for unsupervised analysis. Autoencoders (see, for example, X. Li et al.: "Deep learning enables accurate clustering with batch effect removal in single-cell RNA-seq analysis," Nature Communications 11(2338) (2020)) have been used to perform dimensionality reduction and batch effect noise reduction to obtain meaningful clusters by accounting for nonlinearities in the data, with the aim of discovering specific gene signatures and eliminating batch effect noise. Other studies have used DL approaches to analyze intratumor heterogeneity. For example, Sharma et al.: "Longitudinal single-cell RNA sequencing of patient-derived primary cells reveals drug-induced infidelity in stem cell hierarchy." Nature Communications 9(4931) (2018) identifies cancer-resistant subpopulations by characterizing different tumor progression stages. However, unsupervised analysis approaches have limitations in their use for patient stratification because they do not aim to predict outcomes for unseen cases.

[0025] Some strategies have focused on characterizing individual cells by directly incorporating gene sets rather than traditional differential expression of individual genes. For example, PAGODA (see J. Fan, J. et al.: Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nature Methods 13, pp. 241–244 (2016) for reference) uses gene sets based on gene ontology (GO) and other annotated pathways, as well as some newly identified ones from scRNA-seq expression. It then characterizes cells according to gene sets with overdispersed expression values. DCell (see J. Ma et al.: Using deep learning to model the hierarchical structure and function of a cell. Nature Methods 15, pp. 290–298 (2018) for reference) builds a neural network structure based on hierarchical relationships between annotations, such as GO, specifically tailored for yeast. The network was then trained to predict cell viability based on single- and double-gene deletion genotype data. Expression of a set of "gene modules" was used to hierarchically cluster tumor samples (see F.G. Frost et al.: Pan-cancer RNA-seq data stratifies tumors by some hallmarks of cancer. Journal of Cellular and Molecular Medicine 24, 418-430 (2020) for reference). However, this approach was not designed to stratify unseen samples and cannot make de novo predictions.

[0026] In recent years, various approaches have been applied to the patient stratification problem. While much of this work has focused on general electronic health records (for a recent review, see B. Shickel et al.: Deep EHR: A survey of recent advances in deep learning techniques for electronic health record (EHR) analysis. IEEE Journal of Biomedical and Health Informatics 22(5), pp. 1589–1605 (2018)), other approaches have also been applied. For example, D. Toro-Dominguez et al.: Differential treatments based on drug-induced gene expression signatures and longitudinal systemic lupus erythematosus stratification. Scientific Reports 9(15502) (2019) devised a patient stratification approach based on the correlation of clinical indicators with neutrophil counts over several time points. Frohlich et al.: Premenopausal breast cancer: potential clinical utility of a multi-omics based machine learning approach for patient stratification. EPMA Journal 9, pp. 175-186 (2018) used gradient boosting to stratify breast cancer patients into risk groups from proteomic and metabolomic data, while other machine learning-based patient stratification approaches relied on feature selection approaches (as described in P. Jeyananthan et al.: Classification and regression analysis of lung tumors from multi-level gene expression data. Proceedings of the International Joint Conference on Neural Networks (2019)).Clustering-based approaches (see, for example, Y. Jang et al.: CaPSSA: visual evaluation of cancer biomarker genes for patient stratification and survival analysis using mutation and expression data. Bioinformatics 35(24), pp. 5341-5343 (2019)) and statistical model-based differential expression analysis (see, for example, F.J. Campos-Laborie et al.: DECO: decompose heterogeneous population cohorts for patient stratification and discovery of sample biomarkers using omic data profiling. Bioinformatics 35(19), pp. 3651-3662 (2019)) have also been used for patient stratification. However, all of these approaches base their stratification on the expression of a single gene (or protein, metabolite, etc.). Despite significant advances in sequencing and other technologies, reproducibility across platforms and study sites remains low (for reference, see L. Yu: RNA-Seq reproducibility assessment of the sequencing quality control project. Cancer informatics 19, pp. 1-5 (2020)). Therefore, the generalizability of methods that rely on such individual measurements is unclear.

[0027] According to another approach (see, for reference, K.-T. Kim et al.: Single-cell mRNA sequencing identifies subclonal heterogeneity in anti-cancer drug responses of lung adenocarcinoma cells. Genome Biology 16(127) (2015)), a risk score was derived to predict drug response in lung adenocarcinoma (LUAD) based on the expression of a known panel of genes. In the context of this study, it was also shown that the risk score could be meaningfully applied to individual cells. According to the inventors' assessment, this is the approach most relevant to the proposed patient stratification strategy described in this disclosure. As further explained below, the patient stratification solution of the present invention is experimentally compared to similar approaches and demonstrated to be superior.

[0028] Embodiments of the present invention provide a computational pipeline based on a graph neural network approach that combines scRNA-seq expression data with prior knowledge to perform patient stratification. In contrast to typical differential expression analysis, which focuses on the fold changes of specific genes, embodiments of the present invention achieve patient stratification by evaluating the activity levels of gene modules across all of their cells. Embodiments of the present invention implement a graph neural network-based approach that combines gene regulatory network (GRN) information with gene annotations into a single vector representation, or embedding. The embeddings can then be clustered to obtain gene modules (GMs), which, in the context of this disclosure, are understood to represent groups of similar genes. According to further embodiments, the present invention provides a collaborative filtering approach for removing artifacts and noise in scRNA-seq expression data.

[0029] Generally, therefore, embodiments of the present invention can provide an interpretable patient stratification approach that takes tumor heterogeneity into account using gene modules and cleaned scRNA-seq expression data. Embodiments of the present invention outperform existing approaches. Specifically, as experimentally evaluated by the inventors and as shown below, embodiments of the present invention significantly outperform baselines for stratifying multiple myeloma patients into responders and non-responders to the drug panobinostat.

[0030] According to one embodiment, the present invention provides a method for quantifying cellular activity from high-throughput sequencing data, the method comprising the following steps / aspects: 1. Utilizing a multimodal knowledge base by combining gene regulatory networks, GRNs (graphs whose vertices are genes) with annotations derived from domain knowledge. 2. Embedding process of GRNs and clustering of their vertices to obtain gene modules and their respective cluster centroids. 3. Acquisition of high-throughput data in the form of a sample x gene matrix, where in the case of sc-seq, each cell corresponds to a sample. 4. Combining the sequencing data with the gene embeddings (from step 2) to produce an embedding of the sample. 5. Creation of activation vectors based on the distance between the sample and the gene module centroid. In other words, a distance-based representation in which the distance between the embedded gene clusters (a measure of genetic functionality) and the cell's embedding (a measure of genetic activity) is used to create an activation vector, i.e., a representation of the patient's cell as the activation level of the gene clusters.

[0031] According to one embodiment, the present invention provides a computational pipeline 100 for patient stratification that consists of three main stages / components. A complete flowchart describing the process is outlined in Figure 1. Figures 2-4 show the individual stages / components in more detail.

[0032] First, in step S110 in FIG. 1 and in further detail in FIG. 2, the patient stratification pipeline 100 uses external knowledge from gene regulatory networks (GRNs) and gene annotations to create gene modules (GMs) using a graph neural network. The GMs are created in such a way that each GM contains a group of similar genes. Second, in step S120 in FIG. 1 and in further detail in FIG. 3, the patient stratification pipeline 100 uses bioinformatics and collaborative filtering methods (e.g., singular value decomposition (SVD) filtering) to remove artifacts and noise from the scRNA-seq data. Third, in step S130 in FIG. 1 and in further detail in FIG. 4, the patient stratification pipeline 100 creates a GM activation vector for each cell based on inputs from steps S110 and S120 and trains a machine learning model to predict drug responses for individual cells. Finally, these individual cell predictions are combined to stratify patients into responders and non-responders. The following sections describe each step in detail.

[0033] Learning Gene Embeddings In one embodiment, gene embedding training is performed by using a graph convolutional neural network (GCNN) to create a dense vector representation, or embedding, of each gene based on the GRN and gene annotations. The embeddings are then clustered, and these clusters are adopted as gene modules, GM.

[0034] Each gene in the GRN can be associated with a set of gene annotations, which are treated as bag-of-word features. According to one embodiment, the set of gene annotations can include, for example, a gene ontology vocabulary. The representations are then learned by using the embedding propagation (EP) algorithm, for example, as described in Garcia-Duran, A., Niepert, M.: Learning graph representations with embedding propagation. Advances in Neural Information Processing Systems 30 (2017), the entire contents of which are incorporated herein by reference. Briefly, EP is a GCNN that learns a feature-by-feature embedding for a vertex by minimizing the difference between its embedding and that obtained from neighboring vertices in the GRN.

[0035] More formally, the EP learning framework involves two stages. First, the EP algorithm can learn embeddings for each annotation and for each gene discrimination by propagating messages along the edges of the GRN. In the second step, the EP algorithm can learn gene representations by combining the learned annotations and discrimination embeddings. That is, the gene representations combine information from the GRN structure and annotations.

[0036] The message is a coding function f that maps annotations and gene identities to vector representations. θ The encoding function is parameterized by θ and must be differentiable. Specifically, in one embodiment of the present invention, a neural network consisting of a single fully connected layer, resulting in an embedded lookup table, may be used. In the following, the annotation and identification of a gene g will be referred to as a(g), and the neighbors of g in the GRN will be referred to as n(g). The notation a(n(g)) will be used to denote the annotation and identification of all neighbors of g.

[0037] Furthermore, to represent the current embedding of g, we use h(g)=agg(f θ (x)|x∈a(g)) is used, where agg is an aggregation function, such as choosing the element-wise mean or maximum from a set of embeddings. Similarly, to indicate the aggregated embedding of all neighbors of g in the GRN,

[0038]

number

[0039] is used.

[0040] Mathematically, EP minimizes the following loss function:

[0041]

number

[0042] where d is the Euclidean distance and [x] + is the positive part of x, and γ>0 is the margin hyperparameter.

[0043] Note that in practice, evaluating the inner sum may not be feasible, and thus negative example sampling may be used (e.g., as described in Kotnis, B., Nastase, V.: Analysis of the impact of negative sampling on link prediction in knowledge graphs. KnoProceedings of the 1st Workshop on Knowledge Base Construction, Reasoning and Mining (2018), the entire contents of which are incorporated herein by reference), and a single random gene other than g may be evaluated. θ The parameters of can be updated using a standard backpropagation algorithm until convergence.

[0044] For the remainder of the patient stratification pipeline 100 shown in Figure 1, h(g) is used as the embedding for the gene. Note that the embeddings depend only on the GRN and annotations, and therefore, they can be used for multiple scRNA-seq datasets.

[0045] Creating a gene module According to one embodiment of the present invention, to obtain gene modules, gene embeddings can be clustered using a Bayesian Gaussian Mixture Model (BGMM). For example, cluster centroids can be taken as summaries of gene modules. That is, each gene module is represented by the centroid of its respective cluster. Herein, the centroid for gene module m is denoted as C.

[0046] As a preprocessing step and for numerical stability, the embedding can be scaled so that the values ​​in each dimension have a mean of 0 and a variance of 1. In the context of their analysis, the inventors found that 10 gene modules provided a good tradeoff between accuracy, computational demands, and granularity for interpretation.

[0047] Expression data preparation According to embodiments of the present invention, the patient stratification pipeline 100 may include preparing scRNA-seq for downstream analysis. In this context, first, standard procedures (e.g., as described in "Integrating single-cell transcriptomic data across different conditions, technologies, and species," Nature Biotechnology 36(5), pp. 411-420 (2018), the entire contents of which are incorporated herein by reference) may be used to normalize the number of reads observed per cell:

[0048]

number

[0049] where c is a specific cell, g is a gene (or transcript, etc.), and x c,g is the observed read count of gene g for cell c, and x c,g is the normalized expression value.

[0050] Next, to account for dropout and other technical artifacts and noise, a standard collaborative filtering approach can be used, such as singular value decomposition (SVD) as implemented in the Surprise package (Hug, N.: Surprise: A python library for recommender systems. Journal of Open Source Software 5, 2174 (2020), the entire contents of which are incorporated herein by reference). Intuitively, SVD estimates the expression of cell-gene pairs based on similar cells and genes. Specifically, it solves the following optimization problem to minimize the error between the estimated and observed values:

[0051]

number

[0052] where:

[0053]

number

[0054] and

[0055]

number

[0056] are the observed and predicted (normalized) expression values, respectively, and q c and p g are their coefficients, and b c and b g is the bias. Set R train contains all non-zero values ​​from the sc-seq data matrix. The regularization term λ is a hyperparameter of the model, while the coefficients and biases are learned using an iterative stochastic gradient descent optimization algorithm.

[0057] For example, experimentally, 40 epochs may be chosen while leaving all other parameters at their default values.

[0058] After the model parameters are learned, the expression values ​​are

[0059]

number

[0060] Importantly, all expression values ​​for all cells are updated, and therefore, approaches according to embodiments of the present invention as disclosed herein address not only dropouts but also other sources of noise.

[0061] As a final step, the total expression in each cell can be scaled so that it sums to 1, as shown in (4):

[0062]

number

[0063] As will be appreciated by those skilled in the art, it is equally possible to employ state-of-the-art imputation methods other than the SVD approach described above, but the inventors found that SVD was competitive and generally outperformed other methods in their analysis, and therefore SVD was used for the remainder of the analysis.

[0064] Creation of cell activation vectors In the context of this disclosure, an activation vector (AV) is a representation of a single cell in terms of the activity level of its GM. According to one embodiment of the present invention, an AV can be created in two steps. First, an embedding can be created for each cell based on gene embeddings learned from the GCNN. Second, the distance between the GM centroid and the cell embedding gives the AV for that cell.

[0065] Cell embedding In one embodiment, the cell-wise embedding h(c) can be obtained by linearly combining the gene embeddings learned from GCNN with the gene expression values ​​estimated by collaborative filtering, as shown in (5):

[0066]

number

[0067] The obtained cell embedding values ​​can be scaled by using the same scaler that was applied to the gene embedding matrix before clustering to define the GM.

[0068] Cell Activation Vector The final cell AV can be taken as the cosine similarity (cos) between each gene module centroid and the cell embedding: AV c =[cos(h(c),C1),...,cos(h(c),C M )], (6) Here, M represents the number of gene modules (e.g., M=10, as outlined above).

[0069] Patient stratification The final step of the patient stratification pipeline 100 according to the embodiment of the present invention shown in Figures 1-4 involves training a machine learning model to predict, for a particular drug, which cells will result from responding samples and which will result from non-responding samples. Typically, it is not possible to detect whether a single cell will respond to a drug using scRNA-seq. Therefore, to label individual cells as responders and non-responders, it is first necessary to determine whether each sample responds to the drug. For example, the entire sample may be divided into two parts, one part treated with the drug and the other part used for single-cell sequencing. Based on the response of the part of the sample treated with the drug, the sequenced cells can be labeled as responders or non-responders. Importantly, the sequenced cells were not treated with the drug. Therefore, the label indicates whether the cells will respond to the treatment, not whether they have already responded to the treatment.

[0070] More specifically, this approach constructs a training set based on a cohort of untreated samples. The training set contains the AV for each cell from the samples in this cohort, and each cell's (binary) label is whether it comes from a sample that responded to treatment. A classifier can then be trained to predict whether a cell comes from a responder sample. The trained classifier can then be used to classify all cells from a separate cohort of untreated samples (used as a test set) as either coming from a responder or not. Each sample from this test cohort can then be stratified into a responder or non-responder using a simple voting scheme based on the prediction of each cell from that sample.

[0071] An embodiment of the present invention uses a method algorithm to represent data obtained from high-throughput sequencing in the form of activation vectors, a representation in which each sample is expressed as the activity level of a functional cluster of genes (GMs) contained within it.

[0072] GMs can be obtained by creating multimodal gene regulatory networks in which the relationships between genes are enriched with additional annotations derived from domain knowledge.

[0073] An important aspect of embodiments of the present invention resides in a method for assessing cluster activity, which is the distance between the embedding of each sample from high-throughput data and the centroid of the GM. This assessment is possible thanks to the fact that both sample and gene clusters lie in the same embedding space, thus allowing the distance between the two vectors to be assessed.

[0074] It is important to remember that, according to embodiments, gene modules are obtained solely from domain knowledge, which contributes to improving the explainability of results, for example, when activation vectors are used in machine learning contexts.

[0075] 5 illustrates one embodiment of the algorithm, along with a use case in which the method for quantifying cellular activity as disclosed herein performs the core activity and activation vectors are used for patient stratification from single-cell data. Generally, patient stratification approaches are used to classify patients as responders or non-responders to a particular treatment. Predicting such responses in silico offers several advantages for both patients and clinicians, as it avoids in vivo testing of different drugs for treatment. In the context of this disclosure, it is shown how this process can be applied to a single patient by testing the response of high-throughput sequencing samples with different drugs.

[0076] As shown in FIG. 5, the embodiment includes: multiple interacting components, namely:

[0077] 1. A database 510, located, for example, in a hospital, containing high-throughput expression data from patients and specially cultured cells (cell lines) for different tumors. This database 510 can be configured to contain expression from treated and control samples to have labels for prediction. In addition, it can be configured to contain domain knowledge acquired from different sources, such as annotations and gene regulatory networks.

[0078] 2. A server 520 running an algorithm that represents data as activation vectors and performs a stratification process through machine learning algorithms. Depending on the treatment of interest, data for that particular treatment is collected from a database and used to create an activation vector representation of the data provided by client 540 (see below).

[0079] 3. A machine 530 for performing high-throughput data analysis from samples obtained from patients.

[0080] 4. A client 540 used by the clinical facility to receive the high throughput data from the machine 530 and send them to the server along with a list of drugs for which predictions are desired.

[0081] 5. A therapy plan scheduler 550 that uses the results of the stratification performed in the server 520. It may be configured to receive the identification of the medications to which the patient will respond, and may be used by the clinician to calculate dosages.

[0082] Database 510 The database 510 may be located within the hospital or, if compliant with privacy regulations, may be accessed through a cloud service. In one embodiment, the database 510 contains high-throughput samples of different tumor cells treated with different drugs. For each sample, there will be pre- and post-treatment data. The samples will have multiple labels related to their response to different drugs. The database 510 will also contain external domain knowledge in the form of gene-gene interaction networks (gene regulatory networks, GRNs) and functional annotations for each gene. Additional information about the genes, such as free text, may also be available. These data will be obtained from the high-throughput analysis performed by the machine 530 and used to train an algorithm contained within the server 520 to predict the response of samples uploaded by the client 540.

[0083] Shown below is a snapshot of the high-throughput data contained within database 510, where each row is a cell and each column is the expression level of the gene within it. Each sample is also annotated with a label indicating whether it came from a sample that responded to each drug. Database 510 may store only cells from pre-treatment samples.

[0084] [Table 1]

[0085] [Table 2]

[0086] Server 520 The server 520 contains the activation vector algorithm, i.e., the core method that takes high-throughput data from patients as input and returns a prediction regarding the response for a drug of interest. The specific formula and technical implementation of this algorithm may follow the principles detailed above.

[0087] According to one embodiment, the workflow of the algorithm may have offline and online stages as follows:

[0088] The offline phase may include the following steps:

[0089] Gene module and cluster centroid creation: 1. Merging together the domain knowledge (annotations, GRNs) related to the gene of interest, thus creating a multimodal graph. 2. Embed and cluster the vertices of this graph into a GM. 3. Save the cluster centroids, cluster designations, and gene embeddings.

[0090] Training a drug response prediction model: 1. Construct an AV from the high-throughput sample dataset in database 510. a. As an optional step, dropout values ​​of the high-throughput data on the database may be imputed through a collaborative filtering algorithm. b. Rescale every row (cell) of data values ​​so that they sum across all genes to 1. Linearly combine the rescaled data and the embedded gene vertices to obtain an embedding for each cell. c. Evaluate the distance between each cell and each GM centroid, which results in each cell being represented as an activation vector. 2. Train an ML classifier to predict the label of each cell. 3. Repeat processes 1-2 for each drug of interest.

[0091] The online phase may include the following steps (data for this online phase may come from high-throughput data acquisition performed by machine 530 provided via client component 540): 1. As an optional step, dropout values ​​of high-throughput data on the database may be imputed through a collaborative filtering algorithm. 2. Each row of data values ​​(samples) can be rescaled so that they sum to 1 across all genes. The rescaled data and the embedded gene vertices are linearly combined to obtain an embedding for each sample. 3. Evaluate the distance between each cell and each GM centroid, which results in each sample being represented as an activation vector. 4. Predict the label (responder / non-responder) for each sample using properly trained ML. 5. A voting process uses a confidence score to establish the patient's overall response to that particular drug. Specifically, the percentage of samples (cells) predicted to respond provides confidence that the patient will respond to the drug. Also, according to one embodiment, it may be provided to return the most significant activation vectors used in the prediction to give the clinician additional explanation for the results.

[0092] Predictions of medications that will ensure response are returned to the client 540 and communicated to the clinician within the therapy scheduler 550 .

[0093] High-Throughput Data Acquisition530 This can be any machine that takes a sample from a patient as an input and performs high-throughput sequencing (e.g., single-cell sequencing).

[0094] Client 540 The client 540 can be used by a clinician to upload sequencing data to the server 520 to perform in silico predictions of response. The server 520 can be configured to return a response prediction to the client 540, which can be configured to forward the response to the therapy scheduler 550. The clinician can use the therapy scheduler 550 to visualize the selected treatment.

[0095] Treatment Scheduler 550 The scheduler 550 may be implemented as a front-end application that a clinician can use to load the predicted results and begin scheduling therapy with the selected treatment. Specifically, the scheduler 550 may be configured to automatically generate the appropriate prescription for the patient and, in some cases, forward it to a pharmacy.

[0096] According to a slightly different embodiment, the present invention provides a computer-implemented method and processing system for quantifying cellular activity from high-throughput sequencing data, in which activation vectors are used for patient stratification from bulk sequencing data.

[0097] While the approach of the above-described embodiment is based on single-cell sequencing data, where each biological sample yields many observations (one per cell), the approach is also applicable to bulk sequencing, where each biological sample yields a single observation (such as an average across all cells in the sample). The GM and AV approaches can also be used to create sample embeddings for bulk sequencing data. The only modification required is to step 5 (i.e., the "voting" step) of the online phase (see above), performed by server 520. Because each bulk sequencing sample consists of only a single observation, there is no opportunity for voting. Instead, server 520 can be configured to directly use predictions from the trained ML model within step 4 of the online phase (see above) and use confidence scores to calculate the patient's overall response to that particular drug.

[0098] As already outlined above, currently standard sequencing analysis, such as bulk RNA sequencing, allows for obtaining single reads for various transcripts / genes within a single patient. All of the above-mentioned prior art approaches base their stratification on the expression of a single gene (or protein, metabolite), and therefore require a large number of patients for their models to be valid. In other words, to perform any kind of ML-based prediction (such as patient stratification), a large number of patients must be available, which is nearly impossible. In contrast, high-throughput single-cell sequencing data allows for the acquisition of a large number of samples per patient, so embodiments of the present invention instead rely on a small number of patients. In this case, reads of the order of thousands of cells can be obtained for a single patient. This allows for the possibility of performing predictions using a small number of samples.

[0099] Furthermore, embodiments of the present invention create representations of data that look at gene signatures (gene modules), as opposed to standard techniques that look at the activity differences of single genes. This results in more accurate predictions, as can be seen from Figure 6, which shows the AUC (area under the curve) for predicting cellular responses within a patient stratification process using sc-seq data. When the AV-based method of the present invention is not used (see curves a) "chance," d) "SVD (singular value decomposition) only," and e) "as is") and therefore relies only on expression values ​​from a single gene, the results are very mediocre compared to curves c) "AV only" and d) "full," i.e., SVD+AV. Overall, the approach according to embodiments of the present invention makes ML-based predictions explainable. This is due to the fact that the predicted response within a cell can be decomposed into the contribution of a single activation vector, which has real-world meaning thanks to the domain knowledge used during its construction.

[0100] According to an application scenario, embodiments of the present invention focus on the stratification of patients affected by multiple myeloma (MM). MM is one of the most common forms of hematopoietic malignancies, characterized by the clonal proliferation of malignant B-cell-derived plasma cells within the bone marrow (see Rajkumar, SV: Multiple myeloma: 2012 update on diagnosis, risk-stratification, and management. American Journal of Hematology 87, pp. 77-88 (2012)). Currently, proteasome inhibitors (PIs) are the most effective therapy for MM. Among the most important, bortezomib (Velcade) was the first PI approved for the treatment of MM (in 2003), followed by carfilzomib (Kyprolis) (approved in 2012) and oral PIs including ixazomib (Ninlaro, 2015) and panobinostat (Farydak). In this context, experiments were carried out to perform stratification of patients treated with panobinostat by looking at their gene signatures derived from single-cell sequencing techniques.

[0101] For this study, scRNA-seq expression data from four multiple myeloma patients treated with panobinostat were used (samples collected at BioQuant, Heidelberg, DE). Primary cells from the patient samples were sorted by CD138 to isolate B cells. These cells were treated with either control (DMSO) or panobinostat and ultimately subjected to 5'-end scRNA-seq. A total of 29,997 cells and 26,063 distinct transcripts were quantified from the patient samples. In addition, three cell line models for MM were treated with either DMSO or panobinostat and subsequently sequenced. A total of 6,478 cells and 11,723 transcripts were quantified from the cell lines.

[0102] Of the patient samples, domain experts qualitatively determined that three samples responded to the drug, while one sample did not. Of the cell lines, two responded to the drug and one did not. The DMSO-treated samples were considered to be indicative of the pre-treatment samples, and they form the basis of the patient stratification algorithm as described in this disclosure. The panobinostat-treated samples were not used except to determine the label for the pre-treatment cells. Table 1 details the seven samples used in this study.

[0103] [Table 3]

[0104] As external knowledge, a B cell-specific gene regulatory network was used to obtain information about gene-gene interactions. Gene Ontology (GO) (Generic GO Slim Subset, see http: / / current.geneontology.org / ontology / subsets / goslim_generic.obo) annotations were collected for 18,177 genes. Those genes for which no GO annotation was available were removed from the GRN.

[0105] The top 2000 most expressed genes, based on the number of cells in which they were detected, were retained as the gene set used for the analysis. Next, all genes not contained within a GRN were removed. In total, the final network and analysis were based on 1713 genes. Table 2 provides summary information on the connectivity of the GRNs and GO annotations.

[0106] [Table 4]

[0107] result Genetic implants and gene modules We first qualitatively evaluated the gene embeddings and resulting gene modules (GMs). Figure 7 shows a 2D UMAP (see McInnes, L., Healy, J., Melville, J.: UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (2018). arXiv:1802.03426 [stat.ML]) projection of the embeddings. Clustering of the latent gene representations is based on gene regulatory networks and gene ontology annotations. Each cluster is a gene module.

[0108] Furthermore, we performed an enrichment test (hypergeometric test with Bonferroni correction, p<0.05 was taken as significant) to detect overrepresented gene ontology (GO) annotations within each GM. The GOs encompassed within the clusters are shown in the tables shown in Figure 14a-c. Specifically, we considered only the GO Slim (Generic GO Slim Subset, see http: / / current.geneontology.org / ontology / subsets / goslim_generic.obo) set of annotations.

[0109] Importantly, by manually inspecting the enriched GO vocabulary for each GM, we find that they form biologically meaningful groups. For example, GM4 represents proteosomal activity, while GM1 could be categorized as nuclear activity. This analysis supports our hypothesis that GMs meaningfully combine structured domain knowledge from GRNs and annotations into a single, coherent representation.

[0110] Attribution Next, we validated the performance of our collaborative filtering-based imputation approach. Specifically, we compared our proposed approach with several state-of-the-art imputation approaches, including MAGIC (D. van Dijk et al.: Recovering gene interactions from single-cell data using data diffusion. Cell 174(3), pp. 716-72927, 2018), SAVER (M. Huang et al.: SAVER: gene expression recovery for single-cell RNA sequencing. Nature Methods 15, pp. 539-542, 2018), and scImpute (W.V. Li et al.: An accurate and robust imputation method scImpute for single-cell RNA-seq data. Nature Communications 9(997), 2018). For comparison, we used the GSE45719 dataset (mouse preimplantation embryos at 10 different developmental stages), Series GSE45719 (see https: / / www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE45719). As an evaluation metric, we used the Pearson correlation coefficient (PCC).

[0111] To analyze the accuracy of SVD reconstruction as a function of the amount of dropout (i.e., missing values) in the gene expression matrix, we performed a sensitivity test. Specifically, we downsampled the number of observed genes per cell. The process was repeated with different observation rates. Model performance was evaluated using the root mean squared error (RMSE) metric. The workflow consisted of five steps: 1. Divide the cell-gene pairs into training, validation, and test sets 2. Set the observed percentage (20%, 40%, 60%, 80%, 100%) of the average number of expressed genes in each cell (in the training set). 3. For each cell in the training set, downsample the number of observed genes according to the observation rate. Downsampling within each cell was performed by choosing the number of observed genes for that cell from a standard normal distribution with a mean equal to the observation rate and a standard deviation equal to 10. Next, we dropped out all other expression values ​​from that cell. 4. Optimize the hyperparameters of the SVD model using grid search and the validation set. 5. Finally, the model was used on the test set and the RMSE was calculated again as in the previous step.

[0112] For each observation rate, we repeated the analysis five times to assess performance variation due to randomness. The validation and test sets were kept the same across all observation rates and samples to ensure comparable results for all settings. Even if RMSE is evaluated only for known entities, having a very low value for error makes reconstructed dropout more likely to be meaningful.

[0113] Figure 8 shows imputation benchmarking based on a comparison of the imputed values ​​of two replicate cells. Because the cells are from the same developmental stage, we expect the expression values ​​of most genes to be similar. Therefore, we would expect most imputed values ​​to fall along the main diagonal. As can be seen from Figure 8, our collaborative filtering-based imputation is competitive with state-of-the-art PCC: 0.92 vs. 0.75 (SAVER), 0.82 (scImpute), and 0.95 (MAGIC). More specifically, our SVD approach yields a higher PCC between imputed values ​​between the two cells than most other methods.

[0114] Activation Vector As a first qualitative result, we performed hierarchical clustering of the AV for cells from all untreated samples. Figure 9 shows the hierarchical clustering of activation vectors for all samples treated with the control. Cells from non-responder samples are labeled in light colors. As can be seen from Figure 9, cells tend to form clusters of responders and non-responders. This suggests that machine learning algorithms can be successful in predicting responders and non-responders. On the other hand, cells from individual samples do not all cluster together. Encouragingly, this suggests that the AV captures trends across samples rather than simply reflecting sample-specific artifacts.

[0115] Next, we qualitatively examined the distribution of values ​​for responders and non-responders in the pre-treatment samples for AV. As shown in Figure 10, several of the GMs have highly discriminatory distributions between responders (R) and non-responders (NR). This again suggests that the GMs capture meaningful biological activity and that the classifier can distinguish between responders and non-responders.

[0116] For comparison, we used a standard approach to identify differentially expressed genes between responder and non-responder samples. Figure 11 shows a kernel density estimate of the gene expression distribution for the top genes, i.e., the distribution of expression values ​​for the most significant genes used by the classifier. The distributions are very similar, i.e., there is no discrimination between responders (R) and non-responders (NR), and therefore, using the expression value of a single gene as a discriminant feature may not yield good predictions. This highlights the advantage of using GM for this classification task.

[0117] Machine Learning Evaluation Next, we quantitatively evaluated activation vectors (AVs) for patient stratification. For evaluation, we divided patient and cell line samples into training and testing cohorts.

[0118] We evaluated our approach using two different training scenarios based on how we split samples from patients and cell lines. In the first set of experiments (referred to as "mixed configuration" in the results), both cell line and patient samples were used in both the training and testing datasets. In the second (referred to as "CL configuration"), only cell line samples were used for training, while patient samples were used for testing. In both cases, cells from each sample appeared only in the training set or only in the testing set. This ensures that the machine learning model does not simply memorize what cells from a particular sample look like.

[0119] For both scenarios, we considered four unique computational pipelines: one in which both AV representation and imputation were used ("full"), one in which only AV was used and no imputation was used ("AV only"), one in which only imputation was used and no AV was used ("SVD only"), and a final one in which neither was used ("as is"). For all evaluations, we used 10-fold cross-validation to ensure that adequate error estimates were available.

[0120] The machine learning classifier chosen for this work was logistic regression (LR). We set two hyperparameters: tolerance for early stopping at 10e-3 and maximum number of iterations at 1000. For class weights (unweighted or balanced), coefficient penalty (l-1 or l-2), and inverse regularization strength ({1, 10, 100, 1000}), we selected the best values ​​using a grid search with cross-validation.

[0121] As described above, in splitting the dataset into training and test sets, we ensured that at least one non-responder sample and one responder sample were in both splits. As described above, in some cases, SVD was not used to impute missing gene values. In those cases, missing values ​​were imputed using the average expression of that gene across all cells in which it was observed.

[0122] Figure 12 shows how each experiment was performed: We split the training set into 10 folds using the StratifiedKFold function in sklearn (ensuring that the class ratios of the dataset are preserved within a single fold). For each split of the training set, we further split the training set to create an internal training set and an internal validation set. We then selected LR hyperparameters by training using the internal training set and evaluating using the internal validation set. The hyperparameters for which LR returned the best area under the ROC curve (AuROC) value using the internal validation set were selected as the best model for that training fold. We then performed predictions for the test set using the selected model, and obtained the AuROC value for that fold. We repeated the experiment for the remaining nine folds, and finally, we evaluated the average AUC across the 10 folds.

[0123] Cell Classification Figure 6 shows the results for both the "CL configuration" (left) and the "mixed configuration" (right). More specifically, Figure 6 shows the area under the receiver operating characteristic curve (AuROC) for each of the computational pipelines for both training scenarios. These results clearly demonstrate that the use of activation vectors significantly improves classification performance compared to using baseline methods, with or without imputed expression values. In both training scenarios, combining activation vectors with imputation produced the best results.

[0124] For the "mixed configuration," Figure 6 shows that "SVD only" and "as is" performed particularly poorly (AuROC<0.5). In these cases, the trained model consistently predicted a very high probability (>0.95) that nearly all of the cells in the test sample were from responders, and cells from non-responders were consistently predicted to have an even higher probability (>0.99) of being from responders. This phenomenon results in very low AuROC values.

[0125] Patient stratification Table 3 shows the quantitative results of patient stratification prediction.

[0126] [Table 5]

[0127] The percentages indicate the probability that a sample is predicted as a responder. The actual labels are shown in the second column. Cell line composition is shown to be more robust in predicting the correct behavior of samples. For both dataset construction scenarios, our proposed approach (collaborative filtering-based imputation and the use of AV to encode both domain knowledge and expression data) outperforms other baselines. In both scenarios, AV significantly improves results compared to using gene expression values ​​directly. Although the small number of available samples does not allow us to draw strong conclusions, this provides evidence of the value of including domain knowledge to stratify patients.

[0128] Gene Ontology and Model Interpretation The pipeline presented according to an embodiment of the present invention is interpretable by design, since all predictions are constructed using a set of features (GMs) that can be mapped to specific functions. Specifically, as described above, each GM is associated with a set of GO annotations. Furthermore, the importance of each GM in the classifier is taken as the absolute value of its corresponding coefficient in a logistic regression model. A similar "feature importance" approach (see M.T. Ribeiro et al.: "why should i trust you?": Explaining the predictions of any classifier. In: Proceedings of the 22nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 2016) could also be used for other machine learning models.

[0129] Figure 13 shows the distribution of the importance of each GM across all folds of cross-validation. Specifically, Figure 3 shows boxplots for feature importance for each of the 10 gene modules (GMs) for each of the 10 folds for (a) cell line and (b) mixed configuration. Feature importance is taken as the absolute value of the coefficient for each GM in the trained logistic regression model.

[0130] Figure 13 shows that GM #4 is consistently a significant, i.e., discriminatory, feature of the classifier. GM #4 is highly enriched in proteasomes, consistent with a study by Tanaka et al. (Tanaka, K.: The proteasome: Overview of structure and functions. Proceedings of the Japan Academy, Series B 85, ​​pp. 12-36, 2009), which showed that panobinostat plays an important role in downregulating proteasome activity. Indeed, Figure 10 shows that non-responders have low activity for GM #4 to begin with, and as a result, panobinostat cannot help these patients due to their low proteasome activity to begin with. On the other hand, responders have moderate to high activity levels for GM #4. In these cases, panobinostat can effectively downregulate proteasome activity and result in a response in patients.

[0131] In summary, according to embodiments of the present invention, a pipeline for patient stratification using the activation levels of functional gene groups within cells has been developed. This can be obtained, for example, by combining data derived from single-cell RNA sequencing (scRNA-seq) with external prior knowledge using a graph convolutional neural network. A key component of embodiments of the present invention is the creation of gene modules (GMs) based on external prior knowledge. By focusing on the activation of GMs within cells and predicting their response to treatment, analyses according to embodiments of the present invention are robust to noise from individual genes. In addition, embodiments of the present invention propose a simple collaborative filtering-based approach to account for technical and biological noise in scRNA-seq data, further improving the resilience of this approach. Experimentally, we were able to demonstrate that GM-based prediction yielded an AuROC performance of 0.97 in a set of multiple myeloma patient and cell line samples. In contrast, directly using individual gene expression values ​​yielded only 0.55, which is no better than guesswork.

[0132] Deeper analysis showed that the GM-based approach according to embodiments of the present invention can directly reveal relevant biological functions. This is in contrast to classical differential expression analysis, which only identifies individual genes. For example, in the context of the present invention, GM associated with proteasome activity was shown to be crucial for predicting response to panobinostat, which is consistent with existing literature on panobinostat.

[0133] Many modifications and other embodiments of the inventions described herein will come to mind to one skilled in the art to which these inventions pertain having the benefit of the teachings presented in the foregoing descriptions and the associated drawings. It is to be understood, therefore, that the invention is not limited to the specific embodiments disclosed, and that modifications and other embodiments are intended to be included within the scope of the appended claims. Although specific terms are employed herein, they are used in a generic and descriptive sense only and not for purposes of limitation. [Explanation of symbols]

[0134] 100 Computational Pipelines, Patient Stratification Pipelines 510 Database 520 Server 530 Machinery 540 Client 550 Therapy plan scheduler, treatment scheduler

Claims

1. A computer-implemented method, the method comprising: generating a multimodal knowledge graph by combining gene regulatory networks, GRNs, with gene annotations derived from domain knowledge, where the vertices of the multimodal knowledge graph are genes and the gene annotations enrich the relationships between the genes; creating a plurality of gene modules, GMs, by clustering the gene embeddings of the GRNs, where the embeddings are a single vector representation obtained by combining information of the GRNs with the gene annotations; Embedding samples of high-throughput sequencing data into the multimodal knowledge graph; generating an activation vector for each sample of the high-throughput sequencing data, each sample represented as a distance between the embedding and the centroid of each of the plurality of GMs; A method comprising:

2. 2. The method of claim 1, wherein the high-throughput sequencing data comprises single RNA cell sequencing data, whereby each sample in the high-throughput sequencing data is a single cell.

3. 3. The method of claim 1, wherein embedding the samples of high-throughput sequencing data into the multimodal knowledge graph comprises linearly combining the high-throughput sequencing data with the embeddings of the genes.

4. creating the embeddings for each gene based on the GRN with the gene annotations using a graph convolutional neural network, GCNN; clustering the embeddings and providing the clusters as the GM; The method of any one of claims 1 to 3, further comprising:

5. Before generating the activation vector, The method of any one of claims 1 to 4, further comprising applying a collaborative filtering algorithm to remove dropout values ​​and other sources of noise from the high-throughput sequencing data.

6. 6. The method of any one of claims 1 to 5, further comprising using said activation vectors as input for training a machine learning algorithm to predict the response of individual cells to an applied drug.

7. 7. The method of claim 6, further comprising using the trained machine learning algorithm to predict the label "responder" or "non-responder" for each cell in the high-throughput sequencing data.

8. 8. The method of claim 7, further comprising a voting process that uses a confidence score to establish an overall response of the patient to a particular drug, wherein the percentage of cells predicted to respond to the drug represents a confidence that the patient will respond to the drug.

9. 9. The method of claim 8, further comprising providing as output one or more of the most significant activation vectors used in the prediction as an explanation for the obtained results.

10. A processing system for performing the method of any one of claims 1 to 9, said system comprising one or more processors, said one or more processors: generating a multimodal knowledge graph by combining gene regulatory networks, GRNs, with gene annotations derived from domain knowledge, wherein the vertices of the multimodal knowledge graph are genes and the gene annotations enrich the relationships between the genes; creating a plurality of gene modules, GMs, by clustering the gene embeddings of the GRNs, wherein the embeddings are a single vector representation obtained by combining information of the GRNs with the gene annotations; Embedding samples of high-throughput sequencing data into the multimodal knowledge graph; generating an activation vector for each sample of the high-throughput sequencing data, each sample represented as a distance between the embedding and the centroid of each of the plurality of GMs; a processing system configured to:

11. 11. The processing system of claim 10, further comprising a database (510) containing high-throughput samples of different tumor cells treated with different drugs together with said multimodal knowledge graph.

12. 12. The processing system of claim 11, further comprising a server (520) configured to use a multimodal knowledge graph from the database (510) to train a machine learning algorithm based on the activation vectors to predict the response of individual cells to an applied drug.

13. The system further comprises a client (540) configured for use by a clinician to upload high-throughput sequencing data obtained from a patient to the server (520), the server (520) comprising: using the trained machine learning algorithm to predict the label "responder" or "non-responder" for each cell in the high-throughput sequencing data; and returning a predicted result to the client (540) including the drug that produces a response; 13. The processing system of claim 12, configured to:

14. 14. The processing system of claim 13, further comprising a treatment scheduler (550) as a front-end application configured to load the prediction results and be used by a clinician to select a treatment.

15. one or more processors of the processing system, generating a multimodal knowledge graph by combining gene regulatory networks, GRNs, with gene annotations derived from domain knowledge, wherein the vertices of the multimodal knowledge graph are genes and the gene annotations enrich the relationships between the genes; creating a plurality of gene modules, GMs, by clustering the gene embeddings of the GRNs, wherein the embeddings are a single vector representation obtained by combining information of the GRNs with the gene annotations; Embedding samples of high-throughput sequencing data into the multimodal knowledge graph; generating an activation vector for each sample of the high-throughput sequencing data, each sample represented as a distance between the embedding and the centroid of each of the plurality of GMs; A non-transitory computer-readable medium comprising code for causing

Images