Time-varying characteristic predictor

Abstract

A computer-implemented method is provided for obtaining a predictor for predicting a time-varying property based on gene transcription data. The method includes receiving a dataset including data samples obtained from each cell sample having different values ​​of the time-varying property. Each data sample includes several transcription levels and an actual value of each time-varying property of the cell sample of each data sample, each transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene. An embedded dataset is generated for each data sample including an embedding sample. The dimensionality of the embedding sample is less than the number of the several transcription levels. The embedding dataset is applied as input to the predictor to generate a predicted value of the time-varying property of each embedded sample. The predictor's prediction coefficients are adjusted to reduce an error measure of the prediction error between the predicted value and the actual value of each time-varying property to obtain the predictor. The time-varying property may be age, such as biological age, or progression of a disease, disorder, or condition.

Description

[Technical field]

[0001] The present disclosure relates to a computer-implemented method, and corresponding computer program product, computer-readable medium and system, for obtaining a predictor for predicting time-varying characteristics from gene transcripts.In particular, but not exclusively, the present disclosure relates to a predictor of the biological or chronological age or progression of disease, disorder or condition.Another non-limiting aspect of the present disclosure relates to estimating the contribution of different genes or gene transcripts to prediction. [Background technology]

[0002] Many diseases, including Parkinson's, Alzheimer's and osteoarthritis, have a component of aging. There is growing interest in identifying ways to induce cell and tissue regeneration using novel therapeutics that can unlock the regenerative potential present in very young cells. Over the past five years, great advances have been made in the highly promising science of a field called partial epigenetic reprogramming.

[0003] The only global method known so far to change cells to a younger state was to generate iPS cells (e.g. using Yamanaka factors). Unfortunately, cells that undergo this pluripotent transition also change their identity, so this technique cannot be used to create anti-aging drugs or to extend health lifespan.

[0004] Now that partial epigenetic reprogramming studies have shown that age-reversal factors can be separated from cell identity factors, efforts are underway to translate this process into clinical practice.

[0005] Aging clocks are an elegant way to understand how to promote the cellular rejuvenation process. The first aging clock was developed by Horvath et al. (see, for example, US20160222448A1 and US20190185938A1) and is called the epigenetic clock because it is based on methylation levels. Although epigenetic clocks predict age very accurately, they have some limitations, such as the difficulty of biological inference and the current inability to validate or target individual sites to validate potential therapeutic effects. Therefore, attention has been focused on transcriptome clocks, which predict age based on RNA expression levels. Transcriptome clocks are described, for example, in US10325673B2 and Holzscheck et al. (npj Aging Mech Dis 7, 15 (2021)). However, a major feature of these transcriptome clocks is that they operate at the summarized transcriptional level of the corresponding gene pathways, and therefore, in order to create such a clock, knowledge of the gene pathways is required in advance. The inventors realized that this has many drawbacks, as explained below. Thus, there is a need in the art for a CLOCK (a predictor of aging) that overcomes these limitations. Summary of the Invention

[0006] Aspects of the invention are set out in the accompanying independent claims. Optional features of some embodiments are set out in the dependent claims.

[0007] The present disclosure provides a computer-implemented method for obtaining a predictor for predicting a time-varying characteristic (a characteristic that changes over time) based on gene transcription data (i.e., RNA expression levels). The aging clock is an example of a predictor of a time-varying characteristic (age), but it will be understood that the present disclosure is not limited to age as the time-varying characteristic and is applicable to other time-varying characteristics.

[0008] The method includes receiving a data set including data samples obtained from respective cell samples having different values ​​of a time-varying characteristic. The cell samples may be single cells or a collection of cells whose transcription levels have been pooled to form the data sample. The cell samples may be obtained, for example, from cell culture in a test tube. Alternatively, the cell samples may be obtained from an individual, for example, by biopsy. Obtaining the cell samples is typically not part of the method. Each data sample includes several transcription levels. Each data sample further includes respective actual values ​​of the time-varying characteristic of the cell samples of each data sample. The time-varying characteristic may be biological or chronological age, progression or stage of a disease or condition, such as cancer or a neurodegenerative condition, such as Alzheimer's or Parkinson's disease. Thus, it will be appreciated that although each cell sample has a different value of the time-varying characteristic, each cell sample may all have been obtained at the same time, for example representing different stages in the progression of a disease or condition. The time-varying characteristic may be of one or more organisms or subjects from which the cell samples were obtained.

[0009] Each transcription level is the transcription level of individual gene transcripts, or the pooled transcription level of gene transcripts of individual genes.Therefore, transcription level can be obtained from the corresponding transcription count of individual genes or gene transcripts in each cell sample.For example, in some embodiments, transcription count can be obtained by using transcriptome technology such as RNA-Seq.

[0010] Because this method works on individual gene transcripts or genes, biases associated with pathway definition and selection can be avoided. Moreover, in some embodiments, new genes involved in causing time-varying properties can be discovered. No knowledge of gene pathways or biological activity is required, and therefore, unlike prior approaches in the art, the analysis can use transcription levels derived from transcription counts of gene transcripts in cell samples without knowledge of gene pathways or biological activity.

[0011] The method includes generating an embedded dataset comprising an embedding sample for each data sample from each transcription level. The dimensionality of the embedding sample is less than the number of transcription levels such that the embedding provides dimensionality reduction. In some embodiments, the dimensionality of the embedded sample can be selected based on the respective predictive performance of the embedding datasets having different respective dimensionalities. Advantageously, reducing the dimensionality can improve computational efficiency and reduce the amount of variance due to technical noise. This is particularly advantageous for single cell samples where technical noise may be large compared to the biological signal.

[0012] In some embodiments, the method can include applying a transformation to the dataset to generate the embedded dataset. The transformation can be obtained by operating on the dataset, for example, by operating on a covariance matrix of the dataset. In some embodiments, the transformation can be obtained without knowledge of gene pathways.

[0013] In some embodiments, the embedding may involve a linear transformation of the transcriptional dataset to generate an embedded dataset, and in some specific embodiments, the embedded dataset comprises a subset of the principal components of the transcriptional dataset. In some embodiments, a non-linear mapping may be used.

[0014] The method may, in some embodiments, include applying an inverse mapping to the predictive coefficients to project them to the dimensions of the transcriptional dataset. The inverse mapping maps from the embedded cell samples to the corresponding cell samples. In this way, a measure of the contribution to predicting the value of the time-varying characteristic can be derived for each gene or transcript. In some cases, the (possibly approximate) inverse mapping of the transformation can be used to (at least approximately) project the predictive coefficients to the dimensions of the dataset. In the case of a linear transformation, the inverse mapping is the inverse matrix found by inverting the matrix. In cases such as PCA, the eigenvectors of the matrix of eigenvectors are orthogonal, so the inverse mapping may be the transpose of the linear transformation or the linear transformation itself. In some embodiments, the transformation may be nonlinear, and the inverse operation of the transformation, the inverse mapping from the embedded data samples to the corresponding data samples, may be used to at least approximately project or transform the predictive coefficients to the dimensions of the dataset. The inverse mapping may be an approximation found, for example, by numerical optimization. The inverse mapping of the coefficients can serve as a measure of the importance of the dimension of the transcriptional dataset, i.e., the importance of each corresponding gene or transcript to prediction. Thus, the inverse mapping can be used to guide data-driven discovery of genes or transcripts involved in driving contribution to biological age, chronological age, and / or disease prediction. The coefficients of each gene or transcript can be aggregated in gene set enrichment analysis to guide discovery of biological pathways, processes, and functions that contribute to biological age, chronological age, and / or disease prediction.

[0015] The embedded dataset is then applied as input to a predictor to generate predicted values ​​of the time-varying properties for each embedded sample, and the prediction coefficients of the predictor are adjusted to reduce a measure of prediction error between the predicted and actual values ​​of each of the time-varying properties. In some embodiments, the predictor is also obtained without using gene pathway or biological activity information. In some embodiments, the predictor is initially obtained in this manner and then refined using prior knowledge of gene pathways or biological activities or biological knowledge obtained from the prediction coefficients of the predictor.

[0016] In some embodiments, the embedded dataset may be scaled to have substantially constant variance across dimensions. This enhances the initial contribution of low variance dimensions of the embedded dataset to the adjustment of the prediction coefficients, as opposed to, for example, unweighted PCA regression. The inventors have realized that high variance components are not necessarily related to time-varying properties, but may represent other sources of biological or technical variation. By weighting the variability of all components equally, low variance components have an equal starting point in the regression optimization, facilitating the discovery of biologically relevant components.

[0017] In some implementations, the predictor is a linear predictor. Advantageously, this allows the prediction coefficients to be easily interpretable, for example as described above. The linear predictor may, in some implementations, include a regularization technique to promote sparsity in the prediction coefficients, which may further aid in interpretability since fewer coefficients make a significant contribution to the prediction. For example, in some implementations, the adjustment of the prediction coefficients comprises an elastic net regression. In some implementations, the prediction error may be the median of the absolute prediction errors.

[0018] Some embodiments include receiving a further dataset, the further dataset including further data samples obtained from respective further cell samples having different values ​​of the time-varying characteristic, each further data sample including several further transcription levels and respective further actual values ​​of the time-varying characteristic of the further cell samples for each further data sample. The further transcription levels are derived from further transcription counts of gene transcripts in the further cell samples without using knowledge of gene pathways, as described above. These embodiments further include transforming the dataset and the further dataset into a common dataset including the data sample and the further data sample, where transforming the dataset and the further dataset includes reducing variability of the data sample and the further data sample that are not common to the dataset and the further dataset. In these embodiments, generating the embedded dataset includes generating an embedded sample for each data sample in the common dataset.

[0019] In some embodiments, the labeled dataset can be used to predict the properties of an unlabeled dataset to predict the time-varying properties of the new dataset. In such embodiments, a further dataset is also received, but this time does not include the time-varying property values. Again, converting the dataset and the further dataset into a common dataset including the data samples and the further data samples includes reducing the variability of the data samples and the further data samples that are not common to the dataset and the further dataset, and generating the embedded dataset includes generating an embedded sample for each data sample in the common dataset. In these embodiments, applying the embedded dataset as input to the predictor includes applying only the embedded samples corresponding to the gene transcription data samples as input to generate the respective predicted values ​​of the time-varying properties for the embedded data samples corresponding to the data samples. After obtaining the predictor, these embodiments include applying the embedded samples corresponding to the further data samples to the predictor to predict the respective values ​​of the time-varying properties for the further cell samples.

[0020] The described embodiments may further include generating a report identifying values ​​of the time-varying properties of one or more individuals or subjects from which the cells were derived, and / or indicators of the predictive contribution of the genes or transcripts in the dataset. The report may be in any suitable form, for example stored digitally on one or more storage media, displayed on a display screen, and / or printed on paper or other suitable medium.

[0021] The present disclosure further extends to a computer program product comprising computer code instructions which, when executed on a processor, perform the methods described, in particular the methods as set forth above, such as one or more computer readable media comprising such computer code instructions.The present disclosure further extends to a system comprising a processor and such a computer readable medium, the processor being configured to execute the computer code instructions, and to a system comprising means for performing the methods described, in particular the methods as set forth above.

[0022] The following description of specific embodiments is made with reference to the accompanying drawings to illustrate at least one disclosed invention. The headings are provided for clarity of description only and are not to be used in interpreting the contents of the disclosure. [Brief description of the drawings]

[0023] [Figure 1] FIG. 1 illustrates a computer-implemented method for determining a prediction of a time-varying property of a cell sample. [Diagram 2] FIG. 2 illustrates a computer-implemented method for obtaining a predictor of a time-varying property of a cell sample, which involves merging data sets corresponding to different batches of cell samples. [Diagram 3] FIG. 3 shows a computer-implemented method for obtaining a predictor for values ​​of a time-varying property of a cell sample, which includes merging datasets corresponding to different batches of cell samples and using a predictor obtained from one batch to predict the time-varying property of another batch. [Figure 4] FIG. 4 illustrates an example hardware implementation suitable for carrying out the disclosed methods. [Diagram 5] Figure 5 shows the distribution of median absolute error (MAE) of predicted ages of the aging clock trained with varying numbers of principal components as input. [Figure 6A]FIG. 6A shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6B] FIG. 6B shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6C] FIG. 6C shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6D] FIG. 6D shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6E] FIG. 6E shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6F] FIG. 6F shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 6G] FIG. 6G shows the performance of a gene expression directly trained clock ("Expression Clock") and a clock trained according to the method described herein ("RD Clock") in single cells of a test set from various mouse organs. [Figure 7] FIG. 7 shows the average time taken for one iteration of clock training versus the number of cells used in the training process. [Figure 8] Figures 8A-C show a comparison of performance metrics for different aging clocks. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0024] Referring to FIG. 1, in step 110, a gene transcription dataset is received for a batch of cell samples. The cell samples and datasets may have been obtained in any suitable manner, for example as described above. The dataset is generated from raw gene transcription counts (individual transcript counts or sums per gene), resulting in one count per transcript or gene, cell sample, and measurement time point, resulting in one expression vector of expression levels per cell sample and measurement time point. The data samples may be obtained from counts of single cells (single cell samples) or from pooled counts of multiple cell samples. The dataset is or has been processed to derive expression levels using conventional numerical conditioning of the count data, such as normalizing the data, log-transforming the data, and normalizing the log-transformed data, for example to have a mean of zero and a unitary standard deviation. The overall scale factor may of course be varied at will. What is important is that each transcript level is an individual gene transcript, or a pool of gene transcripts for individual genes, e.g. summed over that gene. As a result, no prior knowledge of gene pathways or biological activity is required in the process to generate the expression vectors. Moreover, the subsequent described processing can be performed without prior knowledge of gene pathways or biological activity, although it will be appreciated that additional processing steps, such as further post-regression adjustment of the predictive coefficients, may use such prior knowledge, as well as biological knowledge, such as that obtained from the regression itself, about how predictive a particular gene is for a time-varying characteristic, which may be derived from the predictive coefficients, for example, as described below.

[0025] The term "gene pathway" refers to a network of genes that function together to carry out a certain biological process. Such biological processes are also called "biological pathways". That is, it is a series of interactions between molecules in a cell that results in a biological effect, e.g. a change in the cell or the production of a product. These molecules are coded by genes, and thus it can be seen that the result of a network of genes in a gene pathway is a biological pathway. Knowledge about gene pathways and biological pathways can be obtained, for example, from the "Hallmark" pathway collection (Liberzon, A. et al. The molecular signatures database hallmark gene set collection. Cell Syst. 1, 417-425 (2015)) and publicly available databases such as the KEGG pathway database (https: / / www.kegg.jp / ).

[0026] The resulting gene transcription dataset is organized (or received) as a matrix E with one transcription vector per cell sample and time point as row vector. The eigenvector matrix W and the diagonal matrix of eigenvectors (which are the variances of E transformed by the basis vectors W) are found using any suitable technique, such as eigenvalue decomposition or, more typically, singular value decomposition.

[0027]

number

[0028] In step 120, the k-column eigenvectors (or principal components) w i matrix v

number

[0029]

number

[0030] k can be selected as an appropriate value, with higher values ​​requiring more computation but potentially containing more biologically relevant information. k=50 has been found to be a suitable maximum in most settings, and in some embodiments k is between 20 and 30, for example. k can also be selected in an iterative manner by comparing the performance of the coefficient fitting described below for different values ​​of k, and selecting the value that achieves the best or at least satisfactory performance. In some embodiments, instead of selecting the component with the largest eigenvalue k, components can be selected based on performance as described above, for example, by selecting components in the middle of a range of eigenvalues ​​or at a particular ordinal number of eigenvalues.

[0031] It should be noted that each data sample in the dataset further comprises an actual measurement of a time-varying property of the cell sample (or the organism from which the cell sample was derived), although the dataset comprises multiple cell samples at multiple time points, one such value per cell sample and time point. The actual value may be measured at each time point, for example by measuring a quantity that correlates with a biomarker indicative of biological age or disease trajectory, or may be otherwise known for the organism, such as disease progression or stage, or may simply be the time point itself, as in the case of chronological age. Senescence clock measurements, such as epigenetic clock measurements, can be used as biomarkers indicative of biological age. Other time-varying properties of the cell sample or the organism from which it is derived can also be used in addition to biological age, chronological age, or stage or progression of a disease, e.g., a neurodegenerative disease such as Alzheimer's or Parkinson's.

[0032] The actual values ​​are organized or received in a column vector y with the same number of rows as E, one for each data sample.

[0033]

number

[0034] In step 130, the embedded data set is applied to a predictor to obtain values ​​y of the time-varying characteristic y. * In step 140, the regression is trained on the embedding dataset X by predicting y and y by adjusting the prediction coefficients and offsets β in a vector β that contains linear weights for the principal components in the regression. * The coefficients are adjusted to minimize the error between the mean squared error

number

number

[0035] A variety of minimization methods can be used, including simple least squares regression. In some implementations, it has been found to be advantageous to use elastic net linear regression (Zou, H., & Hastie, T. (2005) Regularization and variable selection via the elastic net; Journal of the Royal Statistical Society:Series B (Statistical Methodology), 67(2), 301-320; https: / / doi.org / 10.1111 / j.1467-9868.2005.00503.x, incorporated herein by reference. This paper also discusses several alternative regression approaches that can be used in some embodiments). Advantageously, elastic net regression promotes sparsity of predictive coefficients, i.e., coefficients tend to be small for most coefficients, with the mass of coefficients concentrated in the more predictive regressors (in this case, the more predictive principal components). This facilitates interpretability of the principal components (and therefore the corresponding transcript levels) in terms of biological relevance in the processes captured by the time-varying properties.

[0036] The training and tuning of the coefficients may be performed in any suitable manner. To mitigate the overfitting problem, it is advantageous to train the parameters using n-fold cross-validation. Additionally, some data may be retained as pure test data to evaluate the model performance on unknown data. Any linear predictor may be used depending on the particular implementation, and the embedding and regression steps may be combined. One linear predictor that may be used is partial least squares or a variant thereof, including embedding of both E and y. However, the present disclosure is not limited to linear predictors, and other predictors such as feedforward or recurrent neural networks may be used to provide predictors of the values ​​of the time-varying properties. It is noted that linear predictors are advantageous not only due to their algorithmic simplicity and efficiency, but also due to the interpretability of the prediction coefficients, as described below.

[0037] To assess the contribution of each expression level to the prediction of a time-varying trait, i.e., to determine which expression levels are more predictive than others, the prediction coefficients are calculated as contribution coefficients β * can be projected onto

[0038]

number

[0039] where R "unscales" the coefficients to compensate for the scaling by S during the regression. Thus, β * The elements of provide a measure of how predictive the gene or transcript corresponding to each transcription level is for the time-varying characteristic.

[0040] In an optional step 150, new transcription samples can be received and the trained predictor can be used to predict values ​​of the time-varying properties of the new transcription samples. The new transcription samples may be samples obtained from the same experiment / event, or set of experiments / events, used to obtain the samples used for training, where values ​​of the time-varying properties are not available, or the new transcription samples may be newly obtained. For good predictions, the conditions under which the newly obtained samples are obtained must be carefully controlled to match the conditions under which the training samples were obtained, to avoid significant batch effects due to differences in technical noise that degrade the prediction performance. In many cases, this is difficult, and below we describe how to compensate for such batch effects, i.e., by adding new training data to existing training data, or by combining new unlabeled data with one or more training datasets to improve prediction performance.

[0041] Further, in optional step 160, the β for each gene / transcript is calculated to assess predictive potential. * A report can be generated that provides one or both of the elements of the regression model and, if applicable, provides predicted values ​​of the time-varying properties for one or more new data samples. Other elements of the report may be regression coefficients or other measures of goodness of fit, residuals and / or other quantities that facilitate interpretation of the data and predictors.

[0042] The process for training a predictor and making predictions using the combined datasets includes step 210 receiving a first gene transcription dataset E and a second (further) gene transcription dataset

number

number

[0043] where {·|·} is a combination operation on datasets, which in its simplest implementation is the concatenation of two datasets.

[0044] In some embodiments, the combination process consists of a suitable normalization of the individual data sets, e.g., expression levels are normalized by their cosine norms calculated for each cell sample.

number

[0045] In some embodiments, the combining of the data sets includes correcting for differences between the data sets of different batches (typically due to technical noise). In some embodiments, a batch correction vector is subtracted from each data sample of the second batch, or in terms of the batch correction matrix B of the batch correction row vector:

[0046]

number

[0047] An embedded data set X may then be formed in step 220 in a manner similar to step 120 as follows:

[0048]

number

[0049] Where:

number

number

[0050] In some implementations, combining the data sets in step 214, Equation 6 and Equation 7, includes transforming the data sets to a different coordinate system. In a specific implementation, the principal components of the combined data set are determined and the combined data set is divided into two sets of data with the largest eigenvalues ​​K

number

number

[0051] Computing the principal components of the combined dataset involves centering the mean of each dataset being combined (rather than just the average of the combined dataset) and weighting the contribution of each cell sample to the covariance matrix by the inverse of the number of cell samples in each dataset being combined (or equivalently, using the average of the covariance matrices of the datasets being combined as the covariance matrix for the principal components analysis). Principal components are computed for the combined dataset in a conventional manner, such as using eigenvalue decomposition or singular value decomposition.

[0052] Batch correction then proceeds as described above with reference to Equations 6 and 7 for the selected principal components of the combined dataset. In these implementations, the dimensionality of the combined dataset remains orthogonal throughout the batch correction and can be used to form the combined embedded dataset X using Equation 8, although this is not required and the same selected k-dimensional C can be used.

[0053]

number

[0054] Here, V i,i are the non-zero diagonal entries of the covariance matrix V of C. Of course, one can choose C with smaller dimensionality.

[0055] Various approaches for calculating the batch correction vector B are known and can be used in embodiments. In some embodiments, a mutual nearest neighbor (MNN) approach is used, see Haghverdi, L., Lun, ATL, Morgan, MD, & Marioni, JC (2018) Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors; Nature Biotechnology, 36(5), 421-427; https: / / doi.org / 10.1038 / nbt.4091 and https: / / marionilab.github.io / FurtherMNN2018 / theory / description.html, each of which is incorporated herein by reference.

[0056] MNN was the first

number

number

[0057] if

number

number

number

[0058] K is chosen empirically or empirically for each dataset, with a larger number of nearest neighbors increasing robustness to noise and sampling nearest neighbors deeper into each cell sample cloud, but at increased computational cost. In practice, K=20 is a good choice.

[0059] In embodiments where MNNs are found directly based on expression levels without orthogonalization and / or dimensionality reduction, such as PCA, as described above, the MNN batch correction vector for the MNNs is the difference vector

number

[0060] The above results in a set of batch correction vectors for the MNN, i.e., the MNN batch correction vectors. The batch correction vectors for other non-MNN data samples are then derived from the MNN batch correction vectors, for example by combining them with a Gaussian kernel, using another form of weighted averaging, or by using only the nearest MNN batch correction vectors for each cell sample. This results in locally varying batch correction vectors for all data samples, which are then used in Equation 7 above.

[0061] In some embodiments, the cell samples of each batch are projected onto a respective bisection plane perpendicular to the average vector of the MNN batch vectors of each data set (but adjusted for the MNN cell sample's own projection) before applying the MNN batch vector as described above. This ensures that the merged cell samples will intermingle and not simply be brought together as touching clouds, even if K is not large enough to sample nearest neighbors across the expected opposing faces of the batch. Alternatively, the cell samples of the merged data set after batch correction can be projected onto a common bisection plane perpendicular to the average MNN batch correction vector, or this step can be omitted if the value of K is large enough.

[0062] For complete details of the batch correction method along these lines, see Haghverdi et al. (2018) and supplementary information and software packages are available as part of the batchelor R package, also described in Haghverdi et al. (2018). For further related implementations of the above mentioned method of squashing variability along the mean MNN batch correction vector for each dataset before applying the batch correction vectors, see https: / / marionilab.github.io / FurtherMNN2018 / theory / description.html. Alternative batch calibration methods that output a dimensionality-reduced embedding of the calibration dataset can equally be used, for example Seurat v3 implements a canonical correlation analysis before identifying “anchors” in a similar manner as above.

[0063] 3, the steps of receiving the first and second gene transcription datasets 310, 312, generating the combined dataset 314, generating the combined embedding dataset 320, training the predictor 330, and adjusting the prediction coefficients 340 are similar to steps 210, 212, 214, 220, 230, and 240 described above, with the corresponding disclosures applying where appropriate, except that only the first gene transcription dataset contains the actual value y of the time-varying property, and this information is not received (or is ignored) along with the second gene transcription data in step 312. Thus, in steps 330 and 340, the predictor is trained using only data samples from the first dataset for which property values ​​are available, the prediction coefficients are adjusted, and the resulting predictor is used to predict the respective values ​​of the property for data samples of the second dataset. In this way, unknown values ​​of the property can be predicted, for example, for samples obtained from new individuals of an organism for which the time-varying property is unknown. The step 360 of generating a report is similar to step 160 described above and includes predicted values ​​for the samples of the second data set.

[0064] The described embodiment calculates the embeddings using principal component analysis, implemented for example with SVD, and selects the number of principal components for dimensionality reduction. Other methods of finding embeddings are equally applicable in various embodiments, and PCA can be substituted for the embeddings. For example, the embeddings can be found using kernel methods, such as kernel PCA (kPCA) or nonlinear methods such as training an autoencoder (AE). kPCA applies eigenvalue decomposition or SVD to a kernel matrix derived from the data using a kernel function in a similar manner as PCA does for a covariance matrix. The gene prediction coefficients can be found and restored using inverse mapping to find the weights in the gene space in a similar manner to PCA described above. The inverse mapping is found by numerical optimization, so that the gene prediction coefficients are at least approximately restored. The AE is a neural network trained to match the input and output, and is composed of a hidden embedding layer with fewer units than the input and output layers that perform the embedding. The gene prediction coefficients can be at least approximately restored from the embedded prediction coefficients using a trained decoding network between the hidden embedding layer and the output layer. In general, at least approximate gene prediction coefficients can be obtained from the embedded prediction coefficients by applying an inverse mapping of the embedding transformation to the embedded prediction coefficients. The inverse mapping may correspond to a mathematical inverse transformation or other operation that maps from the embedding space to the gene space, i.e., from the embedded data samples to the corresponding data samples. This projection onto the dimensions of the (non-embedded) data set may be approximate (e.g., found by numerical methods or by training a neural network) or mathematically exact (e.g., found by inverse transformation or transposition of a matrix, as in the case of embedding PCA described in detail above).

[0065] Hardware implementation example FIG. 4 illustrates a block diagram of one implementation of a computing device 400 upon which a set of instructions may be executed to cause the computing device to perform any one or more of the methodologies discussed herein. In alternative implementations, the computing device may be connected (e.g., networked) to other machines in a local area network (LAN), an intranet, an extranet, or the Internet. The computing device may operate in the capacity of a server or a client machine in a client-server network environment, or as a peer machine in a peer-to-peer (or distributed) network environment. The computing device may be a personal computer (PC), a tablet computer, a set-top box (STB), a personal digital assistant (PDA), a mobile phone, a web appliance, a server, a network router, a switch, a bridge, or any machine capable of executing a set of instructions (sequential or otherwise) that specify actions to be performed by that machine. Additionally, although only a single computing device is illustrated, the term "computing device" shall be deemed to include any collection of machines (e.g., computers) that individually or jointly execute a set (or sets) of instructions to perform any one or more of the methodologies discussed herein.

[0066] The exemplary computing device 400 includes a processing device 402, a main memory 404 (e.g., read only memory (ROM), flash memory, dynamic random access memory (DRAM), such as synchronous DRAM (SDRAM) or Rambus DRAM (RDRAM)), a static memory 406 (e.g., flash memory, static random access memory (SRAM), etc.), and a secondary memory (e.g., a data storage device 418), which communicate with each other via a bus 430.

[0067] The processing device 402 represents one or more general-purpose processors, such as a microprocessor, a central processing unit, etc. More specifically, the processing device 402 may be a complex instruction set computing (CISC) microprocessor, a reduced instruction set computing (RISC) microprocessor, a very long instruction word (VLIW) microprocessor, a processor implementing other instruction sets, or a processor implementing a combination of instruction sets. The processing device 402 may also be one or more special-purpose processing devices, such as an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), a digital signal processor (DSP), a network processor, etc. The processing device 402 is configured to execute processing logic (instructions 422) for performing the operations and steps discussed herein.

[0068] The computing device 400 may further include a network interface device 408. The computing device 400 may also include a video display unit 410 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)), an alphanumeric input device 412 (e.g., a keyboard or touch screen), a cursor control device 414 (e.g., a mouse or touch screen), and an audio device 416 (e.g., speakers).

[0069] The data storage device 418 may include one or more machine-readable storage media (more specifically, one or more non-transitory computer-readable storage media) 428. Stored in 428 are a set of one or more instructions 422 embodying any one or more of the methodologies or functions described herein. The instructions 422 may also reside, completely or at least partially, within the main memory 404 and / or within the processing device 402 during execution thereof by the computer system 400, with the main memory 404 and the processing device 402 also constituting computer-readable storage media.

[0070] The various methods described above can be implemented by a computer program. The computer program can include computer code arranged to instruct a computer to perform one or more functions of the various methods described above. The computer program and / or code for performing such methods may be provided to an apparatus such as a computer on one or more computer readable media, or more generally on a computer program product. The computer readable medium may be transitory or non-transient. The one or more computer readable media may be, for example, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, or a propagation medium for data transmission, for example, for downloading code via the Internet. Alternatively, the one or more computer readable media may take the form of one or more physical computer readable media, such as semiconductor or solid state memory, magnetic tape, removable computer diskette, random access memory (RAM), read only memory (ROM), rigid magnetic disk, optical disk such as CD-ROM, CD-R / W, DVD, etc.

[0071] In implementation, the modules, components, and other functions described herein may be implemented as distributed components or may be integrated into the functionality of a hardware component such as an ASICS, FPGA, DSP, or similar device.

[0072] A "hardware component" is a tangible (e.g., non-transient) physical component (e.g., a set of one or more processors) that can perform specific operations and may be configured or arranged in a specific physical manner. A hardware component may include dedicated circuitry or logic that is permanently configured to perform specific operations. A hardware component may be or include a dedicated processor, such as a field programmable gate array (FPGA) or an ASIC. A hardware component may also include programmable logic or circuitry that is temporarily configured by software to perform specific operations.

[0073] Accordingly, the phrase "hardware component" should be understood to encompass a tangible entity that is physically configured, permanently configured (e.g., hard-wired), or temporarily configured (e.g., programmed) to operate in a particular manner or perform particular operations as described herein.

[0074] Additionally, the modules and components may be implemented as firmware or functional circuitry within a hardware device. Additionally, the modules and components may be implemented in any combination of hardware devices and software components, or may be implemented solely in software (e.g., code stored or otherwise embodied in a machine-readable or transmission medium).

[0075] Unless otherwise indicated, and as will be apparent from the discussion that follows, discussions throughout this specification using terms such as "receiving," "determining," "comparing," "enabling," "maintaining," "identifying," "obtaining," "receiving," "generating," "applying," "adjusting," "producing," "scaling," "deriving," and the like refer to operations and processes of a computer system, or similar electronic computing device, that operate on and transform data represented as physical (electronic) quantities in the computer system's registers and memory to other data represented as physical quantities in the computer system's memory or registers, or other information storage, transmission, or display device.

[0076] Working Example To demonstrate the aging clock method, we used Tabula Muris Senis (The Tabula Muris Consortium. A single-cell transcriptomic atlas). We analyzed single-cell gene expression data from a mouse model of aging (https: / / doi.org / 10.1038 / s41586-020-2496-, incorporated by reference 1). The data included the transcriptomes of cells from multiple tissues of mice of known age. The data, obtained with a microfluidic ("droplet") method, included four tissues (heart, lung, limb muscle, and spleen) with enough time points to attempt to train the clock. Among these tissues, the most common cell types (annotated by the original authors) were selected and outlined in Table 1. To prevent the influence of sex differences, we limited our analysis to male mice.

[0077] Table 1-Tabula Muris Table summarizing tissues with enough cells present in Senis to reliably train the single-cell aging clock and the cell types that contribute to it.

[0078] [Table 1]

[0079] The median absolute error (MAE) was found to be a suitable loss function for the training process described herein. Figure 5 shows an example distribution of MAE for predicted ages of aging clocks trained with various numbers of principal components from cardiac endothelial cells as input. For clocks trained with fewer principal components, the most recent ones were discarded. Ten iterations were trained for each number of principal components. The dotted line is the threshold for the near-optimal model, and the shaded box is the minimum number of principal components that could generate a clock with this performance. The threshold was calculated by first identifying the number of principal components that generated the clock with the lowest average MAE. The standard error of the MAE generated by the 10 clocks trained with this number of principal components was then added to the average MAE of this set of clocks.

[0080] Clock Test Figure 6A-G shows the performance of a clock trained directly on gene expression and the method described herein (referred to as "RD-Clock" in the figures) in single cells from a test set of mouse heart (A and B), lung (C and D), limb muscle (E and F), and spleen (G), with individual cell types labeled as follows:

[0081] [Table 2]

[0082] Each boxplot shows the distribution of predicted ages for cells from a single donor mouse; grouped boxplots correspond to mice of the same age and are jittered in the x-direction to aid visualization. The upper and lower hinges correspond to the 75th and 25th percentiles, respectively, and the middle hinge indicates the median. Whiskers are 1.5 mm from the outer hinges. * The interquartile range is extended, with points outside this range represented by black dots. The median absolute error (MAE) by cell is shown for each plot by month, as is the Pearson correlation coefficient (Cor.), y=x, represented by the black dashed line.

[0083] When trained and tested within a single dataset (i.e., no batch effects exist between training and test cells), as can be seen in Figures 6A-G, the error (measured by MAE) of the method described herein is similar to that of a clock trained directly on the top 2000 highly variable genes. However, because later principal components have been discarded from the clock described herein, the model is believed to contribute less from technical noise and to be less biased due to "overfitting". Thus, the accuracy of the clock described herein is less inflated than that of a direct gene expression clock.

[0084] Another advantage is the reduction in time required for clock training, which scales with the number of cells used for training, as can be seen in Figure 7. Figure 7 shows the average time for one iteration of clock training (either directly on gene expression ["Expression", squares] or by the method described herein ["RD", circles]) versus the number of cells used for the training process. Points are alphabetized according to the tissue and cell type used for training, and lines are fitted using linear regression. Inset: Clock values ​​are shown on the same x-axis as the main plot, but the y-axis is truncated to aid visualization. Training was performed using an AMD Ryzen 7 5800X 8-core processor (3.80 GHz) with 32 GB RAM.

[0085] This time reduction is significant given that the training process must be repeated thousands of times during training and optimization. For reference, for a "realistic" training set size of ~5000 cells (splenic B cells), the method described here is approximately 60 times faster.

[0086] Clock transfer between data sets The aging clock method can be used to predict the donor age of single cells in datasets with little or no prior age annotation. We demonstrated this again using Tabula Muris Senis, which also contains single cell expression data for the four tissues used previously, collected by another sequencing method based on fluorescence-activated cell sorting (FACS). Male mouse cells were collected at 3, 18, and 24 months of age. However, due to the lack of male mouse samples between 1 and 18 months of age in all tissues profiled by the droplet method, the 3-month-old cells were excluded from further analysis.

[0087] Figure 8A-C compares performance metrics of the aging clock described herein ("RD") and a clock directly trained on gene expression ("Expression"). In each panel, metrics are normalized to those of the clock described herein. A: Time taken per ELN training iteration on one dataset, B: MAE per cell when the clock was trained and tested on one dataset, C: MAE per cell when the clock was trained on one dataset and used for cell age prediction on another dataset. In C, the clock described herein was trained on a corrected PCA matrix generated by the MNN method for droplet cells and tested on FACS cells; the direct expression clock trained in Figure 6A-G was applied directly to FACS cells; the clock was also trained on an expression matrix reconstructed from the MNN-corrected PCA matrix ("Expression-Reconstructed").

[0088] As noted above, direct expression clocks previously trained on droplet data performed worse on average when applied to FACS data than the clocks described herein trained on droplet data and tested on FACS data after batch correction (Figure 8C). The average reduction in MAE for the clocks described herein was 37%. This improvement is likely to be greater for datasets where batch effects are more pronounced.

[0089] Also, overfitting to technical noise contributes to increased error when direct expression clocks are batch transferred compared to the clocks described herein. In general, overfitting reduces the generality of direct expression clocks and also the biological conclusions that can be drawn from them. This means that direct expression clocks trained on a dataset will perform poorly in biological replicates of that dataset, even in the absence of any batch effects (highly unlikely). Since the latter condition is difficult to meet, one way to examine the generality of these clocks is to use the output of batch correction. Since the output of the MNN method is a correction matrix similar to the corrected PCA matrix, a "corrected" gene expression matrix can be reconstructed from this matrix in a similar manner as described herein. It is important to note that due to the forced movement of cells in PCA space, the resulting "expression" matrix is ​​highly distorted and should generally not be used as a mathematical proxy for real gene expression. However, this matrix is ​​the only practical way in which the generality of clocks such as those described herein can be compared to direct expression clocks in the absence of batch effects. Under these conditions, the clock method described herein reduces error (average MAE reduction = 30%, FIG. 8C), suggesting a significant benefit from the lack of overfitting in the clock according to the present disclosure.

[0090] It is to be understood that the above description is intended to be illustrative, and not restrictive. Many other embodiments will be apparent to those skilled in the art upon reading and understanding the above description. Although the present disclosure has been described with reference to certain exemplary embodiments, it will be recognized that the disclosure is not limited to the described embodiments, but can be practiced with modification and alteration within the spirit and scope of the appended claims. Accordingly, the specification and drawings are to be regarded in an illustrative rather than a restrictive sense. The scope of the present disclosure should therefore be determined with reference to the appended claims, along with the full scope of equivalents to which such claims are entitled.

[0091] The disclosed aspects and embodiments include the following numbered clauses: 1. A computer-implemented method for obtaining a predictor for predicting a time-varying characteristic based on gene transcription data, comprising: receiving a dataset comprising data samples obtained from respective cell samples having different values ​​of a time-varying property, each data sample comprising several transcription levels and respective actual values ​​of the time-varying property of the cell sample for each data sample, each transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; generating an embedding data set including an embedding sample for each data sample, the embedding sample having a dimensionality less than a number of transcription levels; applying the embedded dataset as input to a predictor to generate predictions of the time-varying properties for each embedded sample; and adjusting the prediction coefficients to reduce an error measure of the prediction error between the predicted and actual values ​​of each of the time-varying properties, thereby obtaining a predictor.

[0092] 2. The method of claim 1, comprising applying a transformation to the dataset to generate the embedded dataset, further comprising obtaining the transformation by operating on the dataset. 3. The method of claim 2, comprising obtaining the transformation without using knowledge of gene pathways. 4. The method of claim 2 or 3, comprising obtaining the transformation by operating on the covariance matrix of the data set. 5. The method of any of the preceding clauses, comprising scaling the embedding dataset to have a substantially constant variance across dimensions. 6. The method of any of the preceding clauses, comprising applying a linear transformation to the transcription dataset to generate the embedded dataset. 7. The method of claim 6, wherein the embedding dataset comprises a subset of principal components of the transcription dataset. 8. The method of any of the preceding clauses, comprising applying an inverse mapping to the predictive coefficients that maps from the embedded data samples to the data samples to project the predictive coefficients onto the dimensions of the dataset, thereby deriving a measure of the contribution of each gene or transcript to predicting the value of the time-varying property. 9. A method according to clause 8 when dependent on claims 6 and 7, wherein the inverse mapping comprises a matrix inverse of the linear transformation. 10. The method of any of the preceding clauses, wherein the predictor is a linear predictor. 11. The method of claim 10, wherein the linear predictor includes a regularization method to promote sparsity in the prediction coefficients. 12. The method of any of the preceding clauses, wherein adjusting the prediction coefficients includes elastic net regression. 13. The method of any of the preceding clauses, wherein the forecast error is the median of the absolute forecast errors.

[0093] 14. A method according to any preceding clause, further comprising: receiving a further data set comprising further data samples obtained from respective further cell samples having different values ​​of the time-varying property, each further data sample comprising a number of further transcription levels and respective further actual values ​​of the time-varying property of the further cell samples for each further data sample, each further transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; reducing variance of the data samples and the further data samples that are not common to the data set and the further data set by transforming the data set and the further data set into a common data set comprising the data samples and the further data samples, A method, wherein generating the embedded data set includes generating an embedded sample for each data sample in the common data set.

[0094] 15. The method according to any one of clauses 1 to 14, further comprising: receiving a further dataset comprising further data samples obtained from respective further cell samples having different values ​​of the time-varying property, each further data sample comprising a number of further transcription levels, each further transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; reducing variance of the data samples and the further data samples that are not common to the data set and the further data set by transforming the data set and the further data set into a common data set comprising the data samples and the further data samples; generating the embedded data set includes generating an embedded sample for each data sample in the common data set; applying the embedding dataset as input to the predictor includes applying as input only the embedding samples corresponding to the gene transcription data samples to generate predictions of each of the time-varying properties for the embedding data samples corresponding to the data samples; After obtaining the predictor, the method applies embedding samples corresponding to the further data samples to the predictor to predict values ​​of each of the time-varying properties of the further cell samples.

[0095] 16. The method of any of the preceding clauses, wherein the dimensionality of the embedding samples is selected based on the predictive performance of each of the embedding datasets having different respective dimensionalities. 17. The method of any preceding clause, wherein the time-varying properties are those of one or more organisms or subjects from which the cell sample was obtained. 18. The method of clause 16, further comprising generating a report identifying values ​​of the time-varying property for one or more organisms or subjects. 19. The method of any one of clauses 1 to 16, wherein the cells are obtained from a cell culture. 20. The method of any of the preceding clauses, wherein the cell samples are single cell samples each of a single cell. 21. The method of any of the preceding clauses, wherein the time-varying characteristic is biological age. 22. The method of any one of clauses 1 to 20, wherein the time-varying characteristic is chronological age. 23. The method of any one of clauses 1 to 20, wherein the time-varying characteristic is the progression of a condition or disease. 24. The method of clause 23, wherein the condition or disease is a neurodegenerative disease. 25. The method according to clause 24, wherein the neurodegenerative disease is Alzheimer's disease. 26. The method according to clause 24, wherein the neurodegenerative disease is Parkinson's disease. 27. The method of any one of clauses 1 to 20, wherein the time-varying characteristic is a progression state of cancer.

[0096] 28. The method of any of the preceding clauses, wherein the transcription level and, if present, further transcription levels are derived from the transcript counts of gene transcripts in the cell sample without knowledge of the gene pathway. 29. The method of claim 28, comprising generating the embedding dataset without using knowledge of gene pathways. 30. The method of claim 28, comprising applying the embedded dataset to obtain a predictor without using knowledge of gene pathways. 31. The method of any of the preceding clauses, further comprising, after obtaining the predictor, improving the predictor using prior knowledge of gene pathways or biological activities, any other prior biological knowledge, or knowledge obtained from the prediction coefficients. 32. A computer program product comprising computer code instructions which, when executed on a processor, performs the method according to any of the preceding clauses. 33. One or more computer-readable media comprising computer code instructions which, when executed on a processor, perform the method of any one of clauses 1 to 31. 34. A system including a processor according to clause 33 and a computer readable medium, the processor being configured to execute computer code instructions. 35. A system comprising means for carrying out the method according to any one of clauses 1 to 31.

Claims

1. 1. A computer-implemented method for obtaining a predictor for predicting a time-varying trait based on gene transcription data, comprising: receiving a dataset comprising data samples obtained from respective cell samples having different values ​​of a time-varying property, each data sample comprising several transcription levels and for each data sample an actual value of each of the time-varying properties of the cell sample, each transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; generating an embedding data set including an embedding sample for each data sample, the embedding sample having a number of dimensions less than the number of transcription levels; applying the embedded data set as an input to the predictor to generate a prediction of the time-varying property for each embedded sample; adjusting prediction coefficients of said predictor to reduce an error measure of a prediction error between the predicted and actual values ​​of each of said time-varying properties; The method according to claim 1,

2. The method of claim 1 , comprising applying a transformation to the data set to generate the embedded data set, and further comprising operating on the data set to obtain the transformation.

3. 3. The method of claim 2, comprising obtaining said transformation without knowledge of a genetic pathway.

4. The method of claim 2 or 3, comprising obtaining the transformation by operating on a covariance matrix of the data set.

5. 5. The method of claim 1, further comprising scaling the embedded data set to have a substantially constant variance across dimensions.

6. The method of claim 1 , comprising applying a linear transformation to the transcription data set to generate the embedded data set.

7. The method of claim 6 , wherein the embedded data set comprises a subset of principal components of the transcription data set.

8. 8. The method of claim 1, comprising applying an inverse mapping from the embedded data samples to the data samples to project the prediction coefficients onto the dimension of the dataset, thereby deriving a measure of contribution to predicting the value of the time-varying feature for each gene or transcript.

9. The method of claim 1 , wherein the predictor is a linear predictor.

10. receiving a further data set comprising further data samples obtained from respective further cell samples having different values ​​of said time-varying property, each further data sample comprising several further transcription levels and respective further actual values ​​of said time-varying property of said further cell samples for each further data sample, each further transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; transforming the data set and the further data set into a common data set comprising the data samples and the further data samples, thereby reducing variance of the data samples and the further data samples that are not common to the data set and the further data set; 10. The method of claim 1, wherein generating the embedded data set comprises generating an embedded sample for each data sample in the common data set.

11. receiving a further data set comprising further data samples obtained from respective further cell samples having different values ​​of said time-varying property, each further data sample comprising several further transcription levels, each further transcription level being a transcription level of an individual gene transcript or a pooled transcription level of gene transcripts of an individual gene; transforming the data set and the further data set into a common data set comprising the data samples and the further data samples, thereby reducing variance of the data samples and the further data samples that are not common to the data set and the further data set; generating the embedded data set includes generating an embedded sample for each data sample in the common data set; applying the embedded data set as input to the predictor includes applying as input only the embedded data sample corresponding to the gene transcription data sample to generate a prediction of each of the time-varying properties for the embedded data sample corresponding to the data sample; 11. The method of claim 1, wherein after obtaining the predictor, the embedding samples corresponding to the further data samples are applied to the predictor to predict values ​​of each of the time-varying properties of the further cell samples.

12. 13. The method of claim 1, wherein the dimensionality of the embedding samples is selected based on respective predictive performance of embedding data sets having different respective dimensionalities.

13. 13. The method of claim 1, further comprising generating a report identifying values ​​of the time-varying properties for one or more organisms or subjects from which the cell sample was obtained, the time-varying properties being for the one or more organisms or subjects.

14. 14. The method of claim 1, wherein the cell samples are single cell samples of each single cell.

15. 15. The method of any one of claims 1 to 14, wherein the time-varying characteristic is biological age, chronological age, or progression of a condition or disease, optionally wherein the condition or disease is a neurodegenerative disease or cancer, and optionally wherein the neurodegenerative disease is Alzheimer's disease or Parkinson's disease.

16. 16. The method of any one of claims 1 to 15, wherein the transcription level and, if present, the further transcription level are derived from the transcription counts of gene transcripts in a cell sample without knowledge of gene pathways.

17. 17. The method of claim 16, comprising generating the embedded dataset without knowledge of gene pathways and / or applying the embedded dataset without knowledge of gene pathways to obtain the predictor.

Images