Accurate and robust information deconvolution from bulk tissue transcriptomes

Abstract

The present disclosure relates to methods for deconvolving bulk or spatial RNA sequencing data, computer-readable media storing processor-executable instructions adapted to deconvolve bulk or spatial RNA sequencing data, and systems for deconvolving bulk or spatial RNA sequencing data to characterize cell type composition.

Description

[Technical field]

[0001] The present disclosure relates generally to sequencing technology. More specifically, the present disclosure relates to methods for deconvoluting sequencing data of mixed samples, and to associated apparatus and associated computer-readable media storing processor-executable instructions for obtaining mixture compositional characteristics. [Background technology]

[0002] In molecular biology, messenger RNA (mRNA) is a large group of RNA molecules that carry genetic information from deoxyribonucleic acid (DNA) to ribosomes, specifying the amino acid sequence of the protein product of gene expression. The RNA polymerase enzyme transcribes genes into the primary transcript mRNA (also known as precursor mRNA), which is then processed into mature mRNA. This mature mRNA is then translated into polymers of amino acids: proteins. This is as summarized in the central dogma of molecular biology:

[0003] As in DNA, the genetic information in mRNA is encoded in a sequence of nucleotides, arranged into codons, each of which consists of three base pairs. RNA sequencing (RNA-seq) is the process of determining the sequence of nucleotides in a strand of RNA. Each codon encodes a specific amino acid, except for stop codons (which end protein synthesis). This process of translating codons into amino acids requires two other types of RNA: transfer RNA (tRNA), which mediates the recognition of codons to yield the corresponding amino acid, and ribosomal RNA (rRNA), a central component of the ribosome's protein-making machinery.

[0004] Bulk tissue RNA-seq is a widely adopted method applied to understand genome-wide transcriptome changes in various states (e.g., normal or diseased). Since bulk tissues often consist of various cell types, bulk RNA-seq measures the average expression of each gene. This is the sum of cell type-specific gene expression weighted by cell type proportions. Knowing the cell type composition and its proportions in intact tissues is important to understand the biology of tissues. This can lead to characterizing tissue states from differences in cell composition, and adjusting for these changes can better direct downstream analyses. However, bulk tissue RNA-seq data does not directly provide information on cell type composition because the gene expression levels of each cell type mixed in the bulk data are not clear.

[0005] Recent breakthroughs in spatial transcriptomics have enabled characterization of gene expression across the transcriptome at spatially resolved locations within tissue sections. However, reaching single-cell resolution while measuring tens of thousands of genes per transcriptome remains challenging. Some widely used techniques can achieve a resolution of 50–100 μm (equivalent to 3–30 cells depending on tissue type), where transcripts may originate from one or more cell types. Unlike bulk RNA-seq, profiling data at each spot contains significant dropouts because only a small number of cells are sequenced, posing additional challenges to uncover cell type content. Both bulk RNA-seq and spatial transcriptomics data at multi-cell resolution are complex RNA-seq data, and the cell type mix ratio is unknown.

[0006] The rapid development of single-cell RNA-seq (scRNA-seq) technology has enabled cell-type-specific transcriptome profiling. Although such technology can provide information missing from composite RNA-seq data, as cell type composition and ratios can be obtained directly from scRNA-seq data, the technology has low sensitivity and generates unacceptably large noise. This is due to high dropout rates and cell-to-cell variability. As a result, scRNA-seq technology requires a large number of cells (thousands to tens of thousands) to ensure statistical significance of the results. In addition, cells must remain viable during capture. These requirements make scRNA-seq technology expensive, which limits its application in clinical studies with many subjects, and does not allow tissue dissociation and cell capture in real time. Furthermore, scRNA-seq technology is not well suited to characterize cell type ratios within solid tissues, because the cell dissociation and capture steps can be biased toward certain cell types.

[0007] Sequencing at the single-cell level is not always feasible and has its own limitations as we will discuss. Furthermore, there is also a lot of existing bulk RNA-seq data that could benefit from information obtained from cell type composition. Therefore, computational approaches have been developed to deconvolve cell type proportions from bulk tissue RNA-seq data. The deconvolution process is essentially an optimization problem, where the mixture proportions of a finite number of cell types are the parameters to be optimized. The goal is to minimize the difference between the observed gene expression in the bulk tissue RNA-seq data and its corresponding expected value (calculated as the sum of the given cell type-specific expression weighted by the mixture proportion parameters). The best mixture proportion that minimizes the difference is the final output.

[0008] One such computational method is disclosed in Wang et al., "Bulk tissue cell type deconvolution with multi-subject single-cell expression reference," Nature Communications (published online January 22, 2019). The authors introduce the "multi-subject single-cell deconvolution" (MuSiC) method (code available), which uses cross-subject scRNA-seq to estimate cell type proportions in bulk RNA-seq data. More specifically, MuSiC is a weighted non-negative least squares regression (W-NNLS) and does not require preselected marker genes. MuSiC uses cross-subject changes that reflect gene stability to weight genes. Through an iterative estimation procedure, it automatically imposes larger weights on stable genes and smaller weights on fluctuating genes. Because this is a linear regression-based method, genes that show large cross-subject changes have low influence and therefore have a small impact on the regression. On the other hand, the most influential genes are those with high stability weights. MuSiC is one of many alternative computational methods available.

[0009] In addition, most methods limit data to a given set of cell type-specific genes, and their output varies with different selections of such gene sets, reducing the objectivity and robustness of the results. For example, CIBERSORT is well-known, but has been published as lacking sensitivity (see the World Wide Web at nature.com / articles / s41467-018-08023-x). Moreover, most existing methods are only suitable for relatively simple applications, e.g., peripheral blood mononuclear cells (PBMCs) and pancreas, where only a handful of cell types need to be considered or the differences between cell types are quite large. Their performance in complex tissues with dozens of different cell types or subtle differences in cell subtypes is questionable.

[0010] In view of the foregoing, there is a need for improved methods for accurate and robust deconvolution from bulk tissue and spatial transcriptomes at multi-cell resolution. Summary of the Invention

[0011] This disclosure provides methods (including computer-implemented methods), computer programs, computer systems, and devices for deconvolution of bulk RNA sequencing data, with the goal of satisfying the need for obtaining accurate and robust cell type fraction estimates from bulk tissue transcriptomes.

[0012] The present disclosure provides a method for deconvolving bulk RNA sequencing data with predetermined cell type specific expression obtained from single cell RNA-seq of cell types associated with a bulk tissue, the method comprising: i) selecting from the single cell RNA-seq data a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) calculating cell type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data from the single cell RNA-seq data using cell type annotations; and iii) calculating from the single cell RNA-seq data a subset of most variably expressed genes for each cell type from the count-based sequencing data matrix, the subset of most variably expressed genes, and the cell type annotations. The method includes any one or more of: fitting a cross-sample distribution for each of the most variably expressed genes and a bulk distribution for each of the most variably expressed genes from the normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; v) defining a loss function between the bulk distribution and the mixed single-cell distribution; and vi) applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. In some embodiments, the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix.In some embodiments, the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix. In some embodiments, the cross-sample distribution for each cell type and for each subset of the most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution. In some embodiments, the bulk distribution for each subset of the most variably expressed genes from the normalized bulk matrix, and for the subset of the most variably expressed genes, is a bulk Gaussian distribution.

[0013] The disclosure also provides a method for deconvolution of bulk RNA sequencing data, comprising six exemplary steps: i) taking input from three sources (bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations) and selecting a subset of the most variably expressed genes from a matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; and ii) estimating mean and variance parameters of expression per gene per cell type. iii) calculating cross-cell type specificity for genes; iv) estimating cross-sample gene variability from the composite data or single cell samples depending on multiple sample availability; v) estimating scaling coefficients per gene using both the composite data and single cell data; and vi) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in bulk or spatial RNA sequencing data, thereby providing a method for inferring cell type percentages in bulk RNA sequencing data.

[0014] The present disclosure also provides a computer-readable medium storing processor-executable instructions adapted to cause one or more computing devices to deconvolve bulk RNA sequencing data according to any of the methods described herein.

[0015] The disclosure also provides a system that includes one or more processors and a memory having processor-executable instructions that, when executed by the one or more processors, cause an apparatus to deconvolve bulk RNA sequencing data according to any of the methods described herein.

[0016] It is to be understood that both the foregoing general description and the following detailed description are exemplary of the present disclosure and are not restrictive.

[0017] The disclosure is best understood when the following description is read in conjunction with the accompanying drawings. The patent or application document contains at least one drawing figure executed in color. Copies of this patent or patent application publication with color drawings may be obtained from the U.S. Patent and Trademark Office upon request and payment of the necessary fee. The drawings include the following figures: [Brief description of the drawings]

[0018] [Figure 1] FIG. 1 shows an overview of the workflow of one embodiment of the disclosed method. [Diagram 2] A, B, C illustrate three different hypothetical gene expression versus cell type patterns as a basis for selecting the most informative genes in one embodiment of the disclosed method. [Diagram 3] FIG. 1 illustrates the calculation of gene variance across all cell types and selection of the top 2,500 variable genes in one embodiment of the disclosed method. [Figure 4] FIG. 1 illustrates the calculation of overall or global mean variance within cell type mean variance in one embodiment of the disclosed method. [Diagram 5]FIG. 13 illustrates the estimation of cell type specific variance and mean by fitting Gaussian distributions in one embodiment of the disclosed method. [Figure 6] FIG. 1 illustrates estimation of bulk data cross-sample variance and mean by fitting a Gaussian distribution in one embodiment of the disclosed method. [Figure 7] FIG. 1 illustrates a comparison between a mixture distribution of single cell data and a distribution of bulk cell data in one embodiment of the disclosed method. [Figure 8] FIG. 2 illustrates the results of fitting a Gaussian distribution following application of step 3 of the embodiment of the disclosed method shown in FIG. 1 for an illustrative example. [Figure 9] FIG. 2 illustrates the results of fitting a Gaussian distribution following application of step 4 of the embodiment of the disclosed method shown in FIG. 1 for an illustrative example. [Figure 10] FIG. 2 shows, in an illustrative example, the results of applying step 5 of one embodiment of the disclosed method shown in FIG. 1, i.e., defining a model using ratio parameters, weights, and distributions for each gene learned from single cell and example bulk cell data. [Figure 11] FIG. 1 shows an outline of the workflow of another embodiment of the disclosed method (the AdRoit method). [Figure 12A] FIG. 12 illustrates two options for selecting the most informative genes during the first step of the method disclosed in FIG. 11. [Figure 12B] FIG. 12 shows a hypothetical example illustrating the types of cells that may be selected using the method disclosed in FIG. [Figure 13] FIG. 12 illustrates the second step of the method disclosed in FIG. 11, which estimates mean and variance parameters by fitting a negative binomial distribution to each gene in a cell type. [Figure 14] FIG. 12 shows a hypothetical example demonstrating the effect of per-gene scaling factors applied during the fifth step of the method disclosed in FIG. 11. [Figure 15A]Summary of human islet cell composition in 18 subjects. [Figure 15B] t-SNE graph showing that the four cell types are distinct from each other. [Figure 16A] 1 is a graph reflecting a comparison of the estimation accuracy against the true percentage in the AdRoit method's estimation for all cell types from 18 subjects. [Figure 16B] 13 is a graph reflecting a comparison of the estimation accuracy against the true percentage in the MuSiC method's estimation for all cell types from 18 subjects. [Figure 16C] 1 is a graph reflecting a comparison of the estimation accuracy against the true percentage in the NNLS method's estimates for all cell types from 18 subjects. [Figure 17] FIG. 16 is a table listing four separate statistical measures (mAD, RMSD, Spearman and Pearson correlation) calculated for each of the three graphs in FIGS. 16A, 16B, and 16C. [Figure 18A] Summary of human trabecular meshwork cell composition of eight donors. [Figure 18B] t-SNE graphs showing differences as well as similarities between cell types. The data was used to evaluate the disclosed method against other conventional methods. [Figure 19A] 1 is a graph reflecting a comparison of the estimated precision to the true percentage in the AdRoit method results for eight donors. [Figure 19B] 1 is a graph reflecting a comparison of estimated precision to true percentage in MuSiC method results for eight donors. [Figure 19C] 1 is a graph reflecting a comparison of the estimated precision to the true percentage in the NNLS method results for eight donors. [Figure 20] FIG. 19 is a table listing four separate statistical measures (mAD, RMSD, and Spearman and Pearson correlations) calculated for each of the three graphs in FIGS. 19A, 19B, and 19C. [Figure 21]A comparison of how much the estimates deviate from the truth among the three methods. Each dot represents a donor and each row is a cell type in the human trabecular meshwork. [Figure 22] FIG. 13 reflects estimated and true data calculated using both the AdRoit and MuSiC methods for human trabecular meshwork cell types. [Figure 23] FIG. 1 shows a receiver operating characteristic (ROC) curve showing that the area under the curve (AUC) for detecting human trabecular meshwork cell types is significantly higher for the AdRoit method than for the MuSiC method, indicating that AdRoit is more sensitive. [Figure 24A] Summary of the cellular composition of five mice. [Figure 24B] 1 is a t-SNE plot of cell types found in the mouse dorsal root ganglion single cell data used, which was later used to evaluate the disclosed method against other conventional methods. [Figure 25A] 1 is a graph reflecting a comparison of the estimated precision versus true cell percentage in the AdRoit results for five mice. [Figure 25B] 13 is a graph reflecting a comparison of the estimated precision versus true cell percentage in the MuSiC method results for five mice. [Figure 25C] 1 is a graph reflecting a comparison of the estimated precision versus true cell percentage in the NNLS method results for five mice. [Figure 26] FIG. 13 is a graphical representation comparing the results of the AdRoit, MuSiC, and NNLS methods on mouse data, using mAD, RMSD, Pearson and Spearman correlation as statistical measures. [Figure 27] FIG. 13 shows that AdRoit method-based estimates of cell type percentages for real human islet bulk RNA-seq data are highly reproducible for repeated samples from the same donor. [Figure 28]FIG. 1 shows that cell type percentages of human islet data estimated using the Adroit method are consistent with RNA-Fish measurements of cell type percentages. [Figure 29] FIG. 1 shows that beta cell fraction estimated using the Adroit method has a significant negative linear relationship with donor HbA1C levels (including both healthy and T2D cells). [Diagram 30] FIG. 1 shows that beta cell fraction estimated using the Adroit method in T2D patients is significantly lower than in healthy subjects. [Diagram 31] FIG. 11 compares the estimations achieved by the Stereoscopic and AdRoit methods on a simulated spatial spot containing five different PEP cell subtypes. [Diagram 32] Figure 16. Comparison of performance when the percentage of cells is low using simulated data. A set of low percentage PEP cells was simulated and mixed with two other PEP cell types. The results show that the median estimate achieved using the AdRoit method was closer to the true ratio than the estimate achieved using the stereoscope. [Diagram 33] FIG. 10 compares the detection rates of AdRoit and stereoscopy using simulated spatial spots. The simulation involves six different mixtures of cell types. Each type of mixture includes a set of low percentage cell types. The evaluation is to see how much of the low percentage cell types are detected in each given low percentage. [Diagram 34] FIG. 1 illustrates the cell type content estimated by the AdRoit method in each spatial spot of a mouse brain coronal tissue section. [Diagram 35] FIG. 35 is an ISH image of Wfs1, Prox2, and Rarres2 genes from the Allen mouse brain atlas verifying that the cell type locations shown are accurate. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0019] Various terms relating to the aspects of the present disclosure are used throughout the specification and claims. Such terms should be given their ordinary meaning in the art unless otherwise indicated. Other specifically defined terms should be interpreted consistent with the definitions given herein.

