Method and system for inferring cell-type abundance

Abstract

A method of inferring cell-type abundance in a biological sample, comprises: receiving bulk gene expression data of the biological sample; feeding the bulk gene expression data as an input to an artificial neural network trained to apply a non-negative matrix factorization to factorize the bulk gene expression data into a matrix multiplication of a non-negative cell type signature matrix by a non-negative cell type fraction matrix; and receiving from the network, as an output, the cell type fraction matrix, thereby inferring the cell-type abundance in the biological sample.

Description

[0001] METHOD AND SYSTEM FOR INFERRING CELL-TYPE ABUNDANCE

[0002] RELATED APPLICATION

[0003] This application claims the benefit of priority of U.S. Provisional Patent Application No. 63 / 452,353, filed on March 15, 2023, the contents of which are incorporated herein by reference in their entirety.

[0004] FIELD AND BACKGROUND OF THE INVENTION

[0005] The present invention, in some embodiments thereof, relates to computational biology and, more particularly, but not exclusively, to a method and system for inferring cell-type abundance.

[0006] In the context of cellular biology, understanding the specific types and proportions of cells present in a given sample is useful for gaining insights into tissue function, disease mechanisms, and potential therapeutic targets. Traditional methods often rely on experimental techniques like fluorescence-activated cell sorting or single-cell RNA sequencing to identify and characterize individual cell types.

[0007] Computational cell-type deconvolution leverages advanced algorithms and statistical models to infer the relative abundance of different cell types within a heterogeneous sample based on bulk gene expression data. This approach is particularly useful for complex tissues including many cell populations. The aim of the deconvolution is to decipher the gene expression profiles of individual cell types within the sample, allowing to estimate the contribution of each cell type to the overall gene expression signature of the sample.

[0008] Known computational methods rely on a signature matrix of cell-specific expression profiles to predict the cell type abundance. For example, NNLS [Charles L Lawson and Richard J Hanson. Solving least squares problems. SIAM, 1995], CIBERSORT [Newman et al., Nature methods, 12(5):453-457, 2015], CIBERSORTx [Newman etal., Nature biotechnology, 37(7):773- 782, 2019] are computational methods which are based on support vector regression and GEDIT [Nadel etal., GigaScience, 248 10(2), 022021] is a computational method which is based on linear regression, and SCADEN [Menden et al., Science advances, 6(30):eaba2619, 2020] is a computational method which employs deep neural network trained on single-cell RNA sequencing data to engineer discriminative features.

[0009] SUMMARY OF THE INVENTION

[0010] According to some embodiments of the invention there is provided a method of inferring cell-type abundance in a biological sample. The method comprises: receiving bulk gene expression data of the biological sample; feeding the bulk gene expression data as an input to an artificial neural network trained to apply a non-negative matrix factorization to factorize the bulk gene expression data into a matrix multiplication of a non-negative cell type signature matrix by a nonnegative cell type fraction matrix; and receiving from the network, as an output, the cell type fraction matrix, thereby inferring the cell-type abundance in the biological sample.

[0011] According to some embodiments of the invention the artificial neural network comprises a plurality of layers, each corresponding to one iteration of the non-negative matrix factorization, the iteration updates each column of the cell type fraction matrix by applying to the column an entry wise multiplication by an entry wise ratio between a first vector and a second vector.

[0012] According to some embodiments of the invention the first vector is calculated by multiplying a first learned matrix by a respective column of the bulk gene expression data, and the second vector is calculated by multiplying a second learned matrix by the column of the cell type fraction matrix.

[0013] According to some embodiments of the invention the artificial neural network has at least three layers.

[0014] According to some embodiments of the invention the method comprises receiving from the network also the cell type signature matrix.

[0015] According to some embodiments of the invention the method comprises feeding the bulk gene expression data, without feeding any input cell type signature matrix.

[0016] According to some embodiments of the invention the method comprises feeding only the bulk gene expression data.

[0017] According to some embodiments of the invention the method comprises assaying the biological sample so as to provide the bulk gene expression data.

[0018] According to an aspect of some embodiments of the present invention there is provided a system for inferring cell-type abundance in a biological sample. The system comprises: an input circuit for receiving a bulk gene expression data of the biological sample; and a hardware processor configured to execute the method as delineated above and optionally and preferably as further detailed below.

[0019] According to an aspect of some embodiments of the present invention there is provided a method of training an artificial neural network to infer cell-type abundance in a biological sample. The method comprises: obtaining a training bulk gene expression data and a training cell type fraction matrix, and applying supervised learning minimizing a first loss function to an untrained artificial neural network a first plurality of times to provide a first trained network model inferring the training cell type fraction matrix from the training bulk gene expression data. The method also comprises applying unsupervised learning minimizing a second loss function to the first trained network model a second plurality of times to provide a second trained network model inferring the training cell type fraction matrix from the training bulk gene expression data, and recording the second trained network model in a computer readable medium.

[0020] According to some embodiments of the invention the first loss function is calculated using the training cell type fraction matrix.

[0021] According to some embodiments of the invention the method comprises receiving an input cell type signature matrix, wherein the second loss function is calculated using the input cell type signature matrix.

[0022] According to some embodiments of the invention the training bulk gene expression data comprises synthetic data, and the training cell type fraction matrix is a synthetic matrix.

[0023] According to some embodiments of the invention the method comprises: receiving a nonsynthetic training bulk gene expression data and a non-synthetic training cell type fraction matrix; retraining the second network model by applying unsupervised learning thereto to infer the nonsynthetic training cell type fraction matrix from the non-synthetic training bulk gene expression data, thereby provide a retrained network model; and redefining the second trained network model as the retrained network model.

[0024] According to some embodiments of the invention the method comprises synthesizing the training bulk gene expression data and the training cell type fraction matrix.

[0025] According to some embodiments of the invention synthesis of the training cell type fraction matrix comprises generating a set of random vectors according to a multivariate probability distribution, and constructing the training cell type fraction matrix using the set as matrix columns or matrix rows.

[0026] According to some embodiments of the invention the multivariate probability distribution is selected from the group consisting of a Dirichlet distribution, a multinomial distribution, a multivariate normal distribution, multivariate logistic distribution, Wishart distribution, a multivariate hypergeometric distribution, a copula distribution, and a multivariate Poisson distribution.

[0027] According to some embodiments of the invention the method comprises receiving an input cell type signature matrix, wherein the synthesizing the training bulk gene expression data comprises multiplying the synthetic matrix by the input cell type signature matrix.

[0028] According to an aspect of some embodiments of the present invention there is provided a system for training an artificial neural network to infer cell-type abundance in a biological sample. The system comprises: an input circuit for receiving a non-synthetic training bulk gene expression data and a non-synthetic training cell type fraction matrix; and a hardware processor configured to execute the method as delineated above and optionally and preferably as further detailed below.

[0029] According to an aspect of some embodiments of the present invention there is provided a computer software product. The computer software product comprises a computer-readable medium in which program instructions are stored, which instructions, when read by a computer, cause the computer to receive a bulk gene expression data of the biological sample and to execute any of the methods delineated above and optionally and preferably any of the methods as further detailed below.

[0030] Unless otherwise defined, all technical and / or scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the invention pertains. Although methods and materials similar or equivalent to those described herein can be used in the practice or testing of embodiments of the invention, exemplary methods and / or materials are described below. In case of conflict, the patent specification, including definitions, will control. In addition, the materials, methods, and examples are illustrative only and are not intended to be necessarily limiting.

[0031] Implementation of the method and / or system of embodiments of the invention can involve performing or completing selected tasks manually, automatically, or a combination thereof. Moreover, according to actual instrumentation and equipment of embodiments of the method and / or system of the invention, several selected tasks could be implemented by hardware, by software or by firmware or by a combination thereof using an operating system.

[0032] For example, hardware for performing selected tasks according to embodiments of the invention could be implemented as a chip or a circuit. As software, selected tasks according to embodiments of the invention could be implemented as a plurality of software instructions being executed by a computer using any suitable operating system. In an exemplary embodiment of the invention, one or more tasks according to exemplary embodiments of method and / or system as described herein are performed by a data processor, such as a computing platform for executing a plurality of instructions. Optionally, the data processor includes a volatile memory for storing instructions and / or data and / or a non-volatile storage, for example, a magnetic hard-disk and / or removable media, for storing instructions and / or data. Optionally, a network connection is provided as well. A display and / or a user input device such as a keyboard or mouse are optionally provided as well. BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING(S)

[0033] Some embodiments of the invention are herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of embodiments of the invention. In this regard, the description taken with the drawings makes apparent to those skilled in the art how embodiments of the invention may be practiced.

[0034] In the drawings:

[0035] FIG. 1 is a schematic illustration providing high level pseudo-code of a procedure used in experiments performed according to some embodiments of the present invention;

[0036] FIG. 2 is a bar graph showing contributions of different training stages of a training method employed in experiments performed according to some embodiments of the present invention;

[0037] FIGs. 3A and 3B show performance evaluation on simulated data, obtained in experiments performed according to some embodiments of the present invention;

[0038] FIGs. 4 A and 4B show performance evaluation on real data, obtained in experiments performed according to some embodiments of the present invention;

[0039] FIGs. 5A-E are scatter plots of ground truth and values predicted by several cell-type abundance techniques, as obtained in experiments performed according to some embodiments of the present invention;

[0040] FIG. 6 is a flowchart diagram describing a method suitable for inferring cell-type abundance in a biological sample under analysis, according to some embodiments of the present invention;

[0041] FIGs. 7A and 7B are schematic illustration of a neural network architecture according to some embodiments of the present invention;

[0042] FIG. 8 is a flowchart diagram of a method suitable for training an artificial neural network to infer cell-type abundance in a biological sample, according to some embodiments of the present invention;

[0043] FIG. 9 is a flowchart diagram of a procedure suitable for generating synthesized training data, according to some embodiments of the present invention; and

[0044] FIG. 10 is a schematic illustration of a computing platform having a client computer and a server computer according to some embodiments of the present invention.

[0045] DESCRIPTION OF SPECIFIC EMBODIMENTS OF THE INVENTION

[0046] The present invention, in some embodiments thereof, relates to computational biology and, more particularly, but not exclusively, to a method and system for inferring cell-type abundance. Before explaining at least one embodiment of the invention in detail, it is to be understood that the invention is not necessarily limited in its application to the details of construction and the arrangement of the components and / or methods set forth in the following description and / or illustrated in the drawings and / or the Examples. The invention is capable of other embodiments or of being practiced or carried out in various ways.

[0047] Information pertaining the specific types and proportions of cells present in a given sample is useful for gaining insights into tissue function, disease mechanisms, and potential therapeutic targets. Such information can also provide insights regarding the dynamics of cellular interactions, microenvironmental influences, and disease-associated changes. The present embodiments provide a method and a system that infer cell-type abundance in a biological sample, and may therefore be used in many applications that require knowledge regarding cell types and cell proportions within a biological sample containing expressed genes, including, without limitation, biomarker identification, personalized medicine, drug development, aging research, environmental exposure study, immunology, neuroscience, regenerative medicine, and the like.

[0048] Some embodiments of the present invention are described below by means of flowchart diagrams describing methods according to various exemplary embodiments of the present invention. It is to be understood that, unless otherwise defined, the operations described hereinbelow can be executed either contemporaneously or sequentially in many combinations or orders of execution. Specifically, the ordering of the flowchart diagrams is not to be considered as limiting. For example, two or more operations, appearing in the following description or in the flowchart diagrams in a particular order, can be executed in a different order (e.g., a reverse order) or substantially contemporaneously. Additionally, several operations described below are optional and may not be executed.

[0049] Referring now to the drawings, FIG. 6 is a flowchart diagram describing a method suitable for inferring cell-type abundance in one or more biological samples under analysis, according to some embodiments of the present invention.

[0050] At least part of the operations described herein can be implemented by a data processing system, e.g., a dedicated circuitry or a general purpose processor, configured for executing the operations described below. At least part of the operations can be implemented by a cloudcomputing facility at a remote location.

[0051] Computer programs implementing the method of the present embodiments can commonly be distributed to users by a communication network or on a distribution medium such as, but not limited to, a floppy disk, a CD-ROM, a flash memory device and a portable hard drive. From the communication network or distribution medium, the computer programs can be copied to a hard disk or a similar intermediate storage medium. The computer programs can be run by loading the code instructions either from their distribution medium or their intermediate storage medium into the execution memory of the computer, configuring the computer to act in accordance with the method of this invention. During operation, the computer can store in a memory data structures or values obtained by intermediate calculations and pull these data structures or values for use in subsequent operation. All these operations are well-known to those skilled in the art of computer systems.

[0052] Processing operations described herein may be performed by means of processer circuit, such as a DSP, microcontroller, FPGA, ASIC, etc., or any other conventional and / or dedicated computing system.

[0053] The method of the present embodiments can be embodied in many forms. For example, it can be embodied in on a tangible medium such as a computer for performing the method operations. It can be embodied on a computer readable medium, comprising computer readable instructions for carrying out the method operations. In can also be embodied in electronic device having digital computer capabilities arranged to run the computer program on the tangible medium or execute the instruction on a computer readable medium.

[0054] Each of the biological sample(s) under analysis contains expressed genes and can be any sample obtained from an organism (e.g., a mammal, such as, but not limited to, a human subject) or from components (e.g., cells) of an organism. Each sample may be of any relevant biological tissue or fluid. The sample may be a sample derived from a subject. Such samples include, but are not limited to, sputum, blood, blood cells (e.g., white cells), amniotic fluid, plasma, semen, bone marrow, and tissue or fine needle biopsy samples, urine, peritoneal fluid, and pleural fluid, or cells therefrom. The biological sample may alternatively include one or more sections of tissues such as frozen sections taken for histological purposes.

[0055] The method begins at 10 and optionally and preferably continues to 11 at which bulk gene expression data of one or more biological samples are received. The bulk gene expression data of each sample is typically in digital form and comprises, for each sample, a list of genes and a corresponding list of expression levels, where each expression level indicates the level by which the respective gene is expressed in the respective biological sample. Each expression level can be provided as a count indicating the number of genes of the respective type that are expressed in the respective biological sample. In some embodiments of the present invention the gene expression data are provided as a matrix. Alternatively, the method can arrange the data as a matrix. Such a matrix represents the gene expression profile of the biological sample(s), and is interchangeably referred to as "the profile matrix," or "the profile." The (k,j) entry of the profile matrix contains the expression level of the kth gene in the jth biological sample.

[0056] In some embodiments of the present invention 11 is preceded by an operation in which the biological sample(s) is / are assayed so as to provide the bulk gene expression data. This can be done by any technique known in the art, including, without limitation, microarray technology, RNA sequencing, polymerase chain reaction, and the like.

[0057] The method continues to 12 at which the bulk gene expression data is fed as an input to an artificial neural network. In some embodiments of the present invention the artificial neural network is not fed by any input cell type signature matrix. These embodiments are advantageous because they allow inferring the cell-type abundance without relying on a specific cell type signature matrix. For example, the method of the present embodiments can infer abundance of a cell-type or a combination of cell-types that is not included in known signature matrices. In some embodiments of the present invention the network is fed only by the bulk gene expression data and is not fed with any other input.

[0058] The network comprises a plurality of layers (e.g., at least three layers or at least four layers or more), and is preferably trained to apply a non-negative matrix factorization (NNMF) in order to factorize the gene expression profile, V, into a matrix multiplication of a non-negative cell type signature matrix, S, by a non-negative cell type fraction matrix, F.

[0059] The non-negative cell type fraction matrix, F, includes a plurality of entries, each corresponding to a particular type of cell and containing its proportion within the biological sample(s). The proportion of a particular cell type is typically expressed as a fraction but can also be expressed as percentage. The (i,j) entry of F contains the relative abundance of the ith cell type in the jth biological sample. The matrix F is "non-negative" in the sense that each of its entries is either positive or zero. The non-negative cell type signature matrix, S, includes a plurality of entries, each corresponding to a particular type of gene and containing a list of values where each value is an expression level of the respective gene in one of the cell types. Thus, the (k,i) entry of S contains the expression level of the kth gene in the ith cell type. When the entries of the matrices S and F are arranged such that the number of columns in S (e.g., corresponding to cell types) is the same as the number of rows in F, the matrix multiplication SF provides a matrix whose entries correspond to individual genes and contain lists of expression levels of those genes in each sample. Thus, the matrix multiplication reconstructs the profile matrix.

[0060] It is appreciated that when the bulk gene expression data correspond to a single sample, the profile matrix, V, has a single column (or, equivalently a single row) and is therefore a vector. The ith element of this vector contains the expression level of the ith gene in the sample. In this case the cell type fraction matrix, F, also has a single column (or, equivalently a single row) and is therefore a vector, whose ith element contains the relative abundance of the ith cell type in the single sample.

[0061] A schematic architecture of an artificial neural network 20 suitable for the present embodiments is illustrated in FIG. 7A. In the illustrated embodiment, network 20 has a deep unrolled network architecture. Unrolling techniques connect between iterative methods and deep networks by viewing each iteration of an underlying iterative process as a layer of a network, such that concatenating the layers forms a deep neural network where the parameters of the iterative process are transferred into network parameters and are learned by training. FIG. 7B schematically illustrates a single layer 22 (the ith layer) of network 20, which corresponds to one iteration of the non-negative matrix factorization process.

[0062] For clarity of presentation, the illustrations in FIGs. 7A and 7B describe propagation of a particular column / of the matrix F (or vector, when the method receives bulk gene expression data of a single sample) through the layers of network 20. The skilled person would appreciate that a similar procedure is applied to the other column(s) of F (if exist). The version of column / at the entry to the deep unrolled network architecture is denoted / o. The version of / at the ith layer is denoted / . Each layer of the architecture includes an operation, g, which updates the version of / at the current layer and produces a version of / for the next layer. Thus, at the first layer / i is updated producing f , at the second layer f is updated producing f , and so on.

[0063] The operation g typically utilizes learned parameters that are used in order to update / . In some embodiments of the present invention the update performed at the ith layer includes applying to the current version / of the column at this layer an entry-wise multiplication by an entry- wise ratio between a first vector and a second vector.

[0064] Using upper indices to represent elements of vectors, an entry- wise ratio between a vector x, whose elements are x1, x2, ... and vector y, whose elements are y1, y2, ..., is defined as a vector z, whose elements are z^xVy1, z2=x2 / 2, ... etc. An entry-wise multiplication is similarly defined except that the individual elements are multiplied instead of divided.

[0065] In some embodiments of the present invention the first vector can be calculated as Ai+i v where Ai+i is a first learned matrix and is multiplied by a column v of the profile V, which corresponds to column / (or by V itself in cases in which the method receives bulk gene expression data of a single sample and V has single column). The second vector can be calculated as Bi+i , where Bi+i is a second learned matrix and is multiplied by the current version / of the column (see EQ. 2 in the Examples section that follows). Thus, in these embodiments, the learned parameters of the ith layer are the elements of matrices Ai+i and £>i+i . Referring again to FIG. 6, the method proceeds to 13 at which the cell type fraction matrix F is received as an output from the neural network. Since the entries of F contain proportions of cell types within the biological sample(s), the received matrix is the inference of the cell-type abundance in the sample(s). In some embodiments of the present invention the method proceeds to 14 at which the cell type signature matrix, S, is also received from the network. This can be achieved by providing the network with an additional operation (not shown in FIGs. 7 A and 7B) that calculates S based on the version of F at the last layer of the deep unrolled network architecture and based on the input profile V. Specifically, the additional operation finds an approximate solution to the matrix equation V=SF for the matrix S. For example, the additional operation can find S such that the norm of V-SF is minimized. In some embodiments of the present invention this is done by applying non-negative least squares to the profile V and the inferred matrix F.

[0066] In some embodiments of the present invention the method proceeds to 15 at which the information provided by the matrix F and optionally and preferably the matrix S is used for conducting a biological experiment. The type of biological experiment depends on the application for which the method is employed. For example, when the method is employed for an optimization of a therapy or a medicament, the experiment can include applying the therapy or administrating the medicament and monitoring changes in cell-type abundance before and after the therapy or a medicament by re-executing method 10. When the method is employed for population health studies, e.g., an environmental exposure study, the experiment can include exposing a study group or an individual to an environmental condition, and studying differences in cell-type abundance before and after the exposure or between a study group that has been exposed and a control group that has not been exposed, by re-executing method 10. When the method is employed for regenerative medicine, the experiment can be a tissue engineering experiment directed to engineer a tissue having the same or different cell-type abundance as the cell-type abundance inferred by method 10. Also in the field of regenerative medicine, the experiment can include applying a regenerative therapy and monitoring changes in cell-type abundance before and after the therapy by re-executing method 10.

[0067] The biological experiment can alternatively or additionally include re-executing the method at a later time, comparing cell-type abundances obtained at different executions of the method, and determining an influence of a biological or environmental condition or event on the cell-type abundance of one or more subjects.

[0068] The method ends at 16.

[0069] Reference is now made to FIG. 8, providing a flowchart diagram of a method suitable for training an artificial neural network to infer cell-type abundance in a biological sample, according to various exemplary embodiments of the present invention. The biological sample can be of any of the sample types described above. The method can be used for training an artificial neural network, which, once trained, can be used for executing method 10, as described above with reference to FIG. 6. The architecture of the artificial neural network to be trained is selected to apply NNMF in a manner that for a profile, V, the output of the network, once trained, comprises a non-negative cell type fraction matrix, F, that relates to the input profile via the relation V=SF, where S is a signature matrix that may optionally and preferably also be obtained from the network, as further detailed hereinabove.

[0070] The method begins at 30 and proceeds to 31 at which training data are obtained. The training data preferably include at least training bulk gene expression data and a training cell type fraction matrix. These data are of the same type and form as the bulk gene expression data and a cell type fraction matrix described above, except that in method 30 both types of data are known and are used for training. The training bulk gene expression data are preferably bulk gene expression data of a plurality of biological samples. When represented as a matrix, the (k,j) entry of the training bulk gene expression data contains the expression level of the kth gene in the jth biological sample.

[0071] The training data can be data obtained in previously executed biological assays, in which case the data is referred to as "real." Alternatively, at least a portion of the training data can be synthesized data. The training data is optionally and preferably provided in digital form and can be read from a computer readable medium or retrieved over a communication network, such as, but not limited to, the internet. A representative example of a procedure suitable for generating synthesized training data is described below with reference to FIG. 9.

[0072] In some embodiments of the present invention the training data also comprises an input cell type signature matrix, which can be of the same type and form as the cell type signature matrix described above except that it is known rather than inferred. The input cell type signature matrix is optionally and preferably provided in digital form and can be read from a computer readable medium or retrieved over a communication network, such as, but not limited to, the internet. Representative examples of cell type signature matrices that can be used as input according to some embodiments of the present invention include, without limitation, the cell type signature matrices disclosed in Reference Nos. [7], [8],

[0013] , and

[0018] , the contents of which are hereby incorporated by reference.

[0073] It is expected that during the life of a patent maturing from this application many relevant cell type signature matrices will be developed and updated, and the scope of the term input cell type signature matrix is intended to a priori include any matrix or matrices selected from these developed and / or updated matrices.

[0074] Also contemplated, are embodiments in which instead of, or in addition to, receiving an input cell type signature matrix as part of the training data, the training data includes single cell data. In these embodiments, the single cell data can be used to learn an atlas of cell type specific profiles, which can then be used to initialize a cell type signature matrix.

[0075] In some embodiments of the present invention the training comprises both synthetic and non- synthetic (namely real) training data. In these embodiments, the method can optionally and preferably be trained on the synthetic training data and separately on the non- synthetic training data. The non- synthetic training data can include non- synthetic training bulk gene expression data and non-synthetic training cell type fraction matrix.

[0076] The method proceeds to 32 at which a supervised learning minimizing a first loss function is applied to an untrained artificial neural network, having a plurality of learnable parameters. The untrained artificial neural network preferably have a deep unrolled network architecture, such as, but not limited to, the architecture illustrated in FIGs. 7A and 7B, above, in which case the learnable parameters are stored in the matrices Ai and Bi. Since the artificial neural network is untrained, the learnable parameters of the network are yet to be determined, and are therefore assigned with some initial values, which may be arbitrary. For example, the matrices Ai and Bi can all be identity matrices. The supervised learning is optionally and preferably applied a plurality of times in an iterative manner, so that starting from the second application of the supervised learning, the initial parameters of the supervised learning (e.g., the matrix elements of Ai and Bi) are those that were obtained during the previous supervised learning, and are therefore no longer arbitrary.

[0077] The loss function that is minimized by the supervised learning is preferably calculated based on the training cell type fraction matrix. For example, the loss function can be a distance (e.g., an L2 distance) of the training cell type fraction matrix to a cell type fraction matrix that is factorized by the network out of the training bulk gene expression data. Operation 32 thus provides a first trained network model. The first network model is trained in the sense that its parameters have been set during operation 32 to allow the model to infer the training cell type fraction matrix from the training bulk gene expression data.

[0078] The method proceeds to 33 at which an unsupervised learning minimizing a second loss function is applied to the first trained network model. As in the case of 32, the unsupervised learning is optionally and preferably applied a plurality of times in an iterative manner. The initial parameters of the unsupervised learning are those that were obtained during the last application of the supervised learning. Thereafter, the initial parameters of the unsupervised learning (e.g., the matrix elements of Ai and Bi) are those that were obtained during the previous unsupervised learning.

[0079] The difference between the second loss function (minimized by the unsupervised learning applied at 33) and the first loss function (minimized by the supervised learning applied at 32) is that the second loss function is calculated irrespectively of the training cell type fraction matrix. In embodiments in which the training data comprises an input cell type signature matrix, the second loss function can be calculated based on this input cell type signature matrix. For example, the loss function can be a distance (e.g., an L2 distance) of the training bulk gene expression data to the matrix multiplication between the input cell type signature matrix and the cell type fraction matrix that is factorized out by the network of the training bulk gene expression data.

[0080] In embodiments in which the method is trained on the synthetic training data and separately on the non-synthetic training data, operation 32 and 33 are executed using the synthetic training data, and the method continues to 34 at which the second trained network model (obtained at 33) is retrained using the non-synthetic training data. The retraining is preferably by applying unsupervised learning as further detailed hereinabove. The loss function that is minimized during the unsupervised learning at 34 can be the same as the loss function that is minimized during the unsupervised learning at 33. As in the case of 32 and 33, the unsupervised learning at 34 is optionally and preferably applied a plurality of times in an iterative manner, so that starting from the second application of the unsupervised learning, the initial parameters of the supervised learning (e.g., the matrix elements of Ai and Bi) are those that were obtained during the previous unsupervised learning.

[0081] The number of times at which the unsupervised learning at 34 is executed can be smaller than the number of times at which the supervised learning at 32 and the number of times at which the unsupervised learning at 33 are employed. For example, the supervised learning at 32 can be executed about 10,000-100,000 times, e.g., about 60,000 times, the unsupervised learning at 33 can be executed about 10,000-100,000 times, e.g., about 60,000 times, and the unsupervised learning at 34 can be executed about 50-200 times, e.g., about 100 times. Other numbers of executions are also contemplated.

[0082] The method can then proceed to 35 at which at least the second network model is recorded on a computer readable medium. In embodiments in which 34 is not executed, the second trained network model obtained at 33 is recorded. In embodiments in which 34 is executed, the second retrained network model obtained at 34 is recorded. Also contemplated are embodiments in which more than one trained network models are recorded (e.g., the models obtained at 32 and 34, or the models obtained at 33 and 34, or the models obtained at 32, 33 and 34). Following operation 35, the method can re-execute selected operations of method 10 described above with reference to FIG. 6.

[0083] The method ends at 36.

[0084] Reference is now made to FIG. 9, which is a flowchart diagram of a procedure suitable for generating synthesized training data, according to some embodiments of the present invention. The procedure can be executed for generating at least a portion of the training data obtained at operation 31 of method 30.

[0085] The procedure begins at 40 and continues to 41 at which a set of random vectors is generated according to a probability distribution, such as, but not limited to, a multivariate probability distribution. Representative examples of multivariate probability distributions suitable for the present embodiments include, without limitation, a Dirichlet distribution, a multinomial distribution, a multivariate normal distribution, multivariate logistic distribution, Wishart distribution, a multivariate hypergeometric distribution, a copula distribution, and a multivariate Poisson distribution. In experiments performed by the Inventors, a Dirichlet distribution has been employed. The procedure continues to 42 at which a training cell type fraction matrix is constructed using the generated set of random vectors. The vectors can be used as the columns or rows of the training cell type fraction matrix.

[0086] The procedure optionally and preferably continues to 43 at which training bulk gene expression data are generated from the training cell type fraction matrix constructed at 42. For example, the matrix constructed at 42 can be multiplied by the input cell type signature matrix, thereby providing the bulk gene expression data in the form of a matrix.

[0087] The procedure ends at 44.

[0088] FIG. 10 is a schematic illustration of a client computer 130 having a hardware processor 132, which typically comprises an input / output (I / O) circuit 134, a hardware central processing unit (CPU) 136 (e.g., a hardware microprocessor), and a hardware memory 138 which typically includes both volatile memory and non-volatile memory. CPU 136 is in communication with I / O circuit 134 and memory 138. Client computer 130 preferably comprises a graphical user interface (GUI) 142 in communication with processor 132. I / O circuit 134 preferably communicates information in appropriately structured form to and from GUI 142. Client computer 130 can be part of a computing platform having also a server computer 150 which can include a hardware processor 152, an I / O circuit 154, a hardware CPU 156, a hardware memory 158.

[0089] I / O circuits 134 and 154 of client 130 and server 150 computers can operate as transceivers that communicate information with each other via a wired or wireless communication. For example, client 130 and server 150 computers can communicate via a network 140, such as a local area network (LAN), a wide area network (WAN) or the Internet. Server computer 150 can be in some embodiments be a part of a cloud computing resource of a cloud computing facility in communication with client computer 130 over the network 140.

[0090] GUI 142 and processor 132 can be integrated together within the same housing or they can be separate units communicating with each other. GUI 142 can optionally and preferably be part of a system including a dedicated CPU and I / O circuits (not shown) to allow GUI 142 to communicate with processor 132. Processor 132 issues to GUI 142 graphical and textual output generated by CPU 136. Processor 132 also receives from GUI 142 signals pertaining to control commands generated by GUI 142 in response to user input. GUI 142 can be of any type known in the art, such as, but not limited to, a keyboard and a display, a touch screen, and the like.

[0091] Client 130 and server 150 computers can further comprise one or more computer-readable storage media 144, 164, respectively. Media 144 and 164 are preferably non-transitory storage media storing computer code instructions for executing one or more of the methods and procedures as further detailed herein, and processors 132 and 152 execute these code instructions. The code instructions can be run by loading the respective code instructions into the respective execution memories 138 and 158 of the respective processors 132 and 152.

[0092] One or both of storage media 144 and 164 can store a trained artificial neural network, and program instructions which, when read by the respective processor, cause the processor to receive a bulk gene expression data of the biological sample as further detailed hereinabove. The program instructions can also cause the processor to feed the bulk gene expression data as an input to the artificial neural network, as further detailed hereinabove.

[0093] In some embodiments of the present invention the bulk gene expression data are received by computer 130 locally, for example, using storage 144 and / or by means of a measurement and analyses system 146 configured to assay a biological sample so as to provide the bulk gene expression data, as further detailed hereinabove. In these embodiments computer 130 can execute the operations of the method described herein. Alternatively, computer 130 can transmit the received bulk gene expression data to computer 150 via communication network 140, in which case computer 130 executes the operations of the method described herein, and transmits the output of the network (e.g., the cell type fraction matrix and optionally and preferably also the cell type fraction matrix) back to computer 130 for storage in medium 144 and / or generating a displayed output on GUI 142.

[0094] One or both of storage media 144 and 164 can alternatively store an untrained artificial neural network, and program instructions which, when read by the respective processor, cause the processor to receive training data as further detailed hereinabove and execute the method described herein for training the artificial neural network. In some embodiments of the present invention the training data are received by computer 130 locally, for example, using storage 144 and / or by means of a measurement and analyses system 146. In these embodiments computer 130 can execute the operations of the training method described herein. Alternatively, computer 130 can transmit the training data to computer 150 via communication network 140, in which case computer 130 executes the operations of the training method described herein, and transmits the trained network or the parameters thereof back to computer 130 for storage in medium 144.

[0095] As used herein the term “about” refers to ± 10 %

[0096] The terms "comprises", "comprising", "includes", "including", “having” and their conjugates mean "including but not limited to".

[0097] The term “consisting of’ means “including and limited to”.

[0098] The term "consisting essentially of" means that the composition, method or structure may include additional ingredients, steps and / or parts, but only if the additional ingredients, steps and / or parts do not materially alter the basic and novel characteristics of the claimed composition, method or structure.

[0099] As used herein, the singular form "a", "an" and "the" include plural references unless the context clearly dictates otherwise. For example, the term "a compound" or "at least one compound" may include a plurality of compounds, including mixtures thereof.

[0100] Throughout this application, various embodiments of this invention may be presented in a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the invention. Accordingly, the description of a range should be considered to have specifically disclosed all the possible subranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range, for example, 1, 2, 3, 4, 5, and 6. This applies regardless of the breadth of the range.

[0101] Whenever a numerical range is indicated herein, it is meant to include any cited numeral (fractional or integral) within the indicated range. The phrases “ranging / ranges between” a first indicate number and a second indicate number and “ranging / ranges from” a first indicate number “to” a second indicate number are used herein interchangeably and are meant to include the first and second indicated numbers and all the fractional and integral numerals therebetween. It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, may also be provided separately or in any suitable subcombination or as suitable in any other described embodiment of the invention. Certain features described in the context of various embodiments are not to be considered essential features of those embodiments, unless the embodiment is inoperative without those elements.

[0102] Various embodiments and aspects of the present invention as delineated hereinabove and as claimed in the claims section below find experimental support in the following examples.

[0103] EXAMPLES

[0104] Reference is now made to the following examples, which together with the above descriptions illustrate some embodiments of the invention in a non limiting fashion.

[0105] This Example present DEep Cell-type DEconvolution (DECODE), a deep learning method for the task that builds on a deep unfolded non-negative matrix factorization technique. This Example demonstrates that the method of the present embodiments outperforms previous approaches on a range of synthetic and real data sets.

[0106] Introduction

[0107] Biological tissues are composed of a variety of distinct cell types. Identifying the composition of cells in tissues can help generate hypotheses regarding cell-type-specific biological mechanisms with important biomedical applications. For example, patients with a large number of infiltrating T cells are more likely to respond positively to immunotherapy

[0016] .

[0108] Flow cytometry is the main standard for experimental deconvolution of a sample. More recently, single-cell RNA sequencing (scRNA-seq) methods have become available. However, flow cytometry requires prompt and careful processing of samples as well as tissue disaggregation, which may result in the loss of fragile cell types and the distortion of gene expression profiles, and ScRNA-seq methods are expensive for large sample studies. Additionally, the Inventors found that cell types such as neurons, myocytes, and adipocytes are difficult to capture by these technologies due to cell size and morphology.

[0109] Several computational methods are known for predicting cell fractions from bulk expression data. Most methods rely on a signature matrix of cell-specific expression profiles to predict the cell type abundance. Recent comparative analyses of deconvolution methods [2—4, 9, 11, 14, 15, 17, 19] have highlighted state-of-the-art methods for this task including non-negative least squares (NNLS) [5], CIBERSORT

[0013] , CIBERSORTx

[0014] , GEDIT

[0010] and SCADEN [9]. The Inventirs found that most of these methods rely heavily on the input signature matrices, which are global matrices that do not contain information that is specific to the input tissue.

[0110] DECODE employs a deep-learning procedure to predict the cell type abundance matrix from bulk gene expression data and signature matrix. The method is based on a deep unfolding procedure for non-negative matrix factorization (NMF) and combines both supervised learning on synthetic data and unsupervised learning to achieve its task. The Inventors benchmarked DECODE using both simulated and real data sets and showed that it outperforms previous approaches.

[0111] In DECODE signature matrices are only used to initialize the model and to generate training data, but are not explicitly represented by the model itself. DECODE is based on a flexible neural network architecture allowing it to use non-negative matrix factorization techniques for simultaneous prediction of cell fractions and signatures. This is advantageous over known methods that depend on signature matrices, since they do not guarantee that cell fraction vectors will sum to one. DECODE combines supervised and unsupervised training. The Inventors found that this is advantageous as the training is unsupervised but the evaluation is according to true (hidden) cell fractions.

[0112] DECODE generates synthetic data and is trained on that synthetic data. This overcomes the small amount of available training data.

[0113] Methods

[0114] In gene expression deconvolution, the input is a matrix of bulk gene expression across multiple samples and a signature matrix consisting of expression profiles of specific cell types. The goal is to infer a matrix of cell fractions indicating for each sample its cell-type decomposition.

[0115] In this Example this problem is addressed using a deep learning procedure for non-negative matrix factorization that aims to factor the input gene expression matrix into the product of the signature matrix and the cell fraction matrix. Due to scarcity of training data, the parameters deep learning procedure are trained using a combination of synthetically generated data and real data. A high level pseudo-code of the procedure is illustrated in FIG. 1. A bulk expression matrix and a signature of expression profiles are received. Supervised training is applied on synthetically generated data using realistic cell fractions for I iterations. The resulting model is trained in an unsupervised fashion for another I iterations on the synthetic data. The model is then further trained on real data in an unsupervised fashion to produce its predictions.

[0116] Consider a bulk expression matrix V of n genes by m samples whose n rows represent the average of cell-type specific gene expression profiles in a sample, weighted by their abundances in that sample. Suppose there are k cell types and let S be a {\em signature matrix} whose k columns are the gene expression profiles of those cell types. The procedure optionally and preferably infers a matrix F whose m columns are probability vectors denoting the fraction of each cell type in the corresponding sample such that V~ S F. The procedure is based on DNMF, a deep unfolding approach for non-negative matrix factorization (NMF) described in

[0012] , the contents of which are hereby incorporated by reference. DNMF is a deep learning procedure for the decomposition of a non-negative matrix Vnxm into a product of two non-negative matrices Snxk and Fk xmsuch that || V - SF||2 is minimal. In more detail, DNMF contains several ( / ) layers that are updated based on the multiplicative update rule of Lee and Seung [6]. The latter rule iteratively updates S and F as follows:

[0117] (EQ. 1) where and [.] / [.] represent entry-wise multiplication and division, respectively. However, unlike in Lee and Seung, in DNMF this rule is relaxed to allow learning better solutions. Specifically, the procedure contains a layer for each iteration and aims to optimize F without explicitly representing S. Focusing on some column / of F and the corresponding column v of V, the output / of layer z is (up to regularization): where Ai+i and Si+i are learnable matrices that correspond in the original update rule to the matrices S1and STS, respectively (ignoring the dependency between the two latter matrices). Starting from an initial matrix Fo, the procedure rolls it through the / layers of the network, imitating the Lee and Seung iterative procedure, until an output Fi matrix is produced.

[0118] The DNMF procedure has two model variants: supervised and unsupervised. The supervised variant assumes a known fraction matrix Fre iwhich is either experimentally measured or generated in simulations. It is trained with an L loss w.r.t. this matrix: 11 Fi - Freai| I2. When no real matrix is available, DNMF uses an unsupervised variant where the loss is the reconstruction error || V- SFi\ |2. This Example uses these two variants into a combined supervised-unsupervised DNMF model that maps from bulk expression (V) to cell fractions (F).

[0119] Following is a description of the DNMF architecture and training process used in the present Example according to some embodiments of the present invention.

[0120] To account for the fact that each column / of F represents a probability vector, / was normalized after each iteration to sum to one. To make use of the given signature matrix S, a DNMF model M(S) was initialized. In the initialized model, the weights of the first layer were selected according to one iteration of Lee and Seung's update rule, setting AI=STand BI=STS. In addition, Fo was initialized with the result of applying NNLS to V and S. All other parameters are initialized according to

[0012] . The training process included a combined supervised and unsupervised training, where the former was based on synthetic data and the the latter was based on both synthetic and real data.

[0121] The model was first trained using synthetic data in a supervised fashion for one epoch with I iterations (or batches). In the present Example, I was set to be 60,000. Let Mk be the resulting model after the A- th iteration (for clarity, the signature matrix S is omitted from the notation). With this notation, the result of the first I iterations is Mi. Thereafter, Mi was trained for another I iterations in an unsupervised fashion on synthetic data to yield an intermediate model Mu. The obtained model Muwas then trained for additional Irepochs on real data to provide a signaturespecific model. In the present Example, Irwas set to be 100. The process was executed for each of a given list of signature matrices and the model with lowest unsupervised error on the real data was selected for predicting the cell-type abundance.

[0122] Supervised training on synthetic data

[0123] To deal with data scarcity, the procedure was first trained on synthetic data. Specifically, the Inventors collected known ranges of cell frequencies from [1] and used values within these ranges as parameters for a Dirichlet distribution. Random fraction vectors were then generated from this distribution, and a fraction matrix was constructed from the generated vectors. In this Example, the generated vectors were defined as the rows of the fraction matrix. The fraction matrix was multiplied by a signature matrix to produce a synthetic bulk expression matrix.

[0124] In this Example, four signature matrices were tested. These included LM22

[0013] , Human Primary Cell Atlas (HPCA-Blood) [7], Blue-Code [8], and Skin Signatures

[0018] . For each signature matrix, a separate fraction matrix was constructed by regenerating random fraction vectors from the same Dirichlet distribution. The procedure therefore resulted in four different synthetic bulk expression matrices to be tested, one synthetic bulk expression matrix for each signature matrix.

[0125] The synthetic bulk expression matrices were fed sequentially to the training process, viewing each such matrix as representing a batch. A small normally-distributed noise was added to each bulk expression matrix with zero mean and small standard deviation. In this Example, the standard deviation of the for the zth batch was z710,000.

[0126] Unsupervised training on synthetic data

[0127] In the unsupervised training, the model Mi was trained to minimize the reconstruction error 11 V- SFmodei||2. Before training the model, the NNLS supervised loss ennis was computed on the synthetic data. The stopping criterion for the training iteration was either when the number of iterations reached I or when the supervised loss of the jth iteration exceeded the value of ennis, where j<I-

[0128] Hyperparameter tuning

[0129] In this Example, the number of layers was selected to be 4, and the learning rate was selected to be 0.001. No regularization was used in this Example. The numbers of training iterations (7 and Ir) were selected by performing a grid search using an independent stromal dataset from

[0011] . The search was conducted over the ranges 40,000 < I < 80,000 and 50 < Ir< 200. The search revealed the combination of 1=60,000 and Ir=100. FIG. 2 is a bar graph showing contribution of different training stages to the final result, when applied to the dataset in

[0011] , where the bars a, b, and c correspond to no supervised training on synthetic data, no unsupervised training on synthetic data, and unsupervised training on real data only, and the bar DECODE correspond to the full training procedure. As shown, each of the training stages aids in improving the performance.

[0130] Data and preprocessing

[0131] The main data source is

[0011] which has three available datasets, all the results of in-silico simulations. Two of the datasets, PBMC1 and PBMC2, represent 200 simulated mixtures of singlecell peripheral blood mononuclear cell (PBMC) expression profiles (100 mixtures in each dataset). For each mixture, individual cells were randomly chosen and their expression profiles summed. Both contain five common PBMC cell types (B, CD4 T, CD8 K, NK and monocytes). A third independent dataset, STROMAL, contains 100 simulated mixtures of stromal cell types (B, CD4 T, CD8 T, macrophage, mast, endothelial and fibroblast cells). In this Example, the first two were used the for testing and the third was used for hyperparameter tuning.

[0132] In addition, a real, GSE65133, dataset was retrieved from

[0013] . This dataset contains 20 samples of real PBMC cell fractions that were measured by flow cytometry.

[0133] To complement the expression datasets, known signature matrices were also used. These were taken from

[0011] , which contains a comparative analysis of 9 deconvolution methods with respect to 10 signature matrices. Among the top performing methods were CIBERSORT, NNLS and GEDIT which are compared to hereinbelow.

[0134] The accuracy values reported in

[0010] (see Figure 1, therein) were averaged for each of the signature matrices. In this Example, the following signatures, showing the highest accuracy value, were selected: (i) LM22

[0013] containing 22 cell types and 547 genes; (ii) Human Primary Cell Atlas (HPCA-Blood) [7] containing 7 cells and 19,715 genes; (iii) Blue-Code [8] containing 34 cells and 13,299 genes; and (iv) Skin Signatures

[0018] containing 21 cells and 20,307 genes. The data were preprocessed to remove cell types that are not present in either the input expression or signature matrix. In this Example the GEDIT's approach

[0010] was used for this removal. Quantile normalization was applied to both the bulk expression matrix and the signature matrix such that each column follows the same distribution as all other columns. Genes that are missing from either matrix were removed, and a subset of 50 genes with lowest entropy was selected for each cell to focus on. The advantage of this operation is that it ensures that for each cell the genes are expressed in a cell type-specific manner. Entropy was minimized when expression is detected only in a single cell type.

[0135] Comparison of performances

[0136] Several performance measures were used to compare DECODE to four existing cell deconvolution algorithms: CIBERSORTX

[0014] , NNLS [5], GEDfT

[0010] and SCADEN[9]. GEDfT was executed with its R source code. CIBERSORTX was executed from its official website www(dot)cibersortxDOTstanford(dot)edu. NNLS was executed with its R function. Scaden was executed with its Python source code and its default training datasets were kept, since it does not train with a signature matrix.

[0137] To compare the performance of the five deconvolution methods, both root mean squared error (RMSE) and Pearson correlation coefficient were measured, and the real and predicted cell fractions estimates were compared.

[0138] Results

[0139] The technique of the present embodiments was tested on two simulated datasets of PBMC cells from

[0011] . FIGs. 3A and 3B show performance evaluation on the simulated data, showing RMSE performance (FIG. 3A), and Pearson correlation performance (FIG. 3B). As shown the technique of the present embodiments (denoted DECODE in FIGs. 3A and 3B) demonstrates superiority over the other methods with respect to both the (root) mean squared error and the Pearson correlation.

[0140] The technique of the present embodiments was also applied to a real dataset of PBMCs from

[0013] . FIGs. 4A and 4B show performance evaluation on the real data, showing RMSE performance (FIG. 4A), and Pearson correlation performance (FIG. 4B). FIGs. 5A-E are scatter plots of ground truth (horizontal axis) and predicted values (vertical axis) for the technique of the present embodiments (FIG. 5A), and the conventional techniques SCADEN (FIG. 5B), CIBERSORTX (FIG. 5C), GEDIT (FIG. 5D) and NNLS (FIG. 5D), as tested on real data. The straight lines in FIGs. 5A-E are 45° prediction lines. As shown in FIGs. 4A-5E, the technique of the present embodiments (denoted DECODE in FIGs. 4A, 4B, and 5A) demonstrates superiority over the other methods with respect to the (root) mean squared error, the Pearson correlation, and the deviation from the 45° prediction line.

[0141] Table 1, below, summarizes the improvement of the technique of the present embodiments (second and last columns) with respect to the RMSE measure relative to conventional methods (third column), for three datasets.

[0142] Table 1

[0143] It is appreciated that while the methodology described in this Example does not depend on a signature matrix, such a matrix was used for the purpose of initializing the neural network. The present embodiments contemplate the inference of the signature matrix as part of the learning process so as not to depend on receiving a signature matrix it as input.

[0144] Further, due to the scarcity of real data, the technique described in this Example used synthetically generated data for the training. It is appreciated that the training can be alternatively be executed using single cell expression data, as such data will be accumulated during the lifetime of a patent maturing from this application. For example, these data can be used to simulate deconvolution scenarios and thus improve the training process.

[0145] Although the invention has been described in conjunction with specific embodiments thereof, it is evident that many alternatives, modifications and variations will be apparent to those skilled in the art. Accordingly, it is intended to embrace all such alternatives, modifications and variations that fall within the spirit and broad scope of the appended claims.

[0146] It is the intent of the applicant(s) that all publications, patents and patent applications referred to in this specification are to be incorporated in their entirety by reference into the specification, as if each individual publication, patent or patent application was specifically and individually noted when referenced that it is to be incorporated herein by reference. In addition, citation or identification of any reference in this application shall not be construed as an admission that such reference is available as prior art to the present invention. To the extent that section headings are used, they should not be construed as necessarily limiting. In addition, any priority document(s) of this application is / are hereby incorporated herein by reference in its / their entirety. REFERENCES

[0147] [1] Blood fractions. www(dot)miltenyibiotec(dot)com / US-en / resources / macs- handbook / human-cells-and-organs / human-cell-sources / blood-human.html.

[0148] [2] Francisco Avila Cobos, Jose Alquicira-Hemandez, Joseph E Powell, Pieter Mestdagh, and Katleen De Preter. Benchmarking of cell type deconvolution pipelines for transcriptomics data. Nature communications, 11(1): 1-14, 2020.

[0149] [3] Manqi Cai, Molin Yue, Tianmeng Chen, Jinling Liu, Erick Fomo, Xinghua Lu, Timothy Billiar, Juan Celedon, Chris McKennan, Wei Chen, et al. Robust and accurate estimation of cellular fraction from tissue omics data via ensemble deconvolution. Bioinformatics, 38(11):3004— 3010, 2022.

[0150] [4] Haijing Jin and Zhandong Liu. A benchmark for ma-seq deconvolution analysis under dynamic testing environments. Genome biology, 22(1): 1—23, 2021.

[0151] [5] Charles L Lawson and Richard J Hanson. Solving least squares problems. SIAM, 1995.

[0152] [6] Daniel Lee and H Sebastian Seung. Algorithms for non-negative matrix factorization. Advances in neural information processing systems, 13, 2000.

[0153] [7] Neil A Mabbott, J Kenneth Baillie, Helen Brown, Tom C Freeman, and David A Hume. An expression atlas of human primary cells: inference of gene function from coexpression networks. BMC genomics, 14(1): 1-13, 2013.

[0154] [8] Joost HA Martens and Hendrik G Stunnenberg. Blueprint: mapping human blood cell epigenomes. Haematologica, 98(10): 1487, 2013.

[0155] [9] Kevin Menden, Mohamed Marouf, Sergio Oller, Anupriya Dalmia, Daniel Sumner Magruder, Karin Kloiber, Peter Heutink, and Stefan Bonn. Deep learning-based cell composition analysis from tissue expression profiles. Science advances, 6(30):eaba2619,

[0156] 2020.

[0157]

[0010] Brian B Nadel, David Lopez, Dennis J Montoya, Feiyang Ma, Hannah Waddel, Misha M Khan, Serghei Mangul, and Matteo Pellegrini. The Gene Expression Deconvolution Interactive Tool (GEDIT): accurate cell type quantification from gene expression data. GigaScience, 10(2), 02 2021. giab002.

[0158]

[0011] Brian B Nadel, Meritxell Oliva, Benjamin L Shou, Keith Mitchell, Feiyang Ma, Dennis J Montoya, Alice Mouton, Sarah Kim-Hellmuth, Barbara E Stranger, Matteo Pellegrini, and Serghei Mangul. Systematic evaluation of transcriptomics-based deconvolution methods and references using thousands of clinical samples. Briefings in Bioinformatics, 22(6), 08

[0159] 2021. bbab265.

[0012] Rami Nasser, Yonina C Eldar, and Roded Sharan. Deep unfolding for non-negative matrix factorization with application to mutational signature analysis. Journal of Computational Biology, 29(l):45-55, 2022.

[0160]

[0013] Aaron M Newman, Chih Long Liu, Michael R Green, Andrew J Gentles, Weiguo Feng, Yue Xu, Chuong D Hoang, Maximilian Diehn, and Ash A Alizadeh. Robust enumeration of cell subsets from tissue expression profiles. Nature methods, 12(5):453-457, 2015.

[0161]

[0014] Aaron M Newman, Chloe B Steen, Chih Long Liu, Andrew J Gentles, Aadel A Chaudhuri, Florian Scherer, Michael S Khodadoust, Mohammad S Esfahani, Bogdan A Luca, David Steiner, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nature biotechnology, 37(7):773-782, 2019.

[0162]

[0015] Ellis Patrick, Mariko Taga, Ayla Ergun, Bernard Ng, William Casazza, Maria Cimpean, Christina Yung, Julie A Schneider, David A Bennett, Chris Gaiteri, et al. Deconvolving the contributions of cell-type heterogeneity on cortical gene expression. PLOS Computational Biology, 16(8):el008120, 2020.

[0163]

[0016] Yasin §enbabaoglu, Ron S Gejman, Andrew G Winer, Ming Liu, Eliezer M Van Allen, Guillermo de Velasco, Diana Miao, Irina Ostrovnaya, Esther Drill, Augustin Luna, et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger rna signatures. Genome biology, 17(1): 1-25, 2016.

[0164]

[0017] Gavin J Sutton, Daniel Poppe, Rebecca K Simmons, Kieran Walsh, Urwah Nawaz, Ryan Lister, Johann A Gagnon-Bartsch, and Irina Voineagu. Comprehensive evaluation of deconvolution methods for human brain gene expression. Nature Communications, 13(1):1— 18, 2022.

[0165]

[0018] William R Swindell, Andrew Johnston, John J Voorhees, James T Elder, and Johann E Gudjonsson. Dissecting the psoriasis transcriptome: inflammatory- and cytokine-driven gene expression in lesions from 163 patients. BMC genomics, 14(1): 1-20, 2013.

[0166]

[0019] Xuran Wang, Jihwan Park, Katalin Susztak, Nancy R Zhang, and Mingyao Li. Bulk tissue cell type deconvolution with multi-subject single-cell expression reference. Nature communications, 10(l):l— 9, 2019.

Claims

WHAT IS CLAIMED IS:

1. A method of inferring cell-type abundance in a biological sample, the method comprising: receiving bulk gene expression data of the biological sample; feeding said bulk gene expression data as an input to an artificial neural network trained to apply a non-negative matrix factorization to factorize said bulk gene expression data into a matrix multiplication of a non-negative cell type signature matrix by a non-negative cell type fraction matrix; and receiving from said network, as an output, said cell type fraction matrix, thereby inferring the cell-type abundance in the biological sample.

2. The method according to claim 1, wherein said artificial neural network comprises a plurality of layers, each corresponding to one iteration of said non-negative matrix factorization, said iteration updates each column of said cell type fraction matrix by applying to said column an entry wise multiplication by an entry wise ratio between a first vector and a second vector.

3. The method according to claim 2, wherein said first vector is calculated by multiplying a first learned matrix by a respective column of said bulk gene expression data, and said second vector is calculated by multiplying a second learned matrix by said column of said cell type fraction matrix.

4. The method according to any of claims 1-3, wherein said artificial neural network has at least three layers.

5. The method according to any of claims 1-4, comprising receiving from said network also said cell type signature matrix.

6. The method according to any of claims 1-5, wherein said feeding comprises feeding said bulk gene expression data, without feeding any input cell type signature matrix.

7. The method according to any of claims 1-6, wherein said feeding comprises feeding only said bulk gene expression data.

8. The method according to any of claims 1-7, comprising assaying the biological sample so as to provide said bulk gene expression data.

9. A computer software product, comprising a computer-readable medium in which program instructions are stored, which instructions, when read by a computer, cause the computer to receive a bulk gene expression data of the biological sample and to execute the method according to any of claims 1-7.

10. A system for inferring cell-type abundance in a biological sample, the system comprising: an input circuit for receiving a bulk gene expression data of the biological sample; and a hardware processor configured to execute the method according to any of claims 1-7.

11. A method of training an artificial neural network to infer cell-type abundance in a biological sample, the method comprising: obtaining a training bulk gene expression data and a training cell type fraction matrix; applying supervised learning minimizing a first loss function to an untrained artificial neural network a first plurality of times to provide a first trained network model inferring said training cell type fraction matrix from said training bulk gene expression data; applying unsupervised learning minimizing a second loss function to said first trained network model a second plurality of times to provide a second trained network model inferring said training cell type fraction matrix from said training bulk gene expression data; and recording said second trained network model in a computer readable medium.

12. The method according to claim 11, wherein said first loss function is calculated using said training cell type fraction matrix.

13. The method according to any of claims 11 and 12, comprising receiving an input cell type signature matrix, wherein said second loss function is calculated using said input cell type signature matrix.

14. The method according to claim 13, wherein said training bulk gene expression data comprises synthetic data, and said training cell type fraction matrix is a synthetic matrix.

15. The method according to claim 14, comprising: receiving a non-synthetic training bulk gene expression data and a non-synthetic training cell type fraction matrix; retraining said second network model by applying unsupervised learning thereto to infer said non-synthetic training cell type fraction matrix from said non-synthetic training bulk gene expression data, thereby provide a retrained network model; and redefining said second trained network model as said retrained network model.

16. The method according to any of claims 14 and 15, comprising synthesizing said training bulk gene expression data and said training cell type fraction matrix.

17. The method according to claim 16, wherein synthesizing said training cell type fraction matrix comprises generating a set of random vectors according to a multivariate probability distribution, and constructing said training cell type fraction matrix using said set as matrix columns or matrix rows.

18. The method according to claim 17, wherein said multivariate probability distribution is selected from the group consisting of a Dirichlet distribution, a multinomial distribution, a multivariate normal distribution, multivariate logistic distribution, Wishart distribution, a multivariate hypergeometric distribution, a copula distribution, and a multivariate Poisson distribution.

19. The method according to any of claims 16-18, comprising receiving an input cell type signature matrix, wherein said synthesizing said training bulk gene expression data comprises multiplying said synthetic matrix by said input cell type signature matrix.

20. A computer software product, comprising a computer-readable medium in which program instructions are stored, which instructions, when read by a computer, cause the computer to receive a bulk gene expression data of the biological sample and to execute the method according to any of claims 11-19.

21. A system for training an artificial neural network to infer cell-type abundance in a biological sample, the system comprising:an input circuit for receiving a non- synthetic training bulk gene expression data and a nonsynthetic training cell type fraction matrix; and a hardware processor configured to execute the method according to any of claims 11-19.

Images