[0020] Unless otherwise expressly stated, it is never intended that any method or aspect described be construed as requiring that its steps be performed in a particular order. Thus, in the claims or description, if a method claim does not specifically state that the steps are to be limited to a particular order, no order is ever intended to be inferred in any sense. This holds true against any possible unexpressed basis for interpretation, including questions of logic regarding the placement or operational flow of the steps, the apparent meaning derived from grammatical construction or punctuation, or the number or type of aspects described in the specification.

[0021] RNA sequencing technology may provide unprecedented opportunities to learn disease mechanisms and discover new therapeutic targets. Recent spatial transcriptomic methods further enable transcriptome profiling in spatially resolved spots within tissue sections. In controlled experiments, it is often crucial to know the variability of cell composition under therapeutic intervention. Understanding the cell type content in each tissue spot is also important for the interpretation of spatial transcriptomic data. Single-cell RNA-seq has the potential to reveal cell type composition and expression heterogeneity in different cells, but remains costly and is often infeasible when live cell acquisition or sufficient dissociation is not possible. To leverage bulk and spatial RNA-seq data when sequencing at the single-cell level is not feasible, presented herein are methods to estimate the proportion of each cell type in bulk or spatial RNA-seq data using known single-cell sequence data of relevant cell types (e.g., data available in the public domain). The methods described herein are more accurate and robust because they jointly model per-gene technique bias, cell type specificity of genes, and cross-sample variability. Systematic benchmarking evaluation has demonstrated superior sensitivity and specificity to other existing methods, even in neuronal cells, where many closely related subtypes exist.

[0022] The method disclosed herein provides a statistical method to estimate the proportion of each cell type in bulk RNA-seq data using independently obtained expression profiles (e.g., single-cell data) of related cell types (often published) obtained from count-based sequencing techniques. The method is particularly well suited to detect rare (proportion less than about 5%) cell types. One of the assumptions in carrying out the method described herein is that the tissue used for single-cell RNA-seq contains the same or no fewer cell types than those in the bulk or spatial sequencing sample.

[0023] As used herein, the term "about" means that the numerical values ​​described are approximate and that small variations would not significantly affect the practice of the disclosed embodiments. When numerical values ​​are used, unless otherwise indicated by context, the term "about" means that the numerical values ​​can vary by ±10% and still remain within the scope of the disclosed embodiments.

[0024] As used herein, the term "about" means that the numerical values ​​described are approximate and that small variations would not significantly affect the practice of the disclosed embodiments. When numerical values ​​are used, unless otherwise indicated by context, the term "about" means that the numerical values ​​can vary by ±10% and still remain within the scope of the disclosed embodiments.

[0025] As used herein, the term "comprising" may be replaced with "consisting of" or "consisting essentially of," as appropriate in certain embodiments.

[0026] The disclosed method, device and computer readable medium aim to estimate the proportion of cell types from bulk tissue transcriptome accurately and robustly.Use existing count-based sequencing data, such as single-cell RNA sequencing data, as a reference with cell identity annotation.Find out the unique distribution characteristics of highly informative genes for each cell type, and estimate the composition of cell types from bulk tissue or spatial RNA sequencing data without relying on marker selection. Critical to the success of the disclosed methods, apparatus, and computer readable media are: 1) considering the entire distribution of gene expression values, or the mean and variance parameters that define the distribution, when estimating mixture proportions, and not just the mean, 2) giving higher weight to genes that are more distinguishable between cell types, i.e., genes with expression that is highly unique to a particular cell type, 3) giving lower weight to genes that are highly variable across multiple samples, 4) using an adaptive learning approach to estimate per-gene scaling factors to account for platform differences between bulk or spatial RNA sequencing data and single-cell RNA sequencing data, and 5) including regularization terms in the model to minimize the effects of statistical collinearity. These five features, among others, of the disclosed methods, apparatus, and computer readable media combine to provide improvements over existing methods.

[0027] The present disclosure provides a method for deconvolution of bulk RNA sequencing data. In some embodiments, the method includes six exemplary steps: i) selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) calculating cell type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using cell type annotations; and iii) calculating cell type specific weights for each cell type and for each subset of most variably expressed genes, the subset of most variably expressed genes, and the cell type annotations from the count-based sequencing data matrix. The method includes any one or more of the following steps: fitting cross-sample distributions to each subset of most variably expressed genes and the subset of most variably expressed genes from the normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; v) defining a loss function between the bulk distribution and the mixed single-cell distribution; and vi) applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. These steps are illustrated in FIG. 1 (showing single-cell RNA sequencing as count-based sequencing). FIG. 1 shows a schematic of the workflow of one embodiment of the disclosed method. Each step is described in turn below with reference to the inputs, outputs, and purpose or rationale for each step. Each of these processing steps can be performed by a computing device (e.g., a computer). In some embodiments, all process steps are performed by a computer.

[0028] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, third, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fifth step. In some embodiments, the method includes a fifth step and one or more of the first, second, third, fourth, and sixth steps, or any combination of these additional steps. In some embodiments, the method includes a sixth step. In some embodiments, the method includes a sixth step and one or more of the first, second, third, fourth, and fifth steps, or any combination of these additional steps.

[0029] In some embodiments, the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix. In some embodiments, the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix. In some embodiments, the cross-sample distribution for each cell type and for each subset of most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution. In some embodiments, the bulk distribution for each subset of most variably expressed genes and the subset of most variably expressed genes from the normalized bulk matrix is ​​a bulk Gaussian distribution. The methods described herein result in the inference of single-cell distribution ratios for bulk RNA sequencing data.

[0030] In some embodiments, the method further comprises forming a matrix of count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells and normalizing the matrix. In some embodiments, the method further comprises forming a bulk matrix of bulk RNA sequencing counts and normalizing the bulk matrix. In some embodiments, the method further comprises forming a matrix of count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells and normalizing the matrix, and forming a bulk matrix of bulk RNA sequencing counts and normalizing the bulk matrix. In some embodiments, the method further comprises obtaining a cell type annotation. In some embodiments, the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix. In some embodiments, the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix.

[0031] In some embodiments, the method further comprises determining a proportion of RNA from each cell type for which the bulk RNA sequencing data was obtained. In some embodiments, the method further comprises determining a proportion of each cell type for which the bulk RNA sequencing data was obtained. In some embodiments, the method further comprises determining a proportion of RNA from each cell type for which the bulk RNA sequencing data was obtained, and determining a proportion of each cell type for which the bulk RNA sequencing data was obtained.

[0032] Step 1: Select the top "N" highly variable genes The goal of the first step in an exemplary embodiment of the disclosed method (FIG. 1) is to select the most informative genes. We apply this step to single-cell RNA-seq (scRNA-seq) data, but it can be applied to any count-based sequencing data as described herein. In scRNA-seq data, not all genes in the genome are informative for the identity of a particular cell. One reason is that many genes are essential for cell survival and other fundamental cellular functions. These genes are ubiquitously expressed across many different cell types and cannot be biologically distinguished by cell type. Furthermore, due to technical limitations, scRNA-seq typically captures hundreds to thousands of genes from each cell, depending on the platform (10x or C1). The majority of these genes have no readout across all cells. In either case, these genes do not contribute to identifying the cell identity and are therefore not informative in the deconvolution of the bulk data.

[0033] One reason why selecting the most informative genes is important is illustrated in Figures 2A, 2B, and 2C with graphs of expression (ordinate) versus cell type C1, C2, C3, C4, and C5 (abscissa) for each of three hypothetical genes. Figure 2A shows gene 1 as informative because the data within each cell type is relatively consistent and the data can be distinguished among the five cell types. The data in Figure 2B is not useful because the data within each cell type is too variable. The data in Figure 2C is not useful because the data between cell types is not different enough.

[0034] Known analytical methods can be used to select highly variable genes. See, for example, A. Butler et al., "Integrating single-cell transcriptomic data across different conditions, technologies, and species," Nat. Biotechnol. (2018) (A. Butler, Nat. Biotechnol.); and F. Wolf et al., "SCANPY: Large-scale single-cell gene expression data analysis," Genome Biol. (2018). Usually, when single cells are computationally clustered based on scRNA-seq data, the top 2,000 highly variable genes provide a good separation between different cell types. However, it is recommended to select some number more than this 2,000 number of genes, because data processing may induce information loss. On the other hand, a balance must be maintained, because selecting too many genes introduces noise. Thus, in some embodiments, the top 2,500 highly variable genes are selected. More or less genes than that number can be selected depending on the application (e.g., cell type). The preferred number of variable genes to select can be pre-determined by trial and error based on which number provides the best validation. "Predetermined" means pre-determined, and thus the pre-determined characteristics must be determined, i.e., selected or at least known, prior to any event. Preferably, the range of minimum and maximum number of highly variable genes to select is about 1,000 to about 5,000. Genes are selected from the entire transcriptome measurable by RNA-seq technology. There are about 25,000 genes in the human transcriptome. There are about 20,000 genes in the mouse transcriptome.

[0035] Due to well-known variance effects in RNA-seq data, directly calculating the change from the count matrix may overestimate the variance. The method described herein addresses such overestimation by calculating the variance from a variance-stabilizing transformation (VST) data matrix, and selects genes based on the rank of these variances. Figure 3 illustrates a typical calculation of the variance of genes across all cell types and the selection of the top 2,500 highly variable genes. The algorithm for this procedure is easily programmed in the "Seurat" R package disclosed in A. Butler, Nat. Biotechnol. For example, the function "FindVariableFeatures" is used to select the top 2,500 highly variable genes. Of course, other algorithms can be used for the selection of the top 2,500 highly variable genes.

[0036] In this first step, a single-cell expression matrix constitutes the input with rows representing genes and columns representing individual cell types, as exemplified in FIG. 3. It is recommended, but not required, to use a unique molecular identifier (UMI) count matrix (data from 10x platform) or RPKM (data from C1 platform). Cell type annotation is not required in this step. After applying VST to the single-cell count matrix, the standard deviation (represented by the Greek symbol sigma or "σ") is calculated for each row (gene) to obtain the 2,500 most variable genes. Standard deviation is a measure of the degree of deviation relative to the entire group of data. Standard deviation is calculated as follows: 1) calculate the average or mean; 2) for each number, subtract the mean and square the result; 3) calculate the average (variance) of the squared differences; 4) calculate the square root of that average. The output from the first step of the disclosed method is the top "N" number (i.e., 2,500) of highly variable genes. The disclosed method later restricts the calculations to these N genes.

[0037] Step 2: Calculating cell type-specific weights The input to the second step in the exemplary embodiment of the disclosed method is the same single cell count matrix as in the first step, but can be any count-based sequencing data matrix as described herein. However, the second step also requires cell identity information (i.e., cell type annotation) as input, because it calculates cell type-specific variance. One of the goals of the second step is to quantify the importance of genes in defining cell types.

[0038] Figure 4 illustrates a representative calculation of the overall or global average variance and the within-cell-type average variance. The average variance is calculated across cells within each cell type and compared to the overall average variance. For the same variance reasons, the variance over the log counts (zero counts plus one) is calculated. If a gene is stably expressed within a cell type (i.e., meaning the average variance is low) but highly variable across all cells (i.e., meaning the average variance is high), then that gene should be a good specifier for that cell type and should therefore receive a large weight.

[0039] To formally specify the weight ("W") of each cell type, let "N" be the total number of cells and let "n k Let N be the number of cells in cell type k. Thus, N=n1+n2+n3+.... For a particular gene i, σ i 2 Let be the variance over all cells, and σ ik 2 Let be the variance of that gene within cell type k. Thus, the weight for a particular gene and cell type is expressed as

[0040]

number

[0041] Step 3: Fitting cell-type specific Gaussians across subjects The input to the third step in an exemplary embodiment of the disclosed method includes the single cell count matrix and the list of highly variable genes and cell type annotations from step 1, but can be any count-based sequencing data matrix as described herein. For multi-sample single cell data, sample information must also be input.

[0042] A statistical test analyzes a particular set of data to arrive at more general conclusions. There are several approaches to doing this, but the most common is based on the assumption that the data in a population has a particular continuous probability distribution. The most commonly used distribution is the bell-shaped Gaussian distribution (also known as the normal distribution). The normal distribution is often used in the natural and social sciences to represent real-valued random variables whose distribution is unknown. Random variables with a Gaussian distribution are said to be normally distributed and are said to have normal deviations.

[0043] One feature of the disclosed method is that it uses a full distribution when estimating the mixture proportions. The distribution is obtained by fitting the normalized counts to a distribution (e.g., a Gaussian distribution) to estimate the variance and mean for each gene. The "normalizing" process involves adjusting values ​​measured at different scales to a common scale (i.e., eliminating the units of measurement), usually before averaging. Figure 5 shows how cell type-specific variance and mean are estimated by fitting a Gaussian distribution.

[0044] The disclosed method can estimate distribution (e.g., Gaussian distribution) using at least two methods, depending on whether multiple samples are available. When multiple samples are available, cells are pooled within cell types by adding read counts to form megacells for each cell type. Megacells reduce data sparsity and sampling variation due to technical limitations, and therefore better represent the unique transcriptome profile of each specific cell type. Unfortunately, multiple samples are not always available. When multiple samples are not available, the disclosed method estimates variance by randomly dividing cells into multiple subpopulations, pooling cells within each subpopulation, and using cells as if they were from different samples.

[0045] The disclosed method then normalizes the megacell count matrix for each sample. The method essentially follows the standard method for normalizing RNA-seq data as disclosed in A. Butler, Nat. Biotechnol., and M. Love et al., "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2," Genome Biol. (2014). Here, the data is rescaled so that each megacell has the same total number of reads, and the rescaled data is log-transformed. However, the disclosed method makes some relatively small adjustments to the standard method to achieve better results. First, the log transformation ensures an approximation to the distribution (e.g., Gaussian), but greatly reduces the variance on which the disclosed method is based. To mitigate the downside, the total number of reads is scaled to a large number (e.g., 10 7) Second, we typically add 1 to zeros in the matrix to avoid logarithmic errors, but at the same time, can unfavorably alter the distribution, especially for low-expressing genes. Adding 1 gives satisfactory results, but adding a smaller number (e.g., 0.1) helps improve the accuracy of the estimation and also speeds up the convergence of the algorithm described in step 5 of the method (below).

[0046] The output from the third step of the disclosed method includes a normalized expression matrix for the megacell and an estimated mean and variance for each selected gene. In the illustrated example, the output includes five sets of 2,500 Gaussian curves, one for each cell type. At this point, the disclosed method has completed processing the single cell data and now moves on to the bulk data.

[0047] Step 4: Fitting a Gaussian to the bulk data The list of selected genes from step 1 is combined with the multi-sample bulk RNA-seq count matrix to form the input to the fourth step of an exemplary embodiment of the disclosed method. Step 4 is very similar to step 3. First, the total number of reads is rescaled to the same number as in the single cell analysis (e.g., 10 7 ), and then add a small number (e.g., 0.1) to the zero counts to perform a log transformation. In theory, rescaling all reads for the bulk RNA-seq data is unnecessary and gives similar results as unscaled data. In some embodiments, the disclosed method includes rescaling for the practical reason that algorithms converge faster when the numbers are closer to the sum of single cells. Figure 6 shows how bulk data cross-sample variance and mean are estimated by fitting Gaussian distributions. The output from the fourth step of the disclosed method includes the normalized expression matrix and estimated mean and variance across samples for each selected gene.

[0048] Step 5: Specify the loss function The fifth step in an exemplary embodiment of the disclosed method takes as input all the outputs from the previous steps, i.e., the top highly variable genes, the cell type specific weights per gene, the normalization matrix for count-based sequencing data (e.g., single cell RNA sequencing and bulk data), and the distribution mean (e.g., Gaussian mean) and variance estimates for count-based sequencing (e.g., single cell RNA sequencing and bulk data). For machine learning problems, choosing the right loss function is important for parameter estimation. Loss functions are central to modern machine learning. Loss functions move algorithms from theory to practice, transforming neural networks from spurious matrix multiplication to deep learning. At their core, loss functions are simple. They are a way to evaluate how well an algorithm models a dataset. If the prediction is completely off, the loss function outputs a large number. If the prediction is good, the loss function outputs a small number. As we modify parts of the algorithm in an attempt to improve the model, the loss function tells us whether the modifications are likely to be successful or not.

[0049] Consider some of the most common loss functions in use today, from simple to complex. Mean Squared Error (MSE) is a basic loss function that is easy to understand and implement and generally performs reasonably well. To calculate MSE, the difference between the prediction and the ground truth is calculated, squared, and averaged over the entire dataset. Another loss function, the likelihood function, is also relatively simple and is widely used in classification problems. The likelihood function takes the predicted probabilities for each input example and multiplies them together. While the output is not human-interpretable, the likelihood function is useful for comparing models. Log-loss is also a loss function often used in classification problems, and is a logarithmic modification of the likelihood function. Loss functions do more than just give a static representation of how the model is performing. They tell us something about how the algorithm fits the data in the first place. Most machine learning algorithms use some kind of loss function during the optimization process, or in the process of finding the best parameters (weights) for a dataset.

[0050] The algorithm of the disclosed method is designed to find the best set of ratio parameters by minimizing the difference between the mixture distribution of single cell data and that of bulk cell data. Figure 7 illustrates a comparison between the two distributions. One of the goals is to minimize the difference between the sum of single cell data and the bulk cell data. Therefore, in some embodiments, the disclosed method uses the Kullback-Leibler (KL) divergence as its loss function (see S. Kullback & R. Leibler, On Information and Sufficiency (Ann. Math. Stat. 1951)). The KL divergence is particularly suitable for quantifying the similarity between two distributions. Let f1(x) and f2(x) be two probability density functions for a continuous variable X. We define the KL divergence between the two as follows:

[0051]

number

[0052]

number

[0053]

number

number

number

[0054]

number

number

[0055]

number

[0056] It is possible, but μ' and σ 2 Parameterizing μ′ would unnecessarily complicate the model, since the goal of the disclosed method is to estimate θ′. The algorithm finds a set of global rate parameters θ′ that fit across all selected genes. Considering that the estimation errors across all genes are randomized and have negligible global influence on the estimation of θ′, μ′ and σ from bulk and single-cell data can be easily modeled. 2 A rough estimate of θ′ is sufficient for the algorithm to find the best θ′. Therefore, the disclosed model uses the estimated μ′ and σ 2 We use θ′ directly and calculate the probability density treating them as known parameters. θ′ is the only unknown parameter to estimate in the model.

[0057] Step 6: Model estimation The prescribed loss function between the actual bulk data distribution and the mixture distribution of single cells output from the fifth step of the disclosed method is the input to the sixth step. The method employs gradient descent to estimate the ratio parameters. Gradient descent is a first-order iterative optimization algorithm for finding the minimum of a function. To find a local minimum of a function using gradient descent, we take steps proportional to the negative of the gradient (or the approximate gradient) of the function at the current point. If instead we take steps proportional to the positive of the gradient, we approach a local maximum of the function. In that case, the procedure is known as gradient ascent. Gradient descent was first proposed by M. Cauchy in 1847 (see M. Cauchy, Methode Generale Pour la Resolution des Systemes D'equations Simultanees (Ubersetzt von Richard Pulskamp 2010)). Compte rendu des seances l'academie des Sci. (1847). See also D. Bertsekas, Nonlinear Programming (2d Athena Scientific 1999).

[0058] The disclosed method first calculates each θ k We derive the gradient (denoted as G) by taking the derivative of the loss function with respect to . To simplify the inference, we decompose the KL divergence into entropy and cross-entropy terms.

[0059]

number

[0060]

number

[0061]

number

[0062] To perform gradient descent, the method uses a non-negative θ (0) At the same time, we randomly initialize all θ k Update it to look like this:

[0063]

number

[0064]

number

[0065] In some embodiments, the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix. In some embodiments, the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix.

[0066] In some embodiments, the cross-sample distribution for each cell type and for each subset of the most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution. In some embodiments, the bulk distribution for each subset of the most variably expressed genes and the subset of most variably expressed genes from the normalized bulk matrix is ​​a bulk Gaussian distribution.

[0067] In some embodiments, the selecting step includes calculating a standard deviation for each gene in the plurality of genes, determining a threshold number of standard deviations, and selecting a subset of the most variably expressed genes having standard deviations above the threshold number. In some embodiments, the calculating the cell type specific weights includes comparing the overall mean variance to the cell type mean variance for each of the fixed number of cells. In some embodiments, the fitting step includes using the entire distribution when estimating the mixture proportions. In some embodiments, the fitting step further includes obtaining a distribution by fitting the normalized counts to the distribution to estimate the variance and mean for each gene. In some embodiments, the distribution is a Gaussian distribution. In some embodiments, the defining the loss function includes applying the Kullback-Leibler divergence. In some embodiments, the applying the loss function includes employing a gradient descent method.

[0068] The present disclosure also provides another embodiment of the AdRoit method for deconvolution of bulk RNA sequencing data. The disclosed method embodiment described herein may be described as an "Accurate and Robust Method for Inference of Transcriptome Composition" and may be identified by the acronym "AdRoit". The AdRoit method aims to accurately and robustly estimate cell type proportions from composite transcriptome data including bulk RNA-seq and spatial transcriptome data. The method uses relevant pre-existing single-cell RNA-seq data with cell feature annotations as a reference, selects informative genes, and estimates the expression mean and variance of the selected genes per cell type. Furthermore, in one embodiment, the AdRoit method calculates the variability per gene across samples as well as their cell type specificity, whereby the loss function of each gene is weighted differently in the model. Furthermore, the AdRoit method calculates a scaling coefficient per gene to minimize the technique difference between the single cell and the target composite data. Together, the AdRoit method estimates cell type percentages by feeding them into a regularization model and optimizing a weighted sum of per-gene loss functions. Critical to the accuracy and robustness of the method are: 1) selecting the most informative genes to use in the deconvolution task, 2) a per-gene loss function that is appropriately weighted by how specifically one cell type can be made to differ from others and by how stable its expression is across multiple samples, 3) a per-gene scaling factor that normalizes gene expression values ​​from different sequencing platforms (e.g., TPM or read counts from bulk RNA-seq, unique molecular identifiers (UMIs) from single-cell RNA-seq and spatial transcriptome sequencing), and 4) a regularization regression model that avoids collinearity between closely related cell types (e.g., subtypes).

[0069] In some embodiments, the AdRoit method for deconvolution of bulk or spatial RNA sequencing data includes the following exemplary steps: i) taking input from three sources (bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations) and selecting a subset of the most variably expressed genes from a matrix of count-based single-cell sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) estimating the mean and variance parameters of expression per gene per cell type; and iii) calculating the cross-cell type specificity of the genes. The method includes any one or more of the following steps: iv) within each cell type, for each gene, estimating cross-sample gene expression variability based on the average gene expression of multiple cells in each sample depending on multiple sample availability or forming multiple samples by subsampling cells from the same sample; v) estimating scaling factors per gene using both the composite data and the single cell data; and vi) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data. These steps are illustrated in FIG. 11. The figure shows a schematic of the workflow of one embodiment of the disclosed method. Each step is described in turn below with reference to the inputs, outputs, and purpose or rationale for each step. Each of these processing steps can be performed by a computing device (e.g., a computer). In some embodiments, all process steps are performed by a computer. Spatial transcriptomics is a special type of bulk sequencing where the number of cells is very small.

[0070] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, third, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fifth step. In some embodiments, the method includes a fifth step and one or more of the first, second, third, fourth, and sixth steps, or any combination of these additional steps. In some embodiments, the method includes a sixth step. In some embodiments, the method includes a sixth step and one or more of the first, second, third, fourth, and fifth steps, or any combination of these additional steps.

[0071] In some embodiments, the AdRoit method for deconvolution of bulk or spatial RNA sequencing data includes any one or more of the following exemplary steps: i) calculating the cross-cell type specificity of genes; ii) within each cell type, for each gene, estimating the cross-sample gene expression variability based on the average gene expression of multiple cells in each sample depending on the availability of multiple samples, or forming multiple samples by subsampling cells from the same sample; iii) estimating a scaling factor per gene using both the composite data and the single-cell data; and iv) building a weighted and regularized regression model using all known quantities and using the model to estimate the cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data. Each step is described below in turn with reference to the inputs, outputs, and objectives or rationale for each step. Each of these processing steps can be performed by a computing device (e.g., a computer). In some embodiments, all process steps are performed by a computer. Spatial transcriptomics is a special type of bulk sequencing where the number of cells is very small. In some embodiments, the calculation of cross-cell-type specificity of genes is based on estimated mean and variance parameters of expression per gene per cell type from a subset of the most variably expressed genes selected from a matrix of count-based single-cell sequencing data (obtained from three sources: i) bulk or spatial RNA-seq data, ii) single-cell RNA-seq data, and iii) cell type annotations), where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells.

[0072] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, and fourth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, or fourth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, and fourth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, and third step, or any combination of these further steps.

[0073] Step 1: Gene selection The objective of the first step in the second embodiment of the disclosed method is to select the most informative genes. We apply this step to single-cell RNA-seq (scRNA-seq) data, but it can be applied to any count-based sequencing data as described herein. The step starts by taking input from three sources: bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations. Thus, the input data is a single-cell UMI count matrix with cell type annotations associated with each cell. Each column of the matrix corresponds to a cell, and each row of the matrix corresponds to a gene. Each entry in the matrix is ​​a UMI count for a particular gene in a cell. The bulk data to be deconvolved can be transcripts per kilobase million (TPM) or read counts. Each row of the matrix is ​​a gene, and each column of the matrix is ​​a sample. The spatial transcriptome data to be deconvolved is also a UMI count matrix, but each column of the matrix is ​​a spatial spot, and each row of the matrix is ​​a gene. The mathematical notation and terminology used in describing the method steps are defined as follows.

[0074] i - gene index, i=1,...I k-cell type index, k=1,...,K k' - index of the cell type with the highest mean expression of gene i j - index of sample (bulk RNA-seq) or spatial spot n k - number of cells of cell type k X ik - collection of single-cell UMI counts of gene i for all cells of cell type k Y ij - count of gene i in bulk sample or spatial spot j λ ik - the variance parameter for gene i in cell type k p ik - Probability of gene i in cell type k getting one UMI λ ij- the dispersion parameter for gene i in the bulk sample or spatial spot j p ij - Probability for gene i of cell type in bulk sample or spatial spot j μ ik - average expression of gene i in cell type k

number

number

number

number

number

number

[0075] As illustrated in FIG. 12A, the method selects genes in one of two alternatives. The first option is to use a population of genes whose expression is enriched in each cell type in the single cell UMI count matrix. These genes are referred to as marker genes. The second option is to use a population of genes that vary the most among all cells in the single cell UMI count matrix. These genes are referred to as highly variable genes. In this second option, the variance for each gene is calculated after cell number balancing and variance stabilizing transform (VST) normalization, and the gene with the highest variance is selected. Either option provides equally accurate estimates.

[0076] To select marker genes, a predefined marker gene list can be input or a built-in tool can be used. The built-in tool takes as input the single cell UMI count matrix and cell type annotation. For each cell type, the tool calculates the fold change between the average UMI in that cell type and the average UMI in all other cell types, and ranks the genes by descending order of fold change. To unravel complex composite transcriptome data, it is sufficient to select about the top 200 genes from each cell type. Since some genes may mark multiple cell types, it is desirable to select markers that are present in 5 or fewer cell types to ensure specificity. Alternatively, it may be desirable to select markers that are present in a fixed number or fewer cell types or a fraction of the total number of cell types, whichever is smaller, to ensure specificity. To ensure accurate estimation, it is desirable to select a minimum of about 1,000 unique genes in total from the population of marker genes of all cell types.

[0077] Finding marker genes is often time consuming and can require extensive computational resources. However, even when marker genes are not readily available, the method can select highly variable genes. Usually, these genes are also informative for distinguishing cell types. To avoid the risk that selected highly variable genes may be dominated by large cell clusters while underrepresenting small clusters, the cell types in the single-cell UMI count matrix can be balanced by finding the median diameter of all cell clusters. Cells from each cluster can then be sampled to be equal to this diameter. The method then calculates the variance of each gene across cells in the balanced single-cell UMI matrix. Given the well-known nature of overdispersion in RNA-seq data, calculating the variance directly from the count matrix can be prone to errors. Therefore, the method calculates the variance over normalized data by a variance stabilizing transformation (VST). See Anders, S. & Huber, W., "Differential expression analysis for sequence count data", Genome Biol. (2010). The top 2,000 genes with large variance can be selected. The algorithm for selecting highly variable genes is the same as that programmed in the "Seurat" R package disclosed in A. Butler, Nat. Biotechnol.

[0078] Figure 12B shows a hypothetical example illustrating the cell types that can be selected. Four graphs are shown. There is one graph for each of the four genes, and each graph reflects gene expression versus cell type. In the illustrated example, gene 3 shows good variation across cell types, and gene 4 is a strong marker for cell type C1. Any of these genes have good information about various cell types and can therefore be selected for modeling.

[0079] Step 2: Estimate gene mean & variance per cell type

[0080] Modeling single-cell RNA-seq data can be challenging due to cellular heterogeneity and technical sensitivity and noise. The expression of some genes may not be detected by chance, while other genes may be found to be highly distributed. Distributed genes may lead to excessive variability even within the same cell type. In addition, with an ever-increasing number of cells per study, it is computationally challenging to directly estimate cell percentages using all cells as training samples. The disclosed method combats high noise and computational complexity by collecting individual cells at the cell type level. The mean and variance of each gene per cell type can be estimated. This approach reduces data complexity while retaining cell type-specific information.

[0081] While typical analysis of RNA-seq data begins with normalization, the disclosed method does not normalize before estimating the mean. By normalizing across all cell types, all cell types are forced to have the same amount of RNA transcripts as measured by total UMI counts per cell. However, different cell types may have dramatically different amounts of transcripts. For example, the amount of RNA transcripts in neuronal cells is about 10 times that in glial cells. Thus, normalization may erroneously alter the relative abundance of cell types, misleading estimates of cell type percentages. To avoid this problem, the disclosed method uses raw UMI counts to model the mean.

[0082] Studies have shown that UMI counts follow a negative binomial distribution. See Hafemeister, C. & Satija, R., "Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression," Genome Biol. (2019); and Svensson, V., "Droplet scRNA-seq is not zero-inflated," Nature Biotechnology (2020). Thus, the disclosed method fits a negative binomial distribution to single cells of each cell type. Figure 13 illustrates the steps of estimating mean and variance parameters by fitting a negative binomial distribution to each gene in cell type k. The disclosed method later builds a model based on the estimated mean and variance parameters from selected genes. More specifically, X ik of, Let x be the set of single-cell UMI counts of genes i ∈ 1, ..., I for all cells in cell types k ∈ 1, ..., K. The letter I denotes the number of selected genes, and K denotes the number of cell types on a single-cell basis. ik The distribution of follows a negative binomial distribution.

[0083]

number

[0084]

number

[0085]

number

[0086]

number

[0087] Step 3: Calculating the cross-cell type specificity of genes Genes with cell type specific expression patterns better represent a particular cell type and are therefore more important when used to solve cell type composition. In line with this property, the disclosed method weights genes with high specificity over genes with low specificity. Highly specific genes usually have consistently high expression and relatively low variance among cells within a cell type. To calculate the cell type specificity of a gene, the disclosed method first identifies the cell type in which the gene has the highest expression (i.e., the cell type in which it is most specifically expressed), and then defines the specificity of this gene as the mean-to-variance ratio within the cell type. In the model described below, we assign a high weight to a gene when the ratio is high. The disclosed method uses the estimated mean and variance parameters from the negative binomial fitting from step 2 (μ in the previous equation) and the estimated mean and variance parameters from the negative binomial fitting from step 3 (μ in the previous equation). ik and

number

number

number

[0088]

number

[0089] Step 4: Estimating cross-sample genetic variability Gene variability indicates how stable a gene is across samples. The idea of ​​weighting genes based on their variability across samples was presented in a paper by Wang et al. (supra). Wang et al. defined variability as cross-sample variance. By weighting genes with high variability, the authors realized a significant advantage over traditional unweighted methods. Genes with low cross-sample variability are better representative of the population and therefore more reliable when used to learn cell composition. The disclosed method adopts a similar idea to weight the importance of genes. The disclosed method defines variability in a different and better way. Specifically, the disclosed method defines cross-sample gene variability using the variance-to-mean ratio (VMR), where the mean and variance are calculated across samples. VMR scales better than simple variance and can avoid both under-weighting of lowly expressed genes and over-weighting of unstable genes.

[0090] In addition, the disclosed method can be extended to address applications where multiple samples are not available. There are three options available for calculating VMR, depending on whether multiple sample data is available or not. Typically, the composite transcriptome data to be deconvolved has multiple samples. In bulk RNA-seq data, multiple samples are usually included to control for biological variability. In spatial transcriptome data, adjacent spatial dots can be seen as multiple samples. Thus, in the first option, the disclosed method calculates cross-sample genetic variability from the composite transcriptome data. In the second option, the composite data does not have multiple samples, but the single cell data does, and the disclosed method synthesizes multiple composite samples, each of which is the average of all cells belonging to one of the samples in the single cell criterion. Finally, if multiple samples are not available for both types of data, in the third option, the disclosed method repeatedly bootstraps single cells and averages the sampled cells to create multiple, synthesized composite samples.

[0091] After taking multiple composite samples, Y ij Let us denote the count of sequences for gene i in sample j ∈ 1,...,J by

number

number

[0092]

number

number

[0093]

number

number

[0094] Step 5: Estimating scaling coefficients for each gene When linking composite data to single-cell data, rescaling factors are often used to account for library size and platform differences. All existing methods adopt a single rescaling factor for each unit of sample; that is, all genes in a single sample are multiplied by the same factor (see Wang et al., supra, and Andersson et al., “Spatial mapping of cell types by integration of transcriptomics data,” bioRxiv (2019)). This operation is based on the assumption that the impact of platform differences on all genes is the same and scales linearly across different cell types, which is rarely correct. In addition, estimates can easily be affected by outliers in linear models, so estimates of cell proportions can be pushed away from the truth by extremely high outlier genes. Therefore, applying a uniform scaling factor to all genes is inappropriate.

[0095] To overcome this problem, the disclosed method instead estimates a per-gene scaling factor via an adaptive learning strategy and rescales each gene by its corresponding scaling factor. To proceed, the disclosed method first calculates the average gene expression from the composite sample (from step 4 above)

number

number

[0096]

number

[0097]

number

[0098] Figure 14 shows a hypothetical example that demonstrates the effect of per-gene scaling factors. Ideally, an accurate estimate of the slope (i.e., cell percentage) would be the slope of the leftmost line in Figure 14. However, direct fitting would result in the rightmost line, influenced by outlier genes. Outlier genes may be induced due to different platforms affecting genes differently. The disclosed method employs an adaptive learning approach, where a coarse estimate of the slope (i.e., the rightmost line) is first learned, and then outlier genes are moved towards it, with more outlying genes moving more towards the true line (i.e., along the longer arrows). After this adjustment, the new estimated slope (center line) is closer to the true line (leftmost line) and is therefore a more accurate estimate.

[0099] Step 6: Building a weighted and regularized regression model In the sixth step, the disclosed method builds a model that incorporates all the aforementioned coefficients to provide the actual estimate of cell percentage. The method is built on a non-negative least squares regression model, giving higher weight to genes with high cell type specificity and low cross-sample variability. This step is performed by optimizing a weighted sum of squared loss function L. The weights consist of two components:

number

number

[0100] In the case of complex tissues (e.g., neural tissues) where many highly similar subtypes are common, closely related subtypes may have strong collinearity, leading to overestimating some cell types while underestimating or missing others. The disclosed method addresses this issue by including the L2 norm of the estimate as a regularization component. k Denote by σ the unscaled coefficient for cell type k. For a composite transcriptome sample j, the loss function is given as:

[0101]

number

[0102]

number

[0103]

number

[0104]

number

[0105] In some embodiments, the plurality of genes from the normalized matrix of count-based sequencing data comprises at least about 20,000 genes. In some embodiments, the selected subset of the most variably expressed genes comprises about 1,000 to about 5,000 genes. In some embodiments, the selected subset of the most variably expressed genes comprises about 2500 genes.

[0106] In some embodiments, any of the methods described herein can further include determining a proportion of RNA from each cell type for which bulk or spatial RNA sequencing data was obtained. In some embodiments, any of the methods described herein can further include determining a proportion of each cell type for which bulk or spatial RNA sequencing data was obtained. In some embodiments, any of the methods described herein can further include determining a proportion of RNA from each cell type for which bulk or spatial RNA sequencing data was obtained. In some embodiments, any of the methods described herein can further include determining a proportion of each cell type for which bulk or spatial RNA sequencing data was obtained.

[0107] The method of deconvolving bulk or spatial RNA sequencing data with information from count-based sequencing data (e.g., scRNA-seq data) can be used in a variety of ways. In general, the methods described herein provide more robust and accurate estimates of specific cell types within a population of multiple cell types. In addition, the methods described herein can be applied to all count-based sequencing data (i.e., the methods described herein are not limited to scRNA-seq data, but can be applied to other types of count-based sequencing data, e.g., ATAC-seq, to ​​cellular products other than RNA, and to a wide range of mixed samples, e.g., mixtures of different tissues).

[0108] The methods described herein can be used, for example, to estimate the mixing ratio for one or more specific cell types given the gene expression pattern of a single cell type. Bulk tissues usually consist of multiple cell types in different ratios. Using the liver as an example, there are hepatocytes, stellate cells, Kupffer cells, and endothelial cells. The methods described herein can be used to estimate the ratio of these individual cell types in bulk liver tissue. Mixture ratios for one or more specific cell types can be determined, for example, for organs, tissues, cell cultures, etc.

[0109] The methods described herein can also be used, for example, to detect tissue contamination. For example, a biopsy or other tissue sample obtained from a human may have a first cell type that is desired in the biopsy, but may have or be suspected of having a second, undesired cell type. The methods described herein can be used to determine whether the biopsy or tissue sample is contaminated with a second cell type, and if so, the amount of contamination. To illustrate, muscle contamination is often found in RNA-seq data from cardiac tissue. The methods described herein can be used to determine whether cardiac tissue is contaminated with myocytes during dissection and isolation, and to estimate how many myocytes are present in the cardiac tissue sample.

[0110] The methods described herein can also be used, for example, to detect tumor infiltration. For example, using the methods described herein, a biopsy or other tissue sample can be obtained from a particular tumor within a human, and the presence, identity, and / or extent of infiltration of the tumor by non-tumor cells can be determined. In some embodiments, the non-tumor cells infiltrating the tumor are immune cells, such as macrophages, lymphocytes, and natural killer cells. In some embodiments, the lymphocytes are B lymphocytes and / or T lymphocytes. In some embodiments, the lymphocytes are tumor-infiltrating lymphocytes (TILs). A tumor infiltrated by immune cells (e.g., T lymphocytes) is an "inflammatory" tumor. A tumor that contains only a small number or no infiltrating T lymphocytes (and thus is not recognized and does not provoke a strong response by the immune system) is a "non-inflammatory" tumor. Thus, the methods disclosed herein can be used to determine whether a particular tumor within a human is considered, for example, by a medical professional, to be an inflammatory tumor or a non-inflammatory tumor by determining the presence, identity, and / or extent of infiltration of the tumor by immune cells. Inflammatory tumors are more susceptible to immunotherapy than non-inflammatory tumors. Thus, humans with inflammatory tumors are better candidates for immunotherapy than humans with non-inflammatory tumors. Inflammatory and non-inflammatory tumors are described, for example, in Galon et al., Nat. Rev. Drug Disc., 2019, 18, 197-218; Bonaventura et al., Front. Immunol., 2019, 10, 168, 1-10; and Seidel et al., Front. Oncol., 2018, 8, 86. Thus, the methods described herein can be used to stratify patients for their susceptibility to immunotherapy. Methods can also be used to estimate the proportion of infiltrating cells (e.g., immune cells) within a particular tumor - to identify patients that are susceptible to immunotherapy. In some embodiments, the methods described herein further comprise administering immunotherapy to a human with an invasive tumor.

[0111] In some embodiments, the cells from which the bulk RNA sequencing data was obtained include tumor cells, and the method further includes identifying a proportion of immune cells among the tumor cells. In some embodiments, the immune cells include tumor infiltrating lymphocytes. In some embodiments, the immune cells include CD8 positive T lymphocytes. In some embodiments, the immune cells include CD8 positive T lymphocytes and dendritic cells. In some embodiments, the methods described herein further include characterizing the tumor from which the tumor cells were obtained as an inflammatory tumor or a non-inflammatory tumor.

[0112] In some embodiments, characterizing the tumor as an inflammatory tumor, the tumor is present in a subject, and the method further comprises determining whether the subject has below, equal to, or above a threshold level of infiltrating immune cells. In some embodiments, the immune cells comprise CD8 positive T lymphocytes. In some embodiments, the immune cells comprise CD8 positive T lymphocytes and dendritic cells. In some embodiments, the subject has above a threshold level of infiltrating immune cells, and the method further comprises identifying the subject as a candidate for immunotherapy.

[0113] In some embodiments, the immunotherapy comprises adoptive cell therapy. In some embodiments, the adoptive cell therapy comprises chimeric antigen receptor T cell (CAR-T cell) therapy. In some embodiments, the immunotherapy comprises immune checkpoint inhibition therapy. In some embodiments, the immune checkpoint inhibition therapy comprises an antibody that inhibits cytotoxic T lymphocyte-associated antigen-4 (CTLA-4), an antibody that inhibits programmed cell death protein 1 (PD-1), an antibody that inhibits programmed cell death ligand 1 (PD-L1), or an antibody that inhibits lymphocyte-associated gene 3 (LAG3), or any combination thereof. In some embodiments, the immune checkpoint inhibition therapy comprises an antibody that inhibits cytotoxic T lymphocyte-associated antigen-4 (CTLA-4) (e.g., without limitation, ipilimumab and REGN4659). In some embodiments, the immune checkpoint inhibition therapy comprises an antibody that inhibits programmed cell death protein 1 (PD-1) (e.g., without limitation, nivolumab, pembrolizumab, and cemiplimab). In some embodiments, the immune checkpoint inhibitor therapy includes an antibody that inhibits programmed death-ligand 1 (PD-L1) (e.g., without limitation, atezolizumab). In some embodiments, the immune checkpoint inhibitor therapy includes an antibody that inhibits lymphocyte-associated gene 3 (LAG3) (e.g., without limitation, REGN3767).

[0114] The methods described herein can also be used, for example, to characterize or score the tumor microenvironment. For example, using the methods described herein, a biopsy or other tissue sample can be obtained from a particular tumor within a human, and the presence, identity, and / or extent of infiltration of the tumor microenvironment by tumor microenvironment cells can be determined. Tumor microenvironment cells include, but are not limited to, stromal cells and immune cells. Stromal cells include, but are not limited to, fibroblasts (e.g., cancer-associated fibroblasts), cancer-associated adipocytes, pericytes, and endothelial cells (e.g., lymphatic endothelial cells and vascular endothelial cells). Immune cells include, but are not limited to, macrophages, lymphocytes, and natural killer cells. In some embodiments, the lymphocytes are B lymphocytes and / or T lymphocytes. In some embodiments, the T lymphocytes are TILs.

[0115] When a human has a tumor microenvironment infiltrated by such tumor microenvironment cells, the cancer is at a more advanced stage than when the tumor microenvironment is not infiltrated by such tumor microenvironment cells. Thus, the methods disclosed herein can be used to determine whether a tumor microenvironment within a human is considered, for example by a medical professional, to be at an advanced stage of cancer by determining the presence, identity, and / or extent of infiltration of the tumor microenvironment by tumor microenvironment cells. Thus, the methods described herein can be used to stratify patients for their sensitivity to immunotherapy based on the cell type composition of the tumor microenvironment. The methods can also be used to estimate the proportion of infiltrating cells within a particular tumor microenvironment - which can be used to identify patients that are sensitive to immunotherapy. In some embodiments, the methods described herein further include administering immunotherapy to a human having an infiltrated tumor microenvironment.

[0116] In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of tumor cells among the tumor microenvironment cells. In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of immune cells among the tumor microenvironment cells. In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of cancer-associated fibroblasts among the tumor microenvironment cells. In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of cancer-associated adipocytes among the tumor microenvironment cells. In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of lymphatic endothelial cells among the tumor microenvironment cells. In some embodiments, the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises determining a proportion of vascular endothelial cells among the tumor microenvironment cells.

[0117] The methods described herein can also be used to estimate cell type ratios in samples of, for example, islets of Langerhans, which are clusters of endocrine cells inside the pancreas. Pancreatic islets contain five endocrine cell types (α, β, δ, ε, and γ), of which β cells secrete insulin and are gradually lost in humans with type 2 diabetes. The "normal" population of β cells should be approximately 50-60%. Thus, determining the cell type ratios of pancreatic islet cells by the methods described herein can be used to determine the presence of, and track the development and / or treatment of, type 2 diabetes.

[0118] The methods described herein can also be used, for example, to estimate cell type ratios in a sample of kidney cells to detect kidney diseases characterized by a gradual loss of kidney function, such as chronic kidney disease (CKD). Fibrosis is a histological feature common to all CKD models. In addition to neutrophils and podocytes, kidney cells are classified into two major groups: immune cell types (macrophages, fibroblasts, T lymphocytes, B lymphocytes, and natural killer cells) and kidney-specific cell types (proximal tubules (PT), distal tubules, Henle's loops, two cell types that form the collecting duct, and endothelial cells). Of these, PTs are the predominant cell type in the kidney, and the ratio of PT cells has been found to decrease with CKD progression. Distal tubule cells (DCT) have been found to be the second most abundant cell type in the kidney, with an expected ratio of approximately 10-20%. The ratio of DCT cells shows a consistent increase with CKD disease progression. In addition, immune cells (particularly macrophages, but also fibroblasts, B lymphocytes, and T lymphocytes) have been shown to play a central role in the pathogenesis of CKD, consistent with clinical and histological observations showing that tissue inflammation is a consistent feature of renal fibrosis. Thus, determining the cell type ratio of kidney cells by the methods described herein can be used to determine the presence of kidney disease (e.g., CDK) and to track its development and / or treatment.

[0119] The methods described herein can also be used, for example, to detect the presence and degree of activated or differentiated cells in a cell population. For example, within any cell population, a certain percentage of cells may be activated, and another percentage of cells may be inactive. Similarly, within any cell population, a certain percentage of cells may be differentiated, and another percentage of cells may be undifferentiated. Such cell stages (activated vs. inactive and / or differentiated vs. undifferentiated) can be used in some cases to track the development of a particular cell, for example, from a precursor cell to a mature differentiated cell, or to track the change from a normal cell to a diseased cell.

[0120] In some embodiments, the methods described herein are computer-implemented. The methods may be implemented in software, hardware, firmware, or any combination thereof. In some embodiments, the methods are implemented in one or more computer programs executed on a programmable computer system including at least one processor, a processor-readable storage medium (including, for example, volatile and non-volatile memory and / or storage elements, etc.), and input and output devices. The computer system may include one or more physical machines or virtual machines executed on one or more physical machines. In addition, the computer system may include a computer or a cluster of many distributed computers connected by the Internet or other networks.

[0121] Each computer program may be a set of instructions or program code in a code module resident in a random access memory of a computer system. The set of instructions may be stored in another computer memory (e.g., on a hard disk drive or in a removable memory, such as an optical disk, external hard drive, memory card, or flash drive) until required by the computer system, or may be stored on another computer system and downloaded via the Internet or other network. Each computer program may be implemented in a variety of computer programming languages, including, as an example, Python.

[0122] The disclosed methods (including computer-implemented methods), computer programs, computer systems, and devices for deconvolving bulk RNA sequencing data each recite a wealth of steps and elements that, taken as a whole, go far beyond abstract ideas. As a first matter, the methods, programs, systems, and devices each teach a specific rule-based approach to automating the task of deconvolving bulk RNA sequencing data. The methods, programs, systems, and devices each teach an ordered combination with specific requirements dictated by the individual steps and elements. The specific disclosed steps and elements of these rules are not widespread, and their combination is not a well-understood routine conventional activity. Rather, the specific disclosed steps and elements of these rules enable the improvements realized by the disclosed methods, programs, systems, and devices.

[0123] Moreover, one focus of the disclosed methods, programs, systems, and devices is on certain alleged improvements in computer capabilities. They improve the function of the computer itself. Improvements to computers related to this disclosure include software improvements to logical structures and processes. Most of the advances made in computer technology consist of improvements to software that by their very nature can be defined by logical structures and methods rather than specific physical features. The specific steps and elements of the disclosed methods, programs, systems, and devices constitute a particular type of data structure designed to improve how computers store and retrieve data in memory. The disclosed methods, programs, systems, and devices are directed to improving the functionality of computers and improving the technical task of deconvolving bulk RNA sequencing data. It is the incorporation of the disclosed steps and elements that improves the existing technical task, not the use of computers. Improvements to computer-related technology are not limited to improvements in the operation of computers or computer networks per se, but can also include sets of "rules" (essentially mathematical relationships) that improve computer-related technology.

[0124] Furthermore, the disclosed methods, programs, systems, and devices allow computing devices to do things they could not do before, such as deconvolving bulk RNA sequencing data with higher accuracy and detecting cell types with a proportion of less than about 0.5%. The disclosed methods, programs, systems, and devices provide a necessarily computationally based solution to overcome problems that arise specifically in the area of ​​deconvolving bulk RNA sequencing data. As described herein, the disclosed methods, programs, systems, and devices teach a particular approach to address the computational limitations of existing methods, programs, systems, and devices used to deconvolve bulk RNA sequencing data. The disclosed methods, programs, systems, and devices overcome the shortcomings of existing methods, programs, systems, and devices by deconvolving bulk RNA sequencing data in a particular, new, and non-obvious way, at least accurately and efficiently.

[0125] The present disclosure also provides a computer-readable medium storing processor-executable instructions adapted to cause one or more computing devices to deconvolve bulk RNA sequencing data. In some embodiments, the computer-readable medium storing processor-executable instructions causes one or more computing devices to deconvolve bulk RNA sequencing data by: i) selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) calculating cell-type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using cell type annotations; and iii) for each cell type, for each subset of most variably expressed genes from the count-based sequencing data matrix, the subset of most variably expressed genes. , and cell type annotations to define a mixed single-cell distribution with proportion parameters; iv) fitting bulk distributions to each subset of most variably expressed genes and subsets of most variably expressed genes from a normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; v) defining a loss function between the bulk distribution and the mixed single-cell distribution; and vi) applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. Each of these embodiments can be performed using any of the method embodiments disclosed herein.

[0126] In some embodiments, a computer readable medium storing processor executable instructions can be configured to cause one or more computing devices to perform a deconvolution of bulk RNA sequencing data by: i) taking inputs from three sources (bulk or spatial RNA-seq data, single cell RNA-seq data, and cell type annotations) and selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; and ii) calculating mean and variance parameters of the data per gene per cell type. The method is adapted to perform any one or more of the following steps: iii) estimating the cross-cell type specificity of genes, iv) estimating cross-sample gene variability from the composite data or single cell samples depending on the availability of multiple samples, v) estimating scaling coefficients per gene using both the composite and single cell data, and vi) building a weighted and regularized regression model using all known quantities and using said model to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. Each of these embodiments can be performed using any of the method embodiments disclosed herein.

[0127] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, third, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fifth step. In some embodiments, the method includes a fifth step and one or more of the first, second, third, fourth, and sixth steps, or any combination of these additional steps. In some embodiments, the method includes a sixth step. In some embodiments, the method includes a sixth step and one or more of the first, second, third, fourth, and fifth steps, or any combination of these additional steps.

[0128] In some embodiments, a computer-readable medium storing processor-executable instructions is adapted to cause one or more computing devices to deconvolve the bulk RNA sequencing data by: i) calculating cross-cell type specificity of genes; ii) estimating cross-sample gene variability from the composite data or single-cell samples, depending on the availability of multiple samples; iii) estimating scaling coefficients for each gene using both the composite data and the single-cell data; and iv) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk RNA sequencing data, and thereby inferring the percentage of cell types in the bulk RNA sequencing data. In some embodiments, the calculation of the cross-cell type specificity of genes is based on estimated mean and variance parameters of expression per gene per cell type from a subset of the most variably expressed genes selected from a matrix of count-based single-cell sequencing data (obtained from three sources: i) bulk or spatial RNA-seq data, ii) single-cell RNA-seq data, and iii) cell type annotations, where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells. Each of these embodiments can be performed using any of the method embodiments disclosed herein.

[0129] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, and fourth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, and fourth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, and fourth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, and third step, or any combination of these further steps.

[0130] The disclosure also provides a system that includes one or more processors and a memory having processor-executable instructions that, when executed by the one or more processors, cause the apparatus to deconvolve bulk RNA sequencing data according to any of the methods described herein. In some embodiments, the method includes the steps of: i) selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) calculating cell type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using cell type annotations; and iii) performing cross-sequencing for each cell type and for each subset of most variably expressed genes, the subset of most variably expressed genes, and the cell type annotations from the count-based sequencing data matrix. The method includes any one or more of the following steps: fitting a sample distribution to each subset of the most variably expressed genes and the subset of most variably expressed genes from the normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; v) defining a loss function between the bulk distribution and the mixed single-cell distribution; and vi) applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data.

[0131] In some embodiments, the method comprises any one or more of the following steps: i) taking input from three sources (bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations) and selecting a subset of the most variably expressed genes from a normalized matrix of count-based sequencing data, where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; ii) estimating mean and variance parameters of the data per gene per cell type; iii) calculating cross-cell-type specificity of genes; iv) estimating cross-sample gene variability from the composite data or single-cell samples depending on the availability of multiple samples; v) estimating scaling factors per gene using both the composite and single-cell data; and vi) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. The deconvolution method may be performed using all or a subset of the method embodiments disclosed herein, with or without using the exact methods in each embodiment described herein.

[0132] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, fourth, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, third, fifth, and sixth step, or any combination of these further steps. In some embodiments, the method includes a fifth step. In some embodiments, the method includes a fifth step and one or more of the first, second, third, fourth, and sixth steps, or any combination of these additional steps. In some embodiments, the method includes a sixth step. In some embodiments, the method includes a sixth step and one or more of the first, second, third, fourth, and fifth steps, or any combination of these additional steps.

[0133] In some embodiments, the method includes any one or more of the following steps: i) calculating the cross-cell type specificity of genes; ii) estimating cross-sample gene variability from the composite data or single-cell samples depending on the availability of multiple samples; iii) estimating a scaling factor per gene using both the composite data and the single-cell data; and iv) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data. In some embodiments, the calculation of the cross-cell type specificity of genes is based on estimated mean and variance parameters of expression per gene per cell type from a subset of the most variably expressed genes selected from a matrix of count-based single-cell sequencing data (obtained from three sources: i) bulk or spatial RNA-seq data, ii) single-cell RNA-seq data, and iii) cell type annotations), where the matrix of count-based sequencing data includes count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells. The deconvolution method may be performed using all or a subset of the method embodiments disclosed herein, with or without using the exact methods in each embodiment described herein.

[0134] In some embodiments, the method includes a first step. In some embodiments, the method includes a first step and one or more of a second, third, and fourth step, or any combination of these further steps. In some embodiments, the method includes a second step. In some embodiments, the method includes a second step and one or more of a first, third, and fourth step, or any combination of these further steps. In some embodiments, the method includes a third step. In some embodiments, the method includes a third step and one or more of a first, second, and fourth step, or any combination of these further steps. In some embodiments, the method includes a fourth step. In some embodiments, the method includes a fourth step and one or more of a first, second, and third step, or any combination of these further steps.

[0135] The following exemplary embodiments are given.

[0136] Embodiment 1. A method for deconvolving bulk RNA sequencing data, comprising the steps of: selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; calculating cell-type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using cell type annotations; and, for each cell type, and for each subset of most variably expressed genes from the count-based sequencing data matrix, the subset of most variably expressed genes, and the cell type annotations. fitting a cross sample distribution to each subset of most variably expressed genes and the subset of most variably expressed genes from a normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; defining a loss function between the bulk distribution and the mixed single-cell distribution; and applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby inferring the percentage of cell types in the bulk RNA sequencing data.

[0137] Embodiment 2. The method according to embodiment 1, wherein the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix.

[0138] Embodiment 3. The method according to embodiment 1, wherein the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix.

[0139] Embodiment 4. The method according to any one of embodiments 1 to 3, wherein the cross-sample distribution for each cell type and for each subset of the most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution.

[0140] Embodiment 5. The method according to any one of embodiments 1 to 4, wherein each subset of the most variably expressed genes from the normalized bulk matrix and the bulk distribution for the subset of most variably expressed genes is a bulk Gaussian distribution.

[0141] Embodiment 6. The method according to any one of embodiments 1 to 5, further comprising forming a matrix of count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells, and normalizing the matrix.

[0142] Embodiment 7. The method according to any one of embodiments 1 to 6, further comprising forming a bulk matrix of bulk RNA sequencing counts and normalizing the bulk matrix.

[0143] Embodiment 8. The method according to any one of embodiments 1-7, wherein the selecting step comprises calculating the standard deviation for each gene in the plurality of genes, determining a threshold number of standard deviations, and selecting the subset of most variably expressed genes having standard deviations above the threshold number.

[0144] Embodiment 9. A method according to any one of embodiments 1 to 8, wherein the step of calculating the cell type specific weights comprises comparing the overall average variance to the subcellular type average variance for each of a fixed number of cells.

[0145] Embodiment 10. A method according to any one of embodiments 1 to 9, wherein the fitting step includes using a full distribution when estimating the mixture proportions.

[0146] Embodiment 11. The method according to embodiment 10, wherein the fitting step further comprises obtaining a distribution by fitting the normalized counts to a distribution to estimate the variance and mean for each gene.

[0147] Embodiment 12. The method according to embodiment 11, wherein the distribution is a Gaussian distribution.

[0148] Embodiment 13. A method according to any one of embodiments 1 to 12, wherein the step of defining the loss function comprises applying the Kullback-Leibler divergence.

[0149] Embodiment 14. The method according to any one of embodiments 1 to 13, wherein the step of applying the loss function includes employing a gradient descent method.

[0150] Embodiment 15. The method according to any one of embodiments 1 to 14, wherein the plurality of genes from the normalized matrix of count-based sequencing data comprises at least about 20,000 genes.

[0151] Embodiment 16. The method according to any one of embodiments 1 to 15, wherein the selected subset of the most variably expressed genes comprises from about 1,000 to about 5,000 genes.

[0152] Embodiment 17. The method according to embodiment 16, wherein the selected subset of the most variably expressed genes comprises about 2500 genes.

[0153] Embodiment 18. The method according to any one of embodiments 1 to 17, further comprising identifying the proportion of RNA from each cell type for which bulk RNA sequencing data was obtained.

[0154] Embodiment 19. The method according to any one of embodiments 1 to 18, further comprising identifying the proportion of each cell type for which bulk RNA sequencing data was obtained.

[0155] Embodiment 20. The method according to any one of embodiments 1 to 19, wherein the cells from which bulk RNA sequencing data was obtained include tumor cells, and the method further comprises identifying the proportion of immune cells among the tumor cells.

[0156] Embodiment 21 The method according to embodiment 20, wherein the immune cells comprise tumor-infiltrating lymphocytes.

[0157] Embodiment 22. The method according to embodiment 20 or embodiment 21, wherein the immune cells comprise CD8 positive T lymphocytes.

[0158] Embodiment 23. The method according to embodiment 20, wherein the immune cells include CD8 positive T lymphocytes and dendritic cells.

[0159] Embodiment 24. The method according to any one of embodiments 20 to 23, further comprising characterizing the tumor from which the tumor cells were obtained as an inflammatory tumor or a non-inflammatory tumor.

[0160] Embodiment 25. The method according to embodiment 24, wherein the tumor is characterized as an inflammatory tumor, the tumor being present in a subject, and the method further comprises determining whether the subject has below, equal to, or above a threshold level of infiltrating immune cells.

[0161] Embodiment 26 The method according to embodiment 25, wherein the immune cells comprise CD8 positive T lymphocytes.

[0162] Embodiment 27. The method according to embodiment 25, wherein the immune cells include CD8 positive T lymphocytes and dendritic cells.

[0163] Embodiment 28. The method according to any one of embodiments 25 to 27, wherein the subject is above a threshold level of infiltrating immune cells, and the method further comprises identifying the subject as a candidate for immunotherapy.

[0164] Embodiment 29. The method according to embodiment 28, wherein the immunotherapy comprises adoptive cell therapy.

[0165] Embodiment 30. The method according to embodiment 29, wherein the adoptive cell therapy comprises chimeric antigen receptor T cell (CAR-T cell) therapy.

[0166] Embodiment 31. The method according to embodiment 28, wherein the immunotherapy comprises immune checkpoint inhibitor therapy.

[0167] Embodiment 32. The method according to embodiment 31, wherein the immune checkpoint inhibition therapy comprises an antibody that inhibits cytotoxic T-lymphocyte-associated antigen-4 (CTLA-4).

[0168] Embodiment 33. The method according to embodiment 31 or embodiment 32, wherein the immune checkpoint inhibition therapy comprises an antibody that inhibits programmed cell death protein 1 (PD-1).

[0169] Embodiment 34. The method according to any one of embodiments 31 to 33, wherein the immune checkpoint inhibition therapy comprises an antibody that inhibits programmed death-ligand 1 (PD-L1).

[0170] Embodiment 35. The method according to any one of embodiments 31 to 34, wherein the immune checkpoint inhibition therapy comprises an antibody that inhibits lymphocyte-associated gene 3 (LAG3).

[0171] Embodiment 36. The method according to any one of embodiments 1 to 19, wherein the cells from which the bulk RNA sequencing data was obtained comprise tumor microenvironment cells, and the method further comprises identifying a proportion of tumor cells among the tumor microenvironment cells.

[0172] Embodiment 37. The method according to any one of embodiments 1 to 19, wherein the cells from which bulk RNA sequencing data was obtained include tumor microenvironment cells, and the method further comprises identifying the proportion of immune cells among the tumor microenvironment cells.

[0173] Embodiment 38. The method according to any one of embodiments 1 to 19, wherein the cells from which bulk RNA sequencing data was obtained include tumor microenvironment cells, and the method further comprises identifying a proportion of cancer-associated fibroblasts among the tumor microenvironment cells.

[0174] Embodiment 39. A method according to any one of embodiments 1 to 19, wherein the cells from which bulk RNA sequencing data was obtained include tumor microenvironment cells, and the method further comprises identifying a proportion of cancer-associated adipocytes among the tumor microenvironment cells.

[0175] Embodiment 40. The method according to any one of embodiments 1 to 19, wherein the cells from which the bulk RNA sequencing data was obtained include tumor microenvironment cells, and the method further comprises identifying a proportion of lymphatic endothelial cells among the tumor microenvironment cells.

[0176] Embodiment 41. The method according to any one of embodiments 1 to 19, wherein the cells from which bulk RNA sequencing data was obtained include tumor microenvironment cells, and the method further comprises identifying a proportion of vascular endothelial cells among the tumor microenvironment cells.

[0177]

[0046] Embodiment 42. A computer readable medium storing processor executable instructions, the processor executable instructions being for one or more computing devices to perform a deconvolution of bulk RNA sequencing data comprising the steps of: selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; calculating cell type specific weights for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using cell type annotations; and for each cell type, as well as for each subset of most variably expressed genes from the count-based sequencing data matrix, the subset of most variably expressed genes, and cell type annotations to define a mixed single-cell distribution with proportion parameters; fitting a bulk distribution to each subset of most variably expressed genes and the subset of most variably expressed genes from a normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; defining a loss function between the bulk distribution and the mixed single-cell distribution; and applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby resulting in an inference of cell type percentages in the bulk RNA sequencing data.

[0178] Embodiment 43. The computer-readable medium according to embodiment 42, wherein the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix.

[0179] Embodiment 44. The computer-readable medium according to embodiment 42, wherein the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix.

[0180] Embodiment 45. The computer-readable medium according to any one of embodiments 42 to 44, wherein the cross-sample distribution for each cell type and for each subset of the most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution.

[0181] Embodiment 46. The computer-readable medium according to any one of embodiments 42 to 45, wherein each subset of the most variably expressed genes from the normalized bulk matrix and the bulk distribution for the subset of the most variably expressed genes is a bulk Gaussian distribution.

[0182] Embodiment 47. The computer-readable medium according to any one of embodiments 42-46, wherein the selecting step includes calculating a standard deviation for each gene in the plurality of genes, determining a threshold number of standard deviations, and selecting a subset of the most variably expressed genes having standard deviations above the threshold number.

[0183] Embodiment 48. A computer-readable medium according to any one of embodiments 42 to 47, wherein the step of calculating the cell-type specific weights includes comparing the overall average variance to the within-cell-type average variance for each of a fixed number of cells.

[0184] Embodiment 49. A computer-readable medium according to any one of embodiments 42 to 48, wherein the fitting step includes using a full distribution when estimating the mixture proportions.

[0185] Embodiment 50. The computer-readable medium according to embodiment 49, wherein the fitting step further comprises obtaining a distribution by fitting the normalized counts to the distribution to estimate the variance and mean for each gene.

[0186] Embodiment 51. The computer-readable medium according to embodiment 50, wherein the distribution is a Gaussian distribution.

[0187] Embodiment 52. A computer-readable medium according to any one of embodiments 42 to 51, wherein the step of defining the loss function includes applying Kullback-Leibler divergence.

[0188] Embodiment 53. A computer-readable medium according to any one of embodiments 42 to 52, wherein the step of applying the loss function includes employing a gradient descent method.

[0189] Embodiment 54. The computer-readable medium according to any one of embodiments 42 to 53, wherein the plurality of genes from the normalized matrix of count-based sequencing data comprises at least about 20,000 genes.

[0190] Embodiment 55. The computer-readable medium according to any one of embodiments 42 to 54, wherein the selected subset of the most variably expressed genes comprises from about 1,000 to about 5,000 genes.

[0191] Embodiment 56. The computer-readable medium according to embodiment 55, wherein the selected subset of the most variably expressed genes comprises about 2500 genes.

[0192] Embodiment 57. The computer-readable medium according to any one of embodiments 42 to 56, wherein the method further comprises identifying a proportion of RNA from each cell type for which bulk RNA sequencing data was obtained.

[0193] Embodiment 58. The computer-readable medium according to any one of embodiments 42 to 57, wherein the method further comprises identifying a proportion of each cell type for which bulk RNA sequencing data was obtained.

[0194]

[0036] Embodiment 59. A system comprising one or more processors and a memory having processor-executable instructions, which, when executed by the one or more processors, cause the apparatus to perform deconvolution of bulk RNA sequencing data, comprising the steps of: selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; calculating a cell type specific weight for each selected gene in the subset of most variably expressed genes in the normalized matrix of count-based sequencing data using the cell type annotation; and, for each cell type, and for each subset of most variably expressed genes from the count-based sequencing data matrix, calculating a cell type specific weight for each selected gene in the subset of most variably expressed genes. and a memory, the system comprising: fitting cross-sample distributions to each subset of most variably expressed genes and cell type annotations to define a mixed single-cell distribution with proportion parameters; fitting bulk distributions to each subset of most variably expressed genes and the subset of most variably expressed genes from a normalized bulk matrix to define a bulk distribution, the bulk matrix including bulk RNA sequencing counts for each gene in the plurality of genes for a fixed number of cells; defining a loss function between the bulk distribution and the mixed single-cell distribution; and applying the loss function to estimate cell type proportions in the bulk RNA sequencing data, thereby causing the method to infer cell type percentages in the bulk RNA sequencing data.

[0195] Embodiment 60. The system according to embodiment 59, wherein the count-based sequencing data is single-cell RNA sequencing data, the count-based sequencing counts are single-cell RNA sequencing counts, and the count-based sequencing data matrix is ​​a single-cell RNA sequencing data matrix.

[0196] Embodiment 61. A system according to embodiment 59, wherein the count-based sequencing data is ATAC-seq data, the count-based sequencing counts are ATAC-seq counts, and the count-based sequencing data matrix is ​​an ATAC-seq data matrix.

[0197] Embodiment 62. A system according to any one of embodiments 59 to 61, wherein the cross-sample distribution for each cell type and for each subset of the most variably expressed genes from the count-based sequencing data matrix is ​​a cross-sample Gaussian distribution.

[0198] Embodiment 63. A system according to any one of embodiments 59 to 62, wherein each subset of the most variably expressed genes from the normalized bulk matrix and the bulk distribution for the subset of the most variably expressed genes is a bulk Gaussian distribution.

[0199] Embodiment 64. The system according to any one of embodiments 59 to 63, wherein the selecting step includes calculating the standard deviation for each gene in the plurality of genes, determining a threshold number of standard deviations, and selecting the subset of most variably expressed genes having standard deviations above the threshold number.

[0200] Embodiment 65. A system according to any one of embodiments 59 to 64, wherein the step of calculating the cell type-specific weights includes comparing the overall average variance to the within-cell type average variance for each of a fixed number of cells.

[0201] Embodiment 66. A system according to any one of embodiments 59 to 65, wherein the fitting step includes using a full distribution when estimating the mixture ratio.

[0202] Embodiment 67. The system according to embodiment 66, wherein the fitting step further includes obtaining a distribution by fitting the normalized counts to the distribution to estimate the variance and mean for each gene.

[0203] Embodiment 68. The system according to embodiment 67, wherein the distribution is a Gaussian distribution.

[0204] Embodiment 69. A system according to any one of embodiments 59 to 68, wherein the step of defining the loss function includes applying Kullback-Leibler divergence.

[0205] Embodiment 70. A system according to any one of embodiments 59 to 69, wherein the step of applying the loss function includes employing a gradient descent method.

[0206] Embodiment 71. The system according to any one of embodiments 59 to 70, wherein the plurality of genes from the normalized matrix of count-based sequencing data comprises at least about 20,000 genes.

[0207] Embodiment 72. The system according to any one of embodiments 59 to 70, wherein the selected subset of the most variably expressed genes comprises about 1,000 to about 5,000 genes.

[0208] Embodiment 73. The system according to embodiment 72, wherein the selected subset of the most variably expressed genes comprises about 2500 genes.

[0209] Embodiment 74. A system according to any one of embodiments 59 to 73, wherein the method further comprises identifying a proportion of RNA from each cell type for which bulk RNA sequencing data was obtained.

[0210] Embodiment 75. A system according to any one of embodiments 59 to 74, wherein the method further comprises identifying the proportion of each cell type for which bulk RNA sequencing data was obtained.

[0211] Embodiment 76. A method for deconvolution of bulk or spatial RNA sequencing data, comprising any one or more of the following steps: a) taking input from sources including bulk or spatial RNA-seq data, single cell RNA-seq data, and cell type annotations, and selecting a subset of the most variably expressed genes from a normalized matrix of count-based sequencing data; b) estimating mean and variance parameters of expression per gene per cell type; c) calculating cross-cell type specificity for each gene; d) estimating cross-sample gene variability from bulk or spatial RNA-seq data or single cell samples depending on availability of multiple samples; e) estimating scaling factors per gene using both bulk or spatial RNA-seq data and single cell data; and f) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data.

[0212] Embodiment 77. The method according to embodiment 76, wherein the matrix of count-based sequencing data comprises count-based sequencing counts for each gene in the plurality of genes for a fixed number of cells.

[0213] Embodiment 78. A method according to embodiment 76 or embodiment 77, wherein the input is a single-cell UMI count matrix with a cell type annotation associated with each cell.

[0214] Embodiment 79. The method according to embodiment 76 or embodiment 77, wherein the bulk data to be deconvolved is transcripts per kilobase million (TPM) or read counts.

[0215] Embodiment 80. A method according to embodiment 76 or embodiment 77, wherein the spatial data to be deconvolved is a UMI count matrix.

[0216] Embodiment 81. A method according to any one of embodiments 76 to 80, wherein uninformative genes that may introduce noise are filtered out and genes containing important information for distinguishing cell types are selected as a subset of the most variably expressed genes from a normalized matrix of count-based sequencing data.

[0217] Embodiment 82. The method according to embodiment 81, wherein a subset of genes is selected from a population of marker genes whose expression is enriched in each cell type within the single cell UMI count matrix.

[0218] Embodiment 83. A method according to embodiment 82, in which the built-in tool takes as input a single cell UMI count matrix and cell type annotations.

[0219] Embodiment 84. The method according to embodiment 83, wherein for each cell type, the tool calculates the fold change between the average UMI in that cell type and the average UMI in all other cell types and ranks the genes by descending order of fold change.

[0220] Embodiment 85. The method according to embodiment 84, wherein about the top 200 genes are selected from each cell type.

[0221] Embodiment 86. The method according to embodiment 84, wherein a selected marker gene is selected that is present in five or fewer cell types.

[0222] Embodiment 87. The method according to embodiment 84, wherein the selected marker gene is selected to be present in a fixed number of cell types or in a proportion of the total number of cell types, whichever is smaller.

[0223] Embodiment 88. The method according to embodiment 84, wherein a total of about 1,000 unique genes are selected.

[0224] Embodiment 89. The method according to embodiment 81, wherein the subset of genes is selected from a population of highly variable genes that are most variable among all cells in the single cell UMI count matrix.

[0225] Embodiment 90 The method according to embodiment 89, wherein the variance for each gene after cell number balancing and VST normalization is calculated.

[0226] Embodiment 91 The method according to embodiment 90, wherein the gene with the highest variance is selected.

[0227] Embodiment 92. A method according to any one of embodiments 89 to 91, wherein the cell types in the single cell UMI count matrix are balanced by finding the median diameter of all cell clumps and sampling cells from each clump to be equal to this diameter.

[0228] Embodiment 93. A method according to embodiment 92, in which the variance of each gene across cells in the balanced single-cell UMI matrix is ​​calculated.

[0229] Embodiment 94. The method according to embodiment 93, in which the variance for the normalized data is calculated by a variance stabilizing transformation (VST).

[0230] Embodiment 95. The method according to embodiment 94, wherein the top 2,000 genes with large variance are selected.

[0231] Embodiment 96. A method according to any one of embodiments 76 to 95, wherein the RNA-seq data is not normalized before estimating the mean value.

[0232] Embodiment 97. The method according to embodiment 96, wherein the average value is modeled using raw UMI counts.

[0233] Embodiment 98. The method according to embodiment 96 or embodiment 97, wherein a negative binomial distribution is fitted to a single cell of each cell type.

[0234] Embodiment 99. The method according to embodiment 98, in which prediction is made for each selected gene in each cell type.

[0235] Embodiment 100. A method according to any one of embodiments 76 to 99, wherein to calculate the cell type specificity of a gene, the cell type in which the gene i) has the highest expression or ii) has the highest fold change relative to others is identified, and the specificity of this gene is defined as the ratio of the mean to the variance within the cell type.

[0236] Embodiment 101. A method according to embodiment 100, in which the estimated mean and variance parameters from the negative binomial fitting are used to calculate a cell type specificity weight for each gene in the set of selected genes.

[0237] Embodiment 102. The method according to any one of embodiments 76 to 101, wherein the cross-sample genetic variability is calculated using the variance-to-mean ratio (VMR) calculated across the samples.

[0238] Embodiment 103. The method according to embodiment 102, wherein cross-sample genetic variability is calculated from composite transcriptome data.

[0239] Embodiment 104. A method according to embodiment 102, wherein the composite data does not have multiple samples, but the single cell data has multiple samples, and multiple composite samples are synthesized, each of which is the average value of all cells belonging to one of the samples in the single cell criterion.

[0240] Embodiment 105. A method according to embodiment 102, wherein if multiple samples are not available for both the composite data and the single cell data, the method includes generating multiple synthetic samples for the single cell data by averaging expression of a subset of cells.

[0241] Embodiment 106. A method according to any one of embodiments 76 to 105, wherein a scaling factor for each gene is estimated using an adaptive learning strategy and each gene is rescaled by its corresponding scaling factor.

[0242] Embodiment 107. A method according to any one of embodiments 76 to 106, wherein each composite sample is independently estimated by a regression model.

[0243] Embodiment 108. A method for deconvolution of bulk or spatial RNA sequencing data, comprising any one or more of the following steps: a) calculating cross-cell type specificity for each gene within a subset of most variably expressed genes selected from a normalized matrix of count-based sequencing data obtained from sources including bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations; b) estimating cross-sample gene variability from bulk or spatial RNA-seq data or single-cell samples depending on the availability of multiple samples; c) estimating scaling factors per gene using both bulk or spatial RNA-seq data and single-cell data; and d) constructing a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data.

[0244]

[0046] Embodiment 109. A computer-readable medium storing processor-executable instructions, the processor-executable instructions configuring one or more computing devices to perform a deconvolution of bulk or spatial RNA sequencing data by: i) taking input from sources including bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations to select a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; and ii) calculating mean and variance parameters of expression per gene per cell type. The computer-readable medium is adapted to perform a method of inferring cell type percentages in bulk or spatial RNA sequencing data, the method comprising: iii) estimating cross-sample gene variability from the composite data or single cell samples depending on the availability of multiple samples; v) estimating scaling coefficients per gene using both the composite data and the single cell data; and vi) building a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data.

[0245] Embodiment 110. A computer-readable medium storing processor-executable instructions adapted to cause one or more computing devices to deconvolve bulk or spatial RNA sequencing data by a method comprising any one or more of the following steps: i) calculating cross-cell type specificity of genes within a subset of most variably expressed genes selected from a normalized matrix of count-based sequencing data obtained from sources including bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations; ii) estimating cross-sample gene variability from composite data or single-cell samples depending on the availability of multiple samples; iii) estimating scaling coefficients per gene using both the composite data and the single-cell data; and iv) constructing a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data.

[0246]

[0036] Embodiment 111. A system comprising one or more processors and a memory having processor-executable instructions, the processor-executable instructions, when executed by the one or more processors, causing the apparatus to perform deconvolution of bulk or spatial RNA sequencing data by: i) taking inputs from sources including bulk or spatial RNA-seq data, single cell RNA-seq data, and cell type annotations, and selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data, the matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; and ii) selecting a subset of most variably expressed genes from a normalized matrix of count-based sequencing data including count-based sequencing counts for each gene in a plurality of genes for a fixed number of cells; iii) estimating mean and variance parameters of the data from the combined data; iv) estimating cross-sample gene variability from either the combined data or single cell samples depending on the availability of multiple samples; v) estimating scaling coefficients per gene using both the combined data and the single cell data; and vi) constructing a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby resulting in a method of inferring cell type percentages in the bulk or spatial RNA sequencing data.

[0247] Embodiment 112. A system comprising one or more processors and a memory having processor-executable instructions that, when executed by the one or more processors, cause the device to deconvolve bulk or spatial RNA sequencing data in a manner that includes any one or more of the following steps: i) calculating cross-cell type specificity of genes within a subset of most variably expressed genes selected from a normalized matrix of count-based sequencing data obtained from sources including bulk or spatial RNA-seq data, single-cell RNA-seq data, and cell type annotations; ii) estimating cross-sample gene variability from composite data or single-cell samples depending on the availability of multiple samples; iii) estimating scaling coefficients per gene using both the composite data and the single-cell data; and iv) constructing a weighted and regularized regression model using all known quantities and using the model to estimate cell type proportions in the bulk or spatial RNA sequencing data, thereby inferring the percentage of cell types in the bulk or spatial RNA sequencing data.

[0248] In order that the subject matter disclosed herein may be more efficiently understood, the following examples are provided. It should be understood that these examples are merely illustrative and should not be construed as limiting the claimed subject matter in any manner.

[0249] Example 1: Deconvolution of bulk RNA sequencing data for immune cells The following hypothetical examples are included to more clearly demonstrate the general nature of the present disclosure. The examples are exemplary of the disclosure and are not limiting.

[0250] For the single cell matrix, which contains three cell types: macrophages, T cells, and B cells, there are three genes measured by RNA sequencing (Spp1, Trem2, and Serpine2).

[0251] Single cell data input: Raw UMI count data [Table 1] Normalize to 1e+6 total UMI counts per column. [Table 2] Add 1 to every number [Table 3] Calculate the natural logarithm transformation [Table 4] Applying step 1 of the disclosed method, the top 2,500 highly variable genes are selected. In this case, the gene Serpine2 is less variable. The remaining two genes Spp1 and Trem2 are more variable and therefore are included in the top 2,500 genes selected for further analysis.

[0252] In step 2, the cell-type specific weights are calculated. [Table 5] The calculation for the Spp1 gene is as follows:

[0253]

number

[0254] Sample 1: [Table 7] Cells are pooled (summed up) for each cell type to generate: [Table 8] After normalization, log(data + 1) is taken to produce the following. [Table 9] Sample 2: [Table 10] Cells are pooled (summed up) for each cell type to generate: [Table 11] After normalization, log(data + 1) is taken to produce the following. [Table 12] Sample 3: [Table 13] Cells are pooled (summed up) for each cell type to generate: [Table 14] After normalization, log(data + 1) is taken to produce the following. [Table 15] Continuing with step 3, a Gaussian distribution is fitted across each of the three samples for each of the two genes in each of the three cell types. The results of these Gaussian distribution fits are shown in Figure 8. These results conclude the application of step 3.

[0255] Step 4 of the disclosed method is applied to normalize the multi-sample matrix and fit Gaussians to each gene. The results of these Gaussian fits are shown in Figure 9. We assume that bulk RNA-seq data is acquired independently from three samples. The cellular content within each sample is unknown. One of the goals of the disclosed method is to learn cell type ratios for each sample using previous single-cell RNA-seq data as a baseline.

[0256] Applying step 5 of the disclosed method, a model is defined with the proportion parameters, weights, and distributions for each gene learned from the single cell and bulk data. The results of these calculations are shown in FIG.

[0257] Finally, step 6 of the disclosed method is applied to determine the proportions θ of each of the three cell types (k=1, 2, 3) in the bulk tissue. k for each of the three samples, kl The outputs for the first sample are θ1(macrophage)=0.5; θ2(T cell)=0.22; and θ3(B cell)=0.28. The outputs for the second sample are θ1(macrophage)=0.42; θ2(T cell)=0.32; and θ3(B cell)=0.26. The outputs for the third sample are θ1(macrophage)=0.23; θ2(T cell)=0.38; and θ3(B cell)=0.39. Naturally, the sum of the proportions of the three cell types must equal 1.

[0258] Example 2: Method evaluation To evaluate the AdRoit method, two comparisons were made. First, the results of the second embodiment of the AdRoit method were compared to those achieved by the "Multi-Subject Single Cell Deconvolution" (MuSiC) method disclosed in the paper by X. Wang et al. (supra). Second, the results of the second embodiment of the AdRoit method were compared to those achieved by a conventional non-negative least squares (NNLS) regression method.

[0259] Assessment 1: Human islet data The data used for the first evaluation were obtained from human pancreatic islets. Islets of Langerhans are regions of the pancreas that contain its endocrine (i.e., hormone-producing) cells. Single cell data from human pancreatic islets are shown in Figures 15A and 15B. These data were selected for comparison because they contain several (specifically, four) cell types, including two major cell types (alpha and beta cells) and two minor cell types (PP and delta cells), from many (specifically, 18) subjects. Cell percentages vary between different subjects.

[0260] FIG. 15A is an overview of the cellular composition of 18 subjects. FIG. 15B reflects t-distributed stochastic neighbor embedding (t-SNE), a machine learning algorithm for visualization developed by Laurens van der Maaten and Geoffrey Hinton. It is a nonlinear dimensionality reduction technique well suited to embedding high-dimensional data for visualization in a two-dimensional low-dimensional space (e.g., the graph in FIG. 15B). Specifically, t-SNE models each high-dimensional object by a two-dimensional point such that, with a high probability, similar objects are modeled by nearby points and dissimilar objects are modeled by distant points. As is typical of t-SNE graphs, the data is displayed in clumps in FIG. 15B.

[0261] To obtain absolute correct or "true" results to use in evaluating the different methods, bulk data corresponding to the single cell data of human islets shown in Figures 15A and 15B were synthesized. Estimation was performed by using single cell standards without the samples used to synthesize the bulk data. The true bulk data is reflected in Figures 16A, 16B, and 16C. Four separate statistical measures were calculated: mean absolute deviation (mAD), root mean square deviation (RMSD), Pearson correlation coefficient (PCC), and Spearman's rank correlation coefficient.

[0262] In statistics, mAD is the average of the absolute deviations of a set of data about the mean of that data. Mean absolute deviation is also called average deviation. Mean absolute deviation is a way of describing the variation within a data set. The smaller the mAD number, the less variable (i.e., the better) the data set is.

[0263] RMSD or root mean square error (RMSE) (or often root mean square error) is a frequently used measure of the difference between the values ​​predicted by a model or estimator (sample or population values) and the observed values. RMSD represents the square root of the second sample moment of the difference between the predicted and observed values ​​or the mean square of these differences. These deviations are called residuals when calculated on the data sample used for the estimation, and errors (or prediction errors) when calculated outside the sample. RMSD serves to collect the magnitude of errors in predictions for various time points into a single measure of predictive ability. RMSD is scale dependent, so it is a measure of accuracy for comparing prediction errors of different models for a specific data set, not across data sets. For mAD, the smaller the RMSD number, the better.

[0264] In statistics, PCC (also known as Pearson's r, Pearson's product-moment correlation coefficient (PPMCC), or bivariate correlation) is a measure of linear correlation between two variables X and Y. According to the Cauchy-Schwarz inequality, it has values ​​from +1 to -1, where +1 is total positive linear correlation, 0 is no linear correlation, and -1 is total negative linear correlation. PCC is widely used in science. It was developed by Karl Pearson from related ideas introduced by Francis Galton in the 1880s, and a mathematical formula for it was derived and published by Auguste Bravais in 1844.

[0265] In statistics, the Spearman rank correlation coefficient or Spearman's rho (named after Charles Spearman and written with the Greek letter ρ (rho) or r s Rank correlation (often denoted as ) is a nonparametric measure of rank correlation (the statistical dependence between the rankings of two variables). It evaluates how well the relationship between two variables can be described using a monotonic function. The Spearman correlation between two variables is equal to the Pearson correlation between the rank values ​​of these two variables. Pearson correlation evaluates a linear relationship, whereas Spearman correlation evaluates a monotonic relationship (linear or not). For both Pearson correlation and Spearman correlation, the closer the number is to 1, the better.

[0266] 16A, 16B, and 16C reflect a comparison between the results of the AdRoit method (FIG. 16A), the MuSiC method (FIG. 16B), and the NNLS method (FIG. 16C). Therefore, three graphs are shown. The ordinate (Y-axis) for each graph is the estimated proportion of cell types given by the corresponding method. The abscissa (X-axis) for each graph is the true proportion of cell types (from the synthesized bulk data). Four separate statistical measures (summarized above) were calculated for each of the three graphs and are tabulated in FIG. 17. For the mAD and RMSD measures, the smaller the number (i.e., the smaller the deviation), the more accurate the method. The Spearman and Pearson correlations approach 1 for more accurate methods. The data support the conclusion that the AdRoit method is highly accurate and superior to the MuSiC and NNLS methods. AdRoit achieves the highest accuracy when applied to synthesized bulk data using single cell data from human pancreatic islets.

[0267] Assessment 2: Human trabecular meshwork data The data used for the second evaluation was obtained from human trabecular meshwork (TM). The TM is a region of tissue in the eye located around the base of the cornea, near the ciliary body, and is responsible for draining aqueous humor from the eye via the anterior chamber (the cavity at the front of the eye covered by the cornea). Single cell data from human TM are shown in Figures 18A and 18B. These data were selected for comparison because they contain a large number (specifically, 12) cell types from many (specifically, 8) donors. See Patel, G. et al., "Molecular taxonomy of human ocular outflow tissues defined by single-cell transcriptomics," Proc. Natl. Acad. Sci. 117, 12856LP-12867 (2020). The cell types are listed in Figure 18A. Cell percentages vary between different donors.

[0268] Figure 18A summarizes the cellular composition of the eight donors. Figure 18B reflects the t-SNE visualization, and as is typical for t-SNE graphs, the data is displayed in clusters in Figure 18B.

[0269] To obtain absolute correct or "true" results to use in evaluating the various methods, data corresponding to the bulk Human TM single cell data shown in Figures 18A and 18B were synthesized. Estimation was performed by using a single cell standard without the donors used to synthesize the bulk data. The true bulk data is reflected in Figures 19A, 19B, and 19C. Again, four separate statistical measures were calculated: mAD, RMSD, PCC, and Spearman's rank correlation coefficient.

[0270] 19A, 19B, and 19C reflect a comparison between the results of the AdRoit method (FIG. 19A), the MuSiC method (FIG. 19B), and the NNLS method (FIG. 19C). Therefore, three graphs are shown. The ordinate (Y-axis) for each graph is the estimated proportion of the cell type given by the corresponding method. The abscissa (X-axis) for each graph is the true proportion of the cell type (from the combined bulk data). Four separate statistical measures (summarized above) were calculated for each of the three graphs and are tabulated in FIG. 20. Overall, the AdRoit method has the closest estimates to the true cell fractions compared to the MuSiC and NNLS methods. The AdRoit method has the lowest mAD and RMSD, and the highest Pearson and Spearman correlations.

[0271] Figure 21 shows 12 bar graphs, one for each cell type. The dots on the graph represent each of the eight different donors, and the bars represent the 1.5x interquartile range. For each cell type in the TM, the AdRoit method has the smallest deviation from the true cell proportions and is the most accurate estimate across the eight different donors.

[0272] Figure 22 reflects the estimated and true data calculated using both the AdRoit and MuSiC methods. The synthetic bulk data was simulated by using only 6 of the 12 cell types and estimated with reference to the complete list of all 12 cell types. The AdRoit method had fewer false positive estimates for the 6 cell types excluded from the simulation and more accurate estimates for the 6 cell types included in the simulation.

[0273] FIG. 23 is a Receiver Operating Characteristic (ROC) curve showing that the AdRoit method had a significantly higher area under the curve (AUC) than the MuSiC method, indicating better sensitivity and specificity. ROC curves are graphical plots illustrating the diagnostic ability of a binary classification system as its discrimination threshold is varied. The diagnostic method was developed for military radar receiver operators and is therefore so named. ROC curves are constructed by plotting the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings. TPR is also known as sensitivity. FPR is also known as the probability of a false alarm and can be calculated as 1-specificity. ROC analysis provides a tool to select the best possible model and discard the suboptimal ones.

[0274] Evaluation 3: Dorsal root ganglion data The data used for the third evaluation were obtained from mouse dorsal root ganglion (DRG) neurons. DRG single-cell RNA-seq data are shown in Figure 24A and Figure 24B. These data were selected for comparison because they contain many (specifically, 14) cell types (including multiple subtypes of neurons) from five mice. Figure 24A is an overview of the cellular composition of the five mice, listing the cell types. Figure 24B reflects a t-SNE visualization of the data. Cellular proportions vary between different mice.

[0275] To obtain absolute correct or "true" results to use in evaluating the various methods, bulk data corresponding to the DRG single cell data shown in Figures 24A and 24B were synthesized. Estimation was performed by using single cell standards without the samples used to synthesize the bulk data. The true bulk data is reflected in Figures 25A, 25B, and 25C.

[0276] Figures 25A, 25B, and 25C reflect a comparison between the results of the AdRoit method (Figure 25A), the MuSiC method (Figure 25B), and the NNLS method (Figure 25C). Therefore, three graphs are shown. The ordinate (Y axis) for each graph is the estimated proportion of the cell type given by the corresponding method. The abscissa (X axis) for each graph is the true proportion of the cell type (from the combined bulk data). For each individual sample, the mAD, RMSD, and Pearson and Spearman correlations were calculated and compared between the three methods. The results are shown in the graphs in Figure 26. Overall, the AdRoit method has the closest estimates to the true cell fractions compared to the MuSiC and NNLS methods. The AdRoit method has the lowest mAD and RMSD, and the highest Pearson and Spearman correlations. In addition, the AdRoit method estimates were the most stable between samples.

[0277] Evaluation 4: Application of human pancreatic islets The data used for the fourth evaluation were obtained from human pancreatic islets (see Evaluation 1 above). Single cell data from human islets are shown in Figures 15A and 15B. These data include four cell types: alpha, beta, PP, and delta cells. Figure 27 is a graph of cell percentages for each of the four cell types, showing that AdRoit method estimates of cell type percentages for real human islet bulk RNA-seq data are highly reproducible for repeated samples from the same donor.

[0278] Ribonucleic acid molecule-targeted fluorescent in situ hybridization (RNA FISH) is a method for detecting and localizing specific RNA molecules in fixed cells. The detection uses a nucleic acid probe that is complementary to a target RNA sequence in the cell. The probe then hybridizes to its target via standard Watson-Crick base pairing, after which the target can be detected via a fluorescent microscope through direct conjugation of a fluorescent molecule to the probe or through a fluorescent signal amplification method. Recent advances in RNA FISH have increased the specificity and sensitivity of the method, allowing the detection of individual RNA molecules, providing accurate measurements of RNA abundance and localization at the single cell or even subcellular level. While most applications to date have been in fixed cells, advances in probe technology have allowed the detection of single RNA molecules in live cells. Figure 28 shows that cell fraction percentages estimated using the Adroit method are consistent with RNA FISH measurements of cell type percentages.

[0279] Glycated hemoglobin, or HbA1c, is formed when glucose (sugar) in the body attaches to red blood cells. HbA1c testing is used to monitor patients with type 2 diabetes (T2D). In these patients, the body cannot properly use sugar, and sugar tends to attach to blood cells and build up in the blood. Red blood cells are active for about 2-3 months, which is why HbA1c testing is done four times a year. High HbA1c means the patient has too much sugar in their blood and is more likely to develop diabetes complications (e.g., eye and foot problems in patients). In T2D patients, the ideal HbA1c level is 48mmol / mol (6.5%) or less. Figure 29 shows that beta cell fraction percentage estimated using the Adroit method has a significant linear relationship with donor HbA1C levels (including both healthy and T2D cells). FIG. 30 shows that beta cell fraction percentages estimated using the Adroit method are significantly lower in T2D patients than in healthy subjects.

[0280] Rating 5: Simulated spatial spot The data used for the fifth evaluation compares the AdRoit method to stereoscopic estimation. Spatial transcriptomics (ST) is a technique used to spatially analyze RNA-seq data (and therefore all mRNAs) in individual tissue sections. Spatially barcoded reverse transcription oligo(dT) primers are ordered and attached to the surface of a microscope slide in an array of spots, allowing encoding and maintenance of positional information throughout mRNA sample processing and subsequent sequencing. This is in contrast to RNA sequencing of single cells or sequencing of bulk RNA extracted from tissue volumes, where precise spatial information is lost. When tissue cryosections are mounted on spatial transcriptomics slides, barcoded primers bind and capture adjacent mRNAs from the tissue. While the tissue sections are mounted on the slides, reverse transcription of the captured mRNA is initiated and the resulting cDNA incorporates the spatial barcode of the primers. Following mRNA capture and reverse transcription, sequencing libraries are prepared and analyzed by Illumina dye sequencing. The spatial barcode present within each generated sequence allows the data for each individual mRNA transcript to be mapped back to its point of origin within the tissue section.

[0281] Stereoscopic vision (also called stereoscopy or stereo imaging) is a technique for creating or enhancing the illusion of depth in an image by stereoscopic vision as opposed to binocular vision. A stereoscope is a type of image viewer that creates the illusion of a three-dimensional image from two similar two-dimensional images through the use of mirrors or lenses. Complex cellular structures are often rendered in stereoscopic photographs.

[0282] In Fig. 31 we compare the estimates achieved by the stereoscopic and AdRoit methods on a simulated spatial spot containing five different PEP cell subtypes. The true mixture proportions are indicated by the vertical red dashed lines. Three schemes were simulated: (1) scheme 1, left side of Fig. 31, the proportions of the five PEP cell types were the same and equal to 0.2; (2) scheme 2, center of Fig. 31, one PEP cell type was 0.1 and the other four were 0.225; (3) scheme 3, right side of Fig. 31, two PEP cell types were 0.1, two were 0.2, and one was 0.4. For all simulation schemes, the AdRoit estimates converged more consistently to the true simulated proportions than did the stereoscopic estimates.

[0283] Figure 32 illustrates the simulation of a very low percentage of a single type of PEP cell. The percentages were 0.02, 0.04, 0.06, 0.08, and 0.1. The true mixture fraction is indicated by the horizontal dashed red line. The median estimate achieved using the AdRoit method was closer to the true fraction than the estimate achieved using stereo vision. In Figure 33, the estimates using stereo vision and the AdRoit method are compared via a graph of detection rate versus simulated fraction for very low amounts of six different cell types. The AdRoit method was more sensitive in detecting low percentage cells and was consistent across different mixtures of cell types.

[0284] Evaluation 6: Mouse brain spatial transcriptome application example Data for the sixth assessment were obtained from mouse brain cell types. The Allen Mouse Brain Atlas is a genome-wide, high-resolution atlas of gene expression throughout the adult mouse brain. The atlas provides genome-wide in situ hybridization (ISH) imaging data for approximately 20,000 genes in adult mice. Each data set is processed through an informatics analysis pipeline to obtain spatially mapped, quantified expression information. See Lein, E. et al., "Genome-wide atlas of gene expression in the adult mouse brain," Nature 445, 168-176 (2007).

[0285] Figure 34 illustrates how spatial mapping of three cell types using the AdRoit method quantitatively shows the content at each spot. Figure 35 shows ISH images of Wfs1, Prox2, and Rarres2 cell types from the Allen mouse brain atlas. Comparison of Figure 34 and Figure 35 shows that the ISH images in Figure 35 are consistent with the corresponding cell locations in Figure 34.

[0286] Although illustrated and described above with reference to certain specific embodiments and examples, it is nevertheless not intended that the disclosure be limited to the details shown. Rather, various modifications may be made in the details within the scope and range of equivalents of the claims without departing from the spirit of the disclosure.

[0287] In addition to those described herein, various modifications of the described subject matter will be apparent to those skilled in the art from the foregoing description. Such modifications are intended to fall within the scope of the appended claims. Each reference cited in this application (e.g., without limitation, original articles, U.S. and non-U.S. patents, patent application publications, international patent application publications, GenBank accession numbers, etc.) is incorporated herein by reference in its entirety.

Claims

1. A computer-implemented method for deconvolving bulk or spatial RNA sequencing data for a sample, the method comprising: a) obtaining, by a computer, inputs from sources including bulk or spatial RNA-seq data for the sample, single cell RNA-seq data, and cell type annotations; b) generating, by the computer, a single-cell UMI (unique molecular identifier) ​​count matrix based on the bulk or spatial RNA-seq data, the single-cell RNA-seq data, and cell type annotations, where each column of the UMI count matrix corresponds to a cell in the sample or a spatial location in the sample, and each row of the UMI count matrix corresponds to a gene in the sample; c) selecting, by the computer, a subset of the most variably expressed genes from the single cell UMI count matrix based on the population of genes whose expression is enriched in each cell type in the single cell UMI count matrix, or the population of genes whose expression varies above a predetermined level among all cells in the single cell UMI count matrix, or both; d) estimating, by the computer, mean and variance parameters of the genes selected from the single-cell UMI count matrix by fitting a negative binomial distribution to each gene in a cell type; e) calculating, by the computer, a cross-cell type specificity for each of the selected genes based on the estimated mean and variance parameters; f) estimating, by the computer, cross-sample genetic variability from the bulk or spatial RNA-seq data and / or single cell samples based on the ratio of mean gene expression levels and variances of the selected genes; g) estimating, by the computer, scaling coefficients for each gene using both the bulk or spatial RNA-seq data and the single-cell RNA-seq data using a regression-based approach; h) generating, by the computer, a regression model based on the cross-cell type specificity of each selected gene, the cross-sample gene variability from the bulk or spatial RNA-seq data or the single cell sample, and a scaling coefficient for each gene; i) estimating, by the computer, cell type proportions in the bulk or spatial RNA sequencing data using the regression model, wherein the regression model gives higher weight to selected genes with high cell type specificity and low cross-sample variability; Including, As a result, the computer infers the percentage of cell types in the bulk or spatial RNA sequencing data.

2. 2. The method of claim 1, wherein the bulk data to be deconvolved is transcripts per kilobase million (TPM) or read counts. [0010] Claim 3: The subset of most variably expressed genes is selected based on a population of genes whose expression is enriched in each cell type within the single-cell UMI count matrix; The method of claim 1 or 2, wherein an embedded tool takes as input the single cell UMI count matrix and cell type annotations.

4. 4. The method of claim 3, wherein for each cell type, the built-in tool calculates the fold change between the average UMI in that cell type and the average UMI in all other cell types and ranks genes by descending order of fold change.

5. The method of claim 4, wherein about the top 200 genes from each cell type are selected.

6. The method of claim 4, wherein a selection gene is selected that is present in five or fewer cell types.

7. 5. The method of claim 4, wherein a selection gene is selected that is present in a fixed number of cell types or a portion of the total number of cell types, whichever is smaller.

8. The method of claim 4, wherein a total of about 1,000 unique genes are selected.

9. The method of claim 8, wherein the subset of most variably expressed genes is selected based on a population of genes that vary above a predetermined level among all cells in the single-cell UMI count matrix; The method of claim 1 or 2, wherein the variance for each gene after cell number balancing and VST normalization is calculated.

10. The method of claim 9, wherein the gene with the highest variance is selected.

11. The method of claim 9 or 10, wherein the cell types in the single cell UMI count matrix are balanced by finding the median size of all cell clumps, and cells from each clump are sampled to be equal to this size.

12. 12. The method of claim 11, wherein the variance of each gene across the cells in the balanced single-cell UMI matrix is ​​calculated.

13. The method of claim 12, wherein the variance for the single-cell UMI matrix is ​​calculated by a variance stabilizing transformation (VST).

14. The method of claim 13, wherein the top 2,000 genes with large variance are selected.

15. The method of any one of claims 1 to 14, wherein the RNA-seq data is not normalized before estimating the mean value.

16. The method of claim 15 , wherein the average value is modeled using raw UMI counts.

17. 17. The method of claim 15 or 16, wherein predictions are made for each selected gene in each cell type.

18. 2. The method of claim 1, wherein the mean and variance parameters estimated from the negative binomial fitting are used to calculate a cell type specificity weight for each gene in the set of selected genes.

19. 2. The method of claim 1, wherein the cross-sample genetic variability is calculated from composite transcriptome data.

20. 2. The method of claim 1, wherein the bulk or spatial RNA-seq data does not have multiple samples, but the single cell data does have multiple samples, and multiple composite samples are synthesized, each being the average of all cells belonging to one of the samples on a single cell basis.

21. The method of claim 1, wherein if multiple samples are not available for both the bulk or spatial RNA-seq data and the single cell data, the method includes generating, by the computer, multiple synthetic samples for the single cell data by averaging the expression of a subset of cells.

22. The method of any one of claims 1 to 21, wherein the per-gene scaling factor is estimated using an adaptive learning strategy, and each gene is rescaled by its corresponding scaling factor.

23. A method according to any one of claims 1 to 22, wherein composite samples are independently estimated by the regression model.

24. A computer-readable medium storing processor-executable instructions for execution by one or more computing devices, comprising: The processor-executable instructions, when executed by the one or more computing devices, are adapted to cause the one or more computing devices to perform deconvolution of bulk or spatial RNA sequencing data for a sample by performing a method comprising the steps of: The method further comprising: i) obtaining input from sources including bulk or spatial RNA-seq data, single cell RNA-seq data, and cell type annotations for the samples; ii) generating a single cell UMI (unique molecular identifier) ​​count matrix based on the bulk or spatial RNA-seq data, the single cell RNA-seq data, and cell type annotations, where each column of the single cell UMI count matrix corresponds to a cell in the sample or a spatial location in the sample, and each row of the single cell UMI count matrix corresponds to a gene in the sample; iii) selecting from the single cell UMI count matrix a subset of the most variably expressed genes based on the population of genes whose expression is enriched in each cell type in the single cell UMI count matrix, or the population of genes whose expression varies above a predetermined level among all cells in the single cell UMI count matrix, or both; iv) estimating mean and variance parameters of the genes selected from the single-cell UMI count matrix by fitting a negative binomial distribution to each gene in a cell type; v) calculating the cross-cell type specificity of each of the selected genes based on the estimated mean and variance parameters; vi) estimating cross-sample genetic variability from the bulk or spatial RNA-seq data and / or single cell samples based on the ratio of mean gene expression levels and variances of the selected genes; vii) estimating gene-by-gene scaling coefficients using both the bulk or spatial RNA-seq data and the single-cell RNA-seq data using a regression-based approach; viii) generating a regression model based on the cross-cell type specificity of each selected gene, the cross-sample gene variability from the bulk or spatial RNA-seq data or the single cell samples, and a scaling coefficient for each gene; ix) estimating cell type proportions in the bulk or spatial RNA sequencing data using the regression model, wherein the regression model gives higher weight to selected genes with high cell type specificity and low cross-sample variability; wherein the computer-readable medium comprises a computer readable medium for performing a method of inferring a percentage of cell types in the bulk or spatial RNA sequencing data.

25. One or more processors; and a memory having processor-executable instructions for execution by the one or more processors, The processor-executable instructions, when executed by the one or more processors, are adapted to cause the one or more processors to perform a method comprising the steps of: deconvolving bulk or spatial RNA sequencing data for a sample; The method further comprising: i) obtaining input from sources including bulk or spatial RNA-seq data, single cell RNA-seq data, and cell type annotations for the samples; ii) generating a single cell UMI (unique molecular identifier) ​​count matrix based on the bulk or spatial RNA-seq data, the single cell RNA-seq data, and cell type annotations, where each column of the single cell UMI count matrix corresponds to a cell in the sample or a spatial location in the sample, and each row of the single cell UMI count matrix corresponds to a gene in the sample; iii) selecting from the single cell UMI count matrix a subset of the most variably expressed genes based on the population of genes whose expression is enriched in each cell type in the single cell UMI count matrix, or the population of genes whose expression varies above a predetermined level among all cells in the single cell UMI count matrix, or both; iv) estimating mean and variance parameters of the genes selected from the single-cell UMI count matrix by fitting a negative binomial distribution to each gene in a cell type; v) calculating the cross-cell type specificity of each of the selected genes based on the estimated mean and variance parameters; vi) estimating cross-sample genetic variability from the bulk or spatial RNA-seq data and / or single cell samples based on the ratio of mean gene expression levels and variances of the selected genes; vii) estimating gene-by-gene scaling coefficients using both the bulk or spatial RNA-seq data and the single-cell RNA-seq data using a regression-based approach; viii) generating a regression model based on the cross-cell type specificity of each selected gene, the cross-sample gene variability from the bulk or spatial RNA-seq data or the single cell samples, and a scaling coefficient for each gene; ix) estimating cell type proportions in the bulk or spatial RNA sequencing data using the regression model, wherein the regression model gives higher weight to selected genes with high cell type specificity and low cross-sample variability; and thereby a method for inferring the percentage of cell types in the bulk or spatial RNA sequencing data.

Images