Prediction of gene expression

Abstract

Methods of predicting the gene expression of one or more genes in a sample are described. The methods comprise receiving sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases including methylated cytosine and hydroxymethylated cytosine at one or more positions in the genome; and predicting for each gene of the one or more genes and using the sequence data associated with the gene, the value of a gene expression metric associated with the gene, using a machine learning model trained to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a gene expression metric for the gene.

Description

[0001] Prediction of gene expression

[0002] Field of the Disclosure

[0003] The present invention relates to methods of predicting gene expression metrics from 6-letters sequence data and particularly, although not exclusively, to methods of predicting gene expression using sequencing data comprising information about the presence of methylcytosine and hydroxymethylcytosine in gene and / or enhancer regions, and methods of performing genetic tests and screening perturbations using said methods.

[0004] Background

[0005] DNA encodes information in both genetic and epigenetic bases. Epigenetic information is closely associated with transcription, ultimately influencing cell fate, disease etc. DNA methylation, the covalent addition of a methyl group to cytosine residues (primarily at CpG sites) is the most widely studied epigenetic mark. The relationship between DNA methylation and gene expression has been previously studied. For example, Lou et al. (20 / 4) presented an integrated analysis of whole-genome bisulfite sequencing and RNA sequencing data from human samples and cell lines. The authors built machine learning models using DNA methylation levels at different subregions of a gene to explain observed expression levels of the gene. The models constructed classified genes between four expression classes (low, medium low, medium high, high - corresponding to the first, second, third and fourth quartiles of an observed distribution of gene expression) with an accuracy of about 0.75 on average. Pongor et al. (2022) examined global DNA methylation in human small cell lung cancer acquired using high-resolution DNA methylation arrays and studied the relationship between gene body methylation, promoter DNA methylation and gene expression, finding that highly expressed genes are hypomethylated on their promoters and hypermethylated in their gene bodies, whereas non-expressed genes fall in 3 distinct groups: genes methylated on the promoter and gene body, genes with no methylation at all, and genes with promoter hypomethylation and gene body methylation. Thus, the relationship between cytosine methylation and gene expression is still far from clear, and existing predicting methods achieve accuracies that are far below theoretical maximal accuracy associated with biological variability.

[0006] An orthogonal approach to gene expression prediction has been proposed in Avsec et al. (202 / ), where DNA sequence data (i.e. genetic information alone) was used to predict gene expression using a deep learning architecture termed “Enformer”. Enformer is trained to predict human and mouse genomic tracks at 128-bp resolution from 200 kb of input DNA sequence, where the genomic tracks include transcriptional activity measured by CAGE (Cap analysis of gene expression, histone modification (measured by histone modification ChIP sequencing), transcription factor binding (measured by transcription factor ChIP sequencing) and DNA accessibility (measured by DNase-seq or ATA-seq). Thus, the determinants of gene expression at the DNA level are still not fully understood.

[0007] The present invention has been devised in light of the above considerations. Summary of the Disclosure

[0008] The present inventors have postulated that the use of technologies that are able to provide more detailed epigenetic information including the presence of hydroxymethylcytosine as well as the presence of methylcytosine, and in particular technologies that are able to provide such epigenetic information in a single assay, could significantly improve on the prior art gene expression prediction methods. The inventors further postulated that additional benefits would be achievable if the single assay can also provide genetic information about the same DNA molecule (i.e. combining genetic and epigenetic information on a single read) for which epigenetic information is obtained. Such methods can be referred to as “6-letters sequencing” (or “6L-sequencing”). To test this hypothesis, they leveraged an assay (described in Fullgrabe et al. 2023) that enables sequencing the complete genetic sequence and the DNA modifications, 5- methylocytosine (5mC) and 5-hydroxy- methylocytosine (5hmC), from low nanogram amounts of DNA, to provide 6-base genomic data. They trained and evaluated a series of machine-learning models to predict different gene expression metrics from 6-base sequence data. The inventors demonstrated that the approach enabled the prediction of gene expression with higher accuracy than was previously possible. The approach is expected to enable prediction of gene expression with higher accuracy than was previously possible, as well as elucidate the complex interplay between genetic and epigenetic determinants of gene expression.

[0009] Thus, according to a first aspect, the disclosure provides a computer-implemented method of predicting the gene expression of one or more genes in a sample, the method comprising: receiving sequence data associated with the one or more genes obtained from the sample, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine; and predicting, for each gene of the one or more genes and using the sequence data associated with the gene, the value of a gene expression metric associated with the gene, wherein the predicting is performed using a machine learning model trained to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a gene expression metric for the gene.

[0010] The methods according to the present aspect may have any one or more of the following optional features.

[0011] The machine learning model may have been trained using a training dataset comprising: sequence data associated with a plurality of genes, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxy methylated cytosine; and expression data comprising for each of the plurality of genes, one or more measured values of a gene expression metric. The training dataset may comprise sequence data obtained from a first plurality of samples, and expression data obtained from a second plurality of samples. The first and second plurality of samples may be from the same organism, tissue and / or cell type or cell line. The machine learning model may have been trained using training sequence data and expression data from the same organism, tissue and / or cell type or cell line as that of the sample for which gene expression is predicted. The machine learning model may have been trained using a training dataset comprising expression data obtained by RNA sequencing, and the predicted gene expression metric may be a gene expression metric obtainable by RNA sequencing. The RNA sequencing may be an RNA sequencing technology indicative of steady state transcript levels, optionally bulk RNA sequencing. The RNA sequencing may be an RNA sequencing technology indicative of newly synthesised RNA over a predetermined period of time, optionally TT-seq. The gene expression metric may be a metric derived from read counts, optionally RPM (reads per million), RPKM (reads per kilo base per million mapped reads), or TPM (transcripts per million), or a metric derived therefrom by log transformation, such as a log(RPKM) or log(TPM).The machine learning model may be a regression model configured to output a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene. The machine learning model may be a regression model configured to output a value between predetermined bounds, optionally between 0 and 1 , and predicting the value of a gene expression metric associated with the gene further may comprise obtaining a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene using a predetermined function derived from the observed dynamic range of the one or more transcripts in a plurality of samples. The machine learning model may be a classification model, and predicting the value of a gene expression metric associated with the gene may comprise classifying the gene, using the machine learning model, between a plurality of classes associated with respective ranges of values of the gene expression metric, optionally wherein the plurality of classes and associated ranges correspond to respective quartiles of a distribution of values of the gene expression metric obtained from a plurality of samples.

[0012] In embodiments, the machine learning model takes as input a plurality of features derived from sequence data including at least one feature indicative of the presence of hydroxymethylated cytosine, and the method further comprises determining, using the sequence data, values for each of said plurality of features. The plurality of features may include, for each of a plurality of genomic regions associated with the genes, one or more of: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the presence of cytosines that are either methylated or hydroxymethylated in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region. The plurality of features may include: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region. The plurality of features may include the same features for all regions or a different set of features depending on the region, provided that at least one of the regions includes a feature indicative of the presence of hydroxymethylated cytosines in the region and all regions include at least one of: a feature indicative of the presence of hydroxy methylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, and a feature indicative of the presence of methylated or hydroxymethylated cytosines in the region. The epigenetic data may comprise sequence reads, and the feature indicative of the presence of hydroxy methylated cytosines in the region may be a summarised value, over one or more CpGs in the region, of the fraction of reads indicative of the presence of a hydroxymethylated cytosine at the CpG. The epigenetic data may comprise sequence reads, and the feature indicative of the presence of methylated cytosines in the region may be a summarised value, over one or more CpGs in the region, of the fraction of reads indicative of the presence of a methylated cytosine at the CpG. A summarised value may be a mean value.

[0013] In embodiments, the plurality of regions include one or more or all of: one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene, a region corresponding to a promoter like sequence, a region of a predetermined length centred on the transcription start site of the gene, a region corresponding to the 5’UTR of the gene, a region corresponding to the first exon of the gene, a region corresponding to the first intron of the gene, a region combining regions corresponding to all exons of the gene, a region combining regions corresponding to all introns of the gene, a region corresponding to the 3’UTR of the gene, and one or more regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene. In embodiments, the one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene comprise a plurality of non-overlapping regions of equal length between 100bp and 500bp, such as about 200bp, together covering a region up to 2kb, 1.8kb, 1.6kb, 1.5kb, 1 .4kb, 1.2kb or 1 kb from the transcription start site of the gene, and / or wherein the one or regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene comprise a plurality of nonoverlapping regions of equal length between 500bp and 1500bp, such as about 1 kb, together covering a region up to 5kb, 6kb, 7kb or 8kb from the 3’ end of 3’UTR of the gene.

[0014] The machine learning model may be a random forest model, a gradient boosted tree model, or a neural network model.

[0015] Alternatively, the machine learning model may be a deep learning model trained to take as input sequence data associated with a gene and produce as output a gene expression metric for the gene, wherein the sequence data further comprises genetic data associated with the gene. Thus, the machine learning model may be a deep learning model trained to: take as input encoded sequence data for a genomic sequence associated with a gene, the encoded sequence data including, for each position of the genomic sequence, encoded genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxy methylated cytosine; and produce as output a predicted the value of a gene expression metric associated with the gene. In embodiments, the encoded sequence data further comprises for each base in the sequence, information indicating whether the base is on the sense or antisense strand. In embodiments, the encoded sequence data comprises encoded genetic data and encoded epigenetic data for each position of each strand of the sense and antisense strand of the genomic sequence. In embodiments, the encoded sequence data includes, for each position of the sequence, genetic data encoded using a one-hot encoding scheme. In embodiments, the encoded sequence data comprises encoded epigenetic data that includes, for each position of the sequence that is a C: one or more of, optionally at least three of: a number of reads that have mC at the position, a number of reads that have hmC at the position, a number of reads that have C at the position, and a total number of reads at the position, each number of read obtained from a sequencing technology that distinguishes between C, mC and hmC. Instead or in addition to this, the encoded sequence data may comprise encoded epigenetic data that includes, for each position of the sequence that is a C: one or more of, optionally at least two of: a fraction of reads that have mC at the position, a fraction of reads that have hmC at the position, and a fraction of reads that have C at the position, each fraction of read obtained from counts of reads obtained using a sequencing technology that distinguishes between C, mC and hmC. Instead or in addition to this, the encoded sequence data may comprise encoded epigenetic data that includes, for each position of the sequence that is a C: one or more of, optionally all of: left and right boundaries of a confidence interval around an estimate of a fraction of reads that have mC at the position, left and right boundaries of a confidence interval around an estimate of a fraction of reads that have hmC at the position, and left and right boundaries of a confidence interval around an estimate of a fraction of reads that have C at the position, each confidence interval obtained from counts of reads obtained using a sequencing technology that distinguishes between C, mC and hmC. The confidence intervals may be estimated using a sampling method assuming that the fractions follow a Dirichlet distribution. Instead or in addition to this, the encoded sequence data may comprise encoded epigenetic data that includes, for each position of the sequence that is a C: one or more of, optionally all of: a statistical estimate of a fraction of reads that have mC at the position, a statistical estimate of a fraction of reads that have hmC at the position, and a statistical estimate of a fraction of reads that have C at the position, each statistical estimate obtained from counts of reads obtained using a sequencing technology that distinguishes between C, mC and hmC. The encoded epigenetic data may further comprise a standard deviation or variance estimate associated with each of said statistical estimates. The statistical estimates may be means of respective categories of a Dirichlet distribution.

[0016] The deep learning model may be a deep learning sequence model. The deep learning model may be a deep neural network, The machine learning model may be a deep neural network comprising one or more transformer layers and / or one or more convolutional layers. In embodiments, the deep neural network comprises a plurality of convolutional layers with one or more residual connections. In embodiments, the deep neural network comprises one or more convolutional layers each comprising a 1 D convolution followed by an activation function. In embodiments, the deep neural network comprises one or more convolutional layers each comprising a dilated convolution. In embodiments, the deep learning model is a transformer based model or a model comprising one or more convolution layers, one or more transformer layers and a classification or regression layer. In embodiments, the deep learning model takes as input a sequence derived from 6-letter sequence data associated with a gene and the deep learning model has been trained using training data comprising sequences derived from 6-letter sequence data associated with a plurality of genes and expression data associated with said plurality of genes. The input sequence may be a sequence derived from the 6-letter sequence data encoded using for each base of the sequence the following encoding scheme: A=[1 ,0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C=[0,0,0,f1 ,f2,f3], and N=[0, 0,0, 0,0,0], where f1 is the fraction of unmodified C reads in the 6-letter sequence data at the base, f2 is the fraction of methylated C reads in the 6-letter sequence data at the base, and f3 is the fraction of hydroxy methylated C in the 6-letter sequence data at the base. The deep learning model may take as input a sequence comprising data for a predetermined number of bases corresponding to the sequence of the gene, a sequence of a first predetermined length upstream of the gene and a sequence of a second predetermined length downstream of the gene. In embodiments, the machine learning model has been trained using training data comprising, for each of a plurality of training genes: (a) encoded sequence data including: for each position of a genomic sequence associated with a training gene, encoded genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxymethylated cytosine; and (b) measured values of the one or more gene expression metrics. The genomic sequence associated with a gene may comprise a region of a predetermined size comprising the transcription start site of the gene. For example, the region may comprise a region of a predetermined size either side (i.e. upstream and downstream) of the transcription start site of the gene. This may also described as a region of a predetermined length centred on the transcription start site of the gene. The region may comprise a region from a first predetermined length upstream of the transcription start site of the gene to a second predetermined length downstream of the transcription start site of the gene, downstream of the transcribed portion of the gene, or downstream of a region corresponding to the 3’UTR of the gene. The region may comprise a region from a first predetermined length upstream of the 5’URT (which may be a length of 0, i.e. starting at the 5’URT) of the gene to a second predetermined length downstream of the transcription start site of the gene, downstream of the transcribed portion of the gene, or downstream of a region corresponding to the 3’UTR of the gene. Thus, the region may comprise a region upstream of the 5’UTR, the 5’URT region, a region upstream of the TSS, a region from the TSS to the end of the transcribed portion of the gene, the 3’UTR region, a region downstream of the 3’UTR region, a region from the TSS to the end of the 3’UTR of the gene, and / or a region from a first predetermined length upstream of the TSS to a second predetermined length downstream of the 3’UTR.

[0017] In any embodiment, the machine learning model may have been trained using a training dataset comprising expression data comprising for each of the plurality of genes, a single value of a gene expression metric, the single value obtained by summarising one or more values corresponding to respective transcripts associated with the gene and / or respective samples, and wherein the summarised value is a mean value across a plurality of transcripts and / or samples.

[0018] According to a second aspect, there is provided a computer-implemented method of training a machine learning model for predicting the gene expression of one or more genes in a sample, the method comprising: receiving a training data set comprising (i) sequence data associated with a plurality of genes, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxy methylated cytosine, and (ii) expression data comprising for each of the plurality of genes, one or more measured values of a gene expression metric; and training a machine learning model, using said training data, to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a predicted gene expression metric for the gene.

[0019] The methods according to the present aspect may have any one or more of the features described in relation to the first aspect. Methods according to the first aspect can comprise training the machine learning model according to the second aspect.

[0020] In embodiments of any aspect described herein, the sequence data may have been obtained using a sequencing technology from which epigenetic and genetic bases can be called on the same read. In embodiments of any aspect described herein, the sequence data may comprise or consist of reads in 6- letters code. According to a third aspect, there is provided a method of performing a genetic test for a subject, the method comprising: obtaining sequence data associated with a sample previously obtained from a subject, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxy methylated cytosine; and performing the method of any embodiment of the first aspect to obtain a predicted gene expression metric for one or more genes. The sequence data may further comprise genetic data indicative of the presence of one or more genetic bases at one or more positions in the genome. The method may further comprise analysing the predicted gene expression metrics and / or genetic data by obtaining the value of one or more metrics characterising the sample. The one or more metrics may be selected from a diagnostic or prognostic metric.

[0021] The methods according to the present aspect may have any of the features described in relation to any preceding aspect.

[0022] According to a fourth aspect, there is provided a method of determining the effect of one or more perturbations on gene expression in a sample, the method comprising: obtaining sequence data associated with a sample previously exposed to the perturbation, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine; and performing the method of any preceding aspect to obtain a predicted gene expression metric for one or more genes. The method may further comprise analysing the predicted gene expression metrics by comparing the predicted gene expression metrics to respective control values. The method may further comprise exposing one or more samples to one or more perturbations. The one or more perturbations may be selected from: exposure to a drug, exposure to a physico-chemical stress, and exposure to a metabolic stress. The control values may be corresponding values obtained for one or more control samples.

[0023] The methods according to the present aspect may have any of the features described in relation to any preceding aspect.

[0024] According to a further aspect, there is provided a system comprising at least one processor and a non- transitory computer readable medium comprising instructions that, when executed by at least one processor, cause the at least one processor to perform the method of any embodiment of any of the preceding aspects, or any method described herein.

[0025] According to a further aspect, there is provided a non-transitory computer readable medium comprising instructions that, when executed by at least one processor, cause the at least one processor to perform the method of any embodiment of any of the first to sixth aspects, or any method described herein.

[0026] According to a further aspect, there is provided a computer program comprising code which, when the code is executed on a computer, causes the computer to perform the method of any embodiment of any of the first to sixth aspects, or any method described herein.

[0027] The invention includes the combination of the aspects and preferred features described except where such a combination is clearly impermissible or expressly avoided. Summary of the Figures

[0028] Embodiments and experiments illustrating the principles of the invention will now be discussed with reference to the accompanying figures in which:

[0029] Figure 1 is a flowchart illustrating a method of predicting gene expression according to a general embodiment of the disclosure (A), and methods of performing a genetic test and / or screening one or more perturbations according to general embodiments of the disclosure (B).

[0030] Figure 2 is a flowchart illustrating a method of providing a trained model for predicting gene expression, according to embodiments of the disclosure.

[0031] Figure 3 shows schematically a system for implementing methods of the disclosure.

[0032] Figure 4 shows the relationship between RPM (reads per million), RPKM (reads per kilo base per million mapped reads), and TPM (transcripts per million) (top) and the distribution of RPM, RPKM and TPM (bottom) in a dataset used in the examples of the disclosure.

[0033] Figure 5 shows the results of prediction of gene expression using a method as described herein, in terms of relationship and R2metric between predicted expression and observed expression from bulk RNA sequencing (RNA-seq).

[0034] Figure 6 shows the feature importance in terms of R2for the predictive features of the model used to generate the results in Figure 5, separating the contribution of methylated cytosine (referred to herein as 5mC or mC) information alone (light colour) and the added contribution of hydroxymethylcytosine (referred to herein as 5hmC or hmC) information (dark colour).

[0035] Figure 7 the results of prediction of gene expression using a method as described herein, in terms of relationship and R2metric between predicted expression and observed expression from transient transcriptome sequencing (TT-seq).

[0036] Figure 8 shows the feature importance in terms of R2for the predictive features of the model used to generate the results in Figure 7, separating the contribution of methylated cytosine information alone (light colour) and the added contribution of hydroxymethylcytosine (dark colour).

[0037] Figure 9 shows schematically alternative approaches for predicting gene expression from 6 letters sequencing data. TSS=transcription start site; 5’UTR=5’ untranslated region; PLS=promoter-like sequence; 3’UTR=3’ untranslated region.

[0038] Figure 10 shows examples of encoding schemes for encoding 4-base and 6-base sequence data.

[0039] Figure 10A shows an example of one-hot encoding of canonical 4-base data. An additional number can be used to indicate the strand (which can be e.g. 0 for the positive / sense strand and 1 for the negative / antisense strand) that the base is associated with.

[0040] Figure 10B shows examples of encoding of 6-base data which uses counts or fractional counts at C positions. As above, an additional number can be included to indicate the strand that the base is associated with. Figure 10C shows examples of encoding of 6-base data which uses statistics of model-based estimates of fractions of C, mC and 5hmC at C positions. sd=standard deviation. Cl=confidence interval.

[0041] Where the figures laid out herein illustrate embodiments of the present invention, these should not be construed as limiting to the scope of the invention. Where appropriate, like reference numerals will be used in different figures to relate to the same structural features of the illustrated embodiments.

[0042] Detailed Description

[0043] Aspects and embodiments of the present invention will now be discussed with reference to the accompanying figures. Further aspects and embodiments will be apparent to those skilled in the art. All documents mentioned in this text are incorporated herein by reference.

[0044] The methods described herein are computer-implemented unless context specifies otherwise (such as e.g. where measurement steps and / or wet steps are involved). Thus, the methods described herein are typically performed using a computer system or computer device. Any reference to an action such as “obtaining”, “processing”, “determining” may therefore refer to a processor performing the action, or a processor executing instructions that cause the processor to perform the action. Indeed, the methods of the present invention comprising sequence data analysis (requiring the analysis of thousands of reads per sample) and machine learning model training and / or deployment, is such that it cannot be performed in the human mind.

[0045] As used herein, the terms “computer system” of “computer device” includes the hardware, software and data storage devices for embodying a system or carrying out a computer implemented method. For example, a computer system may comprise one or more processing units such as a central processing unit (CPU) and / or a graphical processing unit (GPU), input means, output means and data storage, which may be embodied as one or more connected computing devices. Preferably the computer system has a display or comprises a computing device that has a display to provide a visual output display (for example in the design of the business process). The data storage may comprise RAM, disk drives or other computer readable media. The computer system may include a plurality of computing devices connected by a network and able to communicate with each other over that network. For example, a computer system may be implemented as a cloud computer. The term “computer readable media” includes, without limitation, any non-transitory medium or media which can be read and accessed directly by a computer or computer system. The media can include, but are not limited to, magnetic storage media such as floppy discs, hard disc storage media and magnetic tape; optical storage media such as optical discs or CD-ROMs; electrical storage media such as memory, including RAM, ROM and flash memory; and hybrids and combinations of the above such as magnetic / optical storage media. A computer system may comprise one or more networked computers, including e.g. a cloud computer.

[0046] The methods described herein may be provided as computer programs or as computer program products or computer readable media carrying a computer program which is arranged, when run on a computer, to perform the method(s) described herein. The present disclosure relates broadly to the use of epigenetic data (alone or in combination with genetic data - collectively, “sequence data”) and machine learning to predict one or more gene expression metrics.

[0047] DNA encodes information in both genetic and epigenetic bases. Epigenetic information is closely associated with transcription, ultimately influencing cell fate, disease etc. DNA methylation, the covalent addition of a methyl group to cytosine residues (primarily at CpG sites) is the most widely studied epigenetic mark. DNA methylation is typically measured using whole genome bisulfite sequencing (WGBS, see Frommer et al. 1992), enzymatic methyl-sequencing (EM-seq, see Vaisvila et al. 2021) or DNA methylation microarrays. All of these technologies investigate methylation only and genetic information is acquired separately, typically through high-throughput sequencing (also referred to as next generation sequencing) (Villicana and Bell, 2021). Recently, approaches for the simultaneous sequencing of genetic and epigenetic bases have been proposed. Fullgrabe et al. (2023) described a single base-resolution sequencing methodology that identifies G, C, T and A as well as 5-methylcytosine and 5-hydroxymethylcytosine (either combined, also known as “5-letters sequencing”, or separately, also known as “6-letters sequencing”). Sequence data refers to information indicative of the presence of any of the 4 genetic bases (A, T, C, G) and / or any epigenetic base (e.g. 5mC and / or 5hmC) at one or more positions in the genome. Sequence data encompasses epigenetic data and genetic data. Genetic data refers to information indicative of the presence of any of the 4 genetic bases (A, T, C, G) at one or more positions in the genome. Epigenetic data refers to information indicative of the presence of one or more epigenetic bases at one or more positions in the genome. Epigenetic bases include methylated cytosine (5-mC or 5mC), hydroxymethylated cytosine (5-hmC or 5hmC), 5-carboxycytosine (5-caC), formylcytosine (5-fC), and methyladenosine (6-mA). The present disclosure relates in particular to the use of epigenetic data including information about the presence of 5-methylcytosine and 5-hydroxymethylcytosine (also known as “6-letters sequencing”). Thus, in embodiments, epigenetic bases include or are selected from methylated cytosine (5-mC or 5mC) and hydroxy methylated cytosine (5-hmC or 5hmC). In embodiments, sequence data is acquired using a sequencing technology from which epigenetic and genetic bases can be called on the same read. The technology used may be a technology that identifies A, T, C, G and methylated cytosine, also referred to as “5-letters sequencing”. Alternatively, the technology used may be a technology that identifies A, T, C, G, mC and hmC, also referred to as “6-letters” sequencing. Sequencing technologies from which epigenetic and genetic bases can be called on the same read include the technology provided by biomodal (see biomodal.com / product / ), described in e.g. WO 2022 / 023753 A1 and Fullgrabe et al. (2023), Oxford Nanopore technologies (e.g. using the PromethlON instrument) as described in Simpson et al. (2017), or Pacific Biosciences (e.g. using 5-base HiFi sequencing as described in PacBio 2022, “Measuring DNA methylation with 5-base HFi sequencing”, available at www.pacb.com / wp-content / uploads / application- brief-measuring-dna-methylation-with-5-base-hifi-sequencing.pdf). Sequencing technologies from which epigenetic and genetic bases can be called on the same read in 6 letters code include the technology provided by biomodal (see biomodal.com / product / ), described in e.g. WO 2022 / 023753 A1 and Fullgrabe et al. (2023).

[0048] Gene expression metrics are metrics indicative of the level of expression of a gene, i.e. the amount of transcripts that are produced from a gene. Thus, a gene expression metric may also be referred to as a transcription metric or transcriptional metric. A gene expression metric may be in the form of a transcript level (i.e. level of expression of a transcript) or summarised transcript level across a plurality of transcripts associated with the same gene. A transcript level can be expressed using different metrics depending on the technology used to measure transcript level. For example, when using microarray-based technologies the metric may be an intensity or relative intensity. When using a sequencing-based technology, the metric may be a metric derived from the count of reads that map to a transcript. When using an amplificationbased method such as qRT-PCR the metric may be a cycle number. In embodiments, a gene expression quantifies a transcript level in units that can be obtained by RNA sequencing. These are typically units derived from read counts, such as RPM (reads per million), RPKM (reads per kilo base per million mapped reads), or TPM (transcripts per million), or metrics derived therefrom by log transformation (e.g. log2(RPKM)). RPM is typically defined as the number of reads that map to a given transcript or gene divided by the number of reads mapped to the reference sequence (reference genome or transcriptome), and then multiplied by 106. The RPKM metric is similar to the RPM, but normalised by the gene or transcript length. In practice, the gene / transcript length itself is typically normalised by 1 kbp, so that RPKM = RPM * 1000 / gene length or RPKM = RPM * 1000 / transcript length. TPM can be derived from RPKM as TPM = RPKM * 106 / sum(RPKM). Because they only differ by a normalisation factor, TPM and RPKM are linearly correlated, and can typically be used interchangeably in methods described herein (provided that the same unit is consistently used to train and use a machine learning model) without impacting accuracy. In embodiments, an expression metric is measured in one unit and used to train a model in a different unit. For example, expression may be measured in log2(RPM), which can be normalised by gene length (or transcript length, depending on whether the metric is provided at the gene level or transcript level), then Iog2(103) can be added to the result, effectively obtaining a log2(RPKM) value. The use of RPKM (or TPM) is advantageous over RPM as it removes a bias inherent of RPM towards longer genes having more transcripts simply by virtue of being longer, not necessarily because they are more expressed. In embodiments, a summarised expression metric is obtained across a plurality of transcripts for the same gene, and / or across a plurality of replicates, prior to using the expression metric to train a model as described herein. A summarised expression metric may be obtained as the mean, median or trimmed versions thereof across a plurality of transcripts and / or replicates.

[0049] Gene expression metrics can be measured (for use as ground truth or test data when training a model as described herein) or predicted. Gene expression metrics can be measured using any technology known in the art for determining the level of expression of a gene. In embodiments, measured gene expression metrics are obtained using RNA sequencing. References to RNA sequencing encompass any sequencingbased technology that can quantify the level of a transcript in a sample, including bulk RNA sequencing (sometime simply referred to as “RNAseq” or “RNA-seq”), single cell RNA sequencing, pseudo-bulk RNA sequencing (e.g. single cell RNA sequencing summarised at the level of a sample comprising a plurality of cells for which single cell sequencing data has been obtained), TT-sequencing or CAGE sequencing. TT-sequencing, also referred to as TT-seq (Schwalb et al. 2016) is a protocol that measures transcription for a limited period of time (e.g. 5 minutes) by labelling newly synthesized RNA through metabolic incorporation of 4-thiouridine (4sU) in live cells. Cap Analysis of Gene Expression (CAGE) is a method for promoter identification and transcription profiling (see e.g. Takahashi et al. 2012). The technique involves reverse transcription of RNA, capture of cDNAs through biotinylated RNA cap and 3’ ends, and sequencing of the captured sequences (CAGE tags). By counting the number of CAGE tags for each promoter within a gene, it is possible to determine the RNA expression level and the transcription start site of the transcript. Thus, the RNA sequencing technology used may be one which provides data indicative of steady state transcript levels, such as bulk RNA sequencing, single cell sequencing, CAGE sequencing, or pseudo-bulk RNA sequencing. Alternatively, the RNA sequencing may be an RNA sequencing technology indicative of newly synthesised RNA over a predetermined period of time, such as TT-seq. The present inventors have surprisingly discovered that sequence data and in particular epigenetic data can be used to obtain extremely accurate machine learning models predicting the levels of nascent transcripts.

[0050] The term “machine learning model” refers to a mathematical model that has been trained to predict one or more output values (predicted variables) based on input data (predictive variables, also referred to as values of predictive features). Training refers to the process of learning, using training data, the parameters of the mathematical model that result in a model that can predict outputs values that satisfy an optimality criterion or criteria. In the case of supervised learning, training typically refers to the process of learning, using training data, the parameters of the mathematical model that result in a model that can predict output values with minimal error compared to comparative (known) values associated with the training data (where these comparative values are commonly referred to as “labels” or “ground truth”). The are two major types of supervised learning models: classification models and regression models. Classification models aim to classify observations between a plurality of categories. They are typically suitable when the output to be predicted is a categorical label in the training data (e.g. class of expression level). Regression models aim to predict the value of a continuous variable associated with observations (e.g. a continuous expression level). They are typically suitable when the output to be predicted is a continuous value in the training data. A classification model may provide as output a classification label and / or one or more probabilities of an observation belonging to respective one or more classes. A binary classification model may provide as output a single probability indicating the probability that the observation belongs to a positive class (e.g. highly expressed genes). A predetermined threshold may be applied to the one or more probabilities to assign a class label. The threshold may be determined based on desired characteristics of the classification, such as e.g. a desired level of specificity, sensitivity or any accuracy performance that combines aspects of specificity (precision) and sensitivity (recall), such as accuracy and F1 score or balanced versions thereof that take into account the proportions of observations in the training data in each of the classes.

[0051] The term “machine learning algorithm” or “machine learning method” refers to an algorithm or method that trains and / or deploys a machine learning model. The machine learning models of the present disclosure are trained by supervised learning. An optimality criterion may be the minimisation of a loss function that quantifies the model prediction error based on the observed (ground truth) and predicted values of the predicted variables. Suitable loss functions for use in training machine learning models are known in the art and include the mean squared error, and the mean absolute error. Any of these can be used according to the present disclosure. The mean squared error (MSE) can be expressed as £(•) = MSE Xi) = (x;- x,)2where xtand xtare the predicted and observed (ground truth) values for a data instance (in the present case, a candidate drug combination in the training data), respectively. The mean absolute error (MAE) can be expressed as: £(•) = MAE(xi,xi') = |x;- x where xtand xtare the predicted and observed (ground truth) values for a data instance (in the present case, a candidate drug combination in the training data), respectively. The MAE is believed to be more robust to outlier observations than the MSE. The MAE may also be referred to as “L1 loss function”. However, MSE remains a very commonly used loss functions especially when a strong effect from outliers is not expected, as it can make optimization problems simpler to solve. Regularised loss functions are functions that include a loss function as described above, and one or more terms penalizing model complexity in order to reduce the risk of overfitting. Overfitting is a phenomenon that occurs where a machine learning model is trained to very closely reproduce the features of a training data set, resulting in poorer performance on other datasets that do not have the same characteristics (i.e. poor generalizability). L1 regularisation (also known as “Lasso” in the context of regression) add a regularization term to the loss function that penalizes models based on the sum of absolute value of the coefficients of the model. L2 regularisation (also known as “Ridge” in the context of regression) add a regularization term to the loss function that penalizes models based on the sum of squared value of the coefficients of the model. L1 regularisation can be used as a feature selection method as it minimizes the coefficients associated with less informative predictive features. In embodiments, the machine learning model is a regularized model. In embodiments, the machine learning model is a regularized tree-based model. In embodiments, the machine learning model is a regularized gradient boosted decision tree model. Examples of such models are available in the XGBoost software library (xgboost.ai / ). Such models may be referred to as “XGBoost” models, although any other implementation of regularized gradient boosted models may equally be used.

[0052] A machine learning model as described herein may be selected from: decision trees and variants thereof including regularised and / or gradient boosted decision trees and random forest models, regularised discriminant analysis, logistic regression models, artificial neural networks (ANNs) including multilayer perceptrons (with linear or non-linear activation functions) and deep learning models (e.g. long short-term memory networks (LSTMs), Recurrent neural networks (RNNs), Generative Adversarial Networks (GANs), etc), naive Bayes classifiers, support vector machines (SVM, using linear or non-linear kernels such as radial basis function) and multivariate adaptive regression splines (MARS). The present inventors have found it to be particularly beneficial to use non-linear models, such as decision trees and variants thereof (including in particular random forests and gradient boosted trees), SVM with a non-linear kernel, and ANNs (e.g. multilayer perceptrons) with non-linear activation functions.

[0053] In embodiments, a machine learning model comprises an ensemble of models whose predictions are combined. Alternatively, a machine learning model may comprise a single model. Random forest models and gradient boosted tree models (such as XGBoost) are ensemble models. Ensemble versions of any models can be constructed. Ensemble models are expected to result in better prediction performance than single models. For example, the machine learning model may be a random forest classifier or regressor, or a gradient boosted decision tree or regressor model. A random forest classifier is a model that comprises an ensemble of decision trees and outputs a class that is the average prediction of the individual trees. Decision trees perform recursive partitioning of a feature space until each leaf (final partition sets) is associated with a single value of the target. Gradient boosting is a machine learning method that forms an ensemble of weak prediction models (e.g. decision trees) from which a combined strong prediction is obtained. The algorithm iteratively adds new weak predictors to improve the prediction obtained by combining the outputs of the weak predictors. By contrast, random forest iteratively trains a set number of trees using random subsets of the training data. Ensemble models can also be constructed for deep learning sequence models, such as e.g. by training multiple models for the same task using different random seeds, and combining the outputs of such models, e.g. by averaging.

[0054] The term “deep learning sequence model” refers to a deep learning model that is configured to take as input encoded sequence data, and learn an internal representation (also called embedding) from which a property of interest (e.g. gene expression metric) can be predicted. A deep learning sequence model is typically a deep neural network. Multiple architectures can be used. For example, a transformer based deep learning model architecture can be used due to its computational efficiency in handling large input strings. Alternatively, a convolutional neural network can be used. Further, mixtures of convolutional layers and transformer layers can be used. These may be referred to as a “convformer” model. A transformer layer is a layer that processes input data using only attention mechanisms, as described in Vaswani et al. 2017. A convolutional layer is the main building block of a convolutional neural network, that applies a set of filters (also referred to as kernels) to input data. A model comprising convolutional layers may comprise a plurality of convolutional layers with one or more residual connections. A residual connection adds the input of a layer (or a downsampled version thereof) to the output of a subsequent layer. This may also be referred to as a shortcut connection. A model as described herein may comprise one or more convolutional layers, each comprising a 1 D convolution followed by an activation function (e.g. a ReLU activation function). . A model comprising convolution layers may comprise one or more dilated convolutions. Dilated convolutions are described in Yu & Koltun 2016. Dilated convolutions can aggregate multi-scale contextual information without losing resolution. This is particularly useful to enable base level predictions that take the context of the rest of the sequence into account. In embodiments, the sequence model does not comprise any pooling operations in or between any convolutional layer(s). This preserves the single-base- resolution of the models. A model comprising one or more transformer layers may comprise one or more transformer decoder layers. This may be referred to as a transformer decoder stack. Each transformer decoder layer may comprise a multi-head self-attention mechanism followed by a feedforward network with e.g. ReLU activation. A model comprising one or more transformer layer may use learnable positional encodings added to the input of the first transformer layer. The one or more transformer layers may use a causal attention mask to ensure that each position only attends to previous positions. The one or more transformer layers may maintain the sequence-resolution throughout the transformer processing.

[0055] A deep learning sequence model may further comprise a regression head. The regression head may take as input embeddings produced by transformer and / or convolutional layers and produce a numerical output for each base of the input sequence. Alternatively, when base level resolution is not desired, the regression head may take as input embeddings produced by transformer and / or convolutional layers and produce a numerical output for a plurality of bases of the input sequence (or the entire input sequence). A deep learning model may comprise a classification head. The classification head may take as input embeddings produced by transformer and / or convolutional layers and produce a classification output (e.g. classification label or probability for each of one or more classes) for each base of the input sequence. Alternatively, when base level resolution is not desired, the classification head may take as input embeddings produced by transformer and / or convolutional layers and produce a classification output for a plurality of bases of the input sequence (or the entire input sequence). A deep learning sequence model may comprise a plurality of regression and / or classification heads, each trained to predict a different gene expression metric. In embodiments, an entire deep learning sequence model is trained for prediction of a particular gene expression metric. In other embodiments, one or more layers of the model are trained using a first training data set (e.g. comprising data for a plurality of sample types, such as e.g. different cell lines, cell types, tissues, organisms, etc.) and then further trained or included as part of a model comprising additional trainable layers (also referred to as fine tuning) using a second data set (e.g. comprising data for a single sample type, such as e.g. a single organism, cell type, tissue, disease status, etc.) The second training data set may be a subset of the first training data set. The pretrained layers may be frozen (i.e. parameters may be set after the first round of training) or fine-tuned in the second round of training.

[0056] An encoded sequence (also referred to as “encoded sequence data”) refers to a numerical representation of a sequence, for providing as input to a machine learning model as described herein. For the purpose of machine learning models that are not deep learning sequence models, encoded sequence data may comprise values for a plurality of features derived from the sequence data. For the purpose of a deep learning sequence model, the encoded sequence data preserves base specific and positional information, i.e. it includes information about the base present at each position in a sequence (and evidence for the presence of mC and hmC when using 6-base sequencing data) as well as the order of these bases in the sequence. An encoded sequence can be obtained from sequence information comprising genetic sequence information (i.e. A, C, T, G, N, at each position in a sequence) and epigenetic information (i.e. information indicating the presence / absence or proportion of C, mC or hmC at each position in a sequence). The latter may be set to a default value for all positions where the genetic information does not indicate the presence of a C. Thus, an encoded sequence as described herein may comprise encoded genetic sequence information and encoded epigenetic information. Encoded genetic sequence information can be obtained from genetic sequence information using any method known in the art, such as e.g. one- hot encoding. An example one-hot encoding scheme for genetic sequence information is shown on Fig. 10A. Encoded epigenetic information can be obtained from mC / hmC sequencing information comprising, for each position comprising a C in a sequence, a number of reads that have mC at the position, a number of reads that have hmC at the position, and a number of reads that have C at the position (or two of these and a total number of reads, from which the third value can be obtained). Encoded epigenetic information can comprise one or more of these numbers of reads, or normalised versions thereof (i.e. fractions of reads that have mC, hmC or C at the position, respectively). Encoded epigenetic information can comprise one or more statistics associated with an estimate of the fraction of copies of the genomic region in the sample that have a C, mC and / or hmC, obtained from the mC / hmC sequencing information. For example, the statistics may be selected from: an estimate of the statistically most likely fraction of copies of the genomic region in the sample that have a C, mC and / or hmC, the boundaries of a confidence interval around an estimate of the statistically most likely fraction of copies of the genomic region in the sample that have a C, mC and / or hmC, and a mean and standard deviation (or variance) of an estimated distribution of the fraction of copies of the genomic region in the sample that have a C, mC and / or hmC. For example, this can be a posterior distribution for the respective fraction parameter. The statistics may be estimated using a model based on a Dirichlet process, which can capture uncertainty in multinomial count data. A confidence interval may be a confidence interval at a predetermined level of statistical confidence, such as e.g. 90%, 95%, 96%, or 97%, 98%. An encoded sequence may also comprise, for each base in the sequence, information indicating whether the base is on the sense or antisense strand. This may be encoded as a single bit of information (i.e. a single binary value) for each base in the sequence. An encoded sequence for a genomic region may comprise encoded sequence information for the sense strand and encoded sequence information for the antisense strand. Thus, the total length of the encoded sequence for a genomic region may be equal to 2*L*E, where L is the length of the genomic region, and E is the length of the encoding for a single base. For example, in the scheme illustrated on Fig. 10A, E is 4, for the schemes illustrated on Fig. 10B, E is 7, and for the schemes illustrated on Fig. 10C, E is 10. Each of these schemes can be supplemented with strand information for each base, leading respectively to schemes that have E=5 on Fig. 10A, E=8 on Fig. 10B, and E=11 on Fig. 10C.

[0057] Figure 1A is a flowchart illustrating a method of predicting gene expression according to a general embodiment of the disclosure.

[0058] At step 110, sequence data (also referred to herein as “sequence information”) is obtained for one or more samples. The term “sample” as used herein refers to a sample comprising genomic DNA or material from which genomic DNA can be obtained, such as e.g. cells, tissues etc, as well as samples comprising naturally fragmented genomic DNA such as cell free DNA (including circulating tumour DNA). For example, the sample may be a blood sample or sample derived therefrom (e.g. purified peripheral blood mononuclear cells, plasma), a cerebrospinal fluid sample or sample derived therefrom, a biopsy sample (e.g. a tissue biopsy sample or tumour sample), or a sample of cells (e.g. primary cells, which may be purified and / or cultured, cell lines, etc). The sample may be a sample from a subject, typically a eukaryote organism, such as a vertebrate, mammalian and / or a human subject. Similarly, the sample may be a sample of eukaryote cells, such as vertebrate or mammalian cells (e.g. mouse, rat, or human cells). The sample may be a sample from a human subject or model animal. The sample may be a sample of human or model animal cells.

[0059] The sequence data comprise epigenetic information (also referred to herein as “epigenetic data”) about the sample, i.e. information indicative of the presence of methylated cytosine and hydroxymethylated cytosine at a plurality of genomic positions. The sequence data may further comprise genetic information about the sample, i.e. information indicative of the presence of each of the 4 genetic bases at a plurality of genomic positions. The genetic and epigenetic information may have been obtained from the same assay, such as e.g. a sequencing assay in which genetic and epigenetic bases can be called on the same read. The sequence data may comprise a plurality of sequence reads (e.g. in the form of a SAM or BAM file per sample) or information derived therefrom (such as e.g. information identifying, for each read overlapping a genetic or epigenetic locus, whether the read contains a variant or reference allele / modified or unmodified base). According to any embodiment, the sequence data used for prediction and / or training may have been obtained using a sequencing technology from which epigenetic and genetic bases can be called on the same read. Thus, in any embodiment, the sequence data can comprise or consist of reads in 6-letters code. Reads in 6 letter code refers to reads in which each base is specified as A, T, C, G, mC or hMC (or N, indicating that the base could not be confidently called).

[0060] At step 120, encoded sequence data is obtained from the sequence data obtained at step 110. The encoded sequence data includes one or both of: (i) for each position of the sequence, encoded genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxymethylated cytosine, and (ii) values of one or more features derived from data including, for each position of the sequence, genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxymethylated cytosine, and the one or more features including at least one feature indicative of the presence of hydroxymethylated cytosine.

[0061] At step 130, the encoded sequence data is used to predict one or more gene expression metrics using the encoded sequence data. Thus step 130 can comprise predicting, for each gene of the one or more genes and using the sequence data associated with the gene, the value of a gene expression metric associated with the gene, wherein the predicting is performed using a machine learning model trained to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a gene expression metric for the gene. Multiple gene expression metrics may be predicted using one model, or a plurality of models may be used to predict respective gene expression metrics. Predictions may be made in a base specific manner (i.e. one prediction per base), for example using a deep learning sequence model as described herein, or at a gene level (for example using either a deep learning sequence model or a feature based machine learning model as described herein). The machine learning model may have been trained as described herein, for example by reference to Figure 2. The machine learning model may have been trained using training sequence data and expression data from the same organism, tissue and / or cell type or cell line as that of the sample for which gene expression is predicted. The gene expression metric may be a metric derived from read counts, such as RPM (reads per million), RPKM (reads per kilo base per million mapped reads), or TPM (transcripts per million), or a metric derived from one of these metrics by log transformation, optionally wherein the gene expression metric is a log(RPKM) or log(TPM). The machine learning model may be a regression model configured to output a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene. Alternatively, the machine learning model may be a regression model configured to output a value between predetermined bounds, such a between 0 and 1 . In such embodiments, predicting the value of a gene expression metric associated with the gene at step 130 may further comprise obtaining a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene using a predetermined function derived from the observed dynamic range of the one or more transcripts in a plurality of samples. For example, the value of the gene expression metric may be predicted based on a linear mapping between the predetermined bounds of the output value and the bounds of the observed dynamic range of the one or more transcripts in a plurality of samples. For example, an output value between 0 and 1 may be linearly mapped to the observed dynamic range of the one or more transcripts. In other embodiments, the machine learning model is a classification model, and predicting the value of a gene expression metric associated with the gene at step 130 comprises classifying the gene, using the machine learning model, between a plurality of classes associated with respective ranges of values of the gene expression metric. For example, the plurality of classes and associated ranges can correspond to respective quartiles of a distribution of values of the gene expression metric obtained from a plurality of samples. In embodiments of a method of predicting gene expression, the machine learning model takes as input a plurality of features derived from sequence data including at least one feature indicative of the presence of hydroxymethylated cytosine, and the method further comprises determining at step 120, using the sequence data, values for each of said plurality of features. The plurality of features may include, for each of a plurality of genomic regions associated with the genes, one or more of: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the presence of cytosines that are either methylated or hydroxymethylated in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region. For example, the plurality of features may include: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region. The plurality of features can include the same features for all regions or a different set of features depending on the region, provided that at least one of the regions includes a feature indicative of the presence of hydroxy methylated cytosines in the region and all regions include at least one of: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, and a feature indicative of the presence of methylated or hydroxymethylated cytosines in the region. As the skilled person understands, a feature indicative of the presence of hmC or a feature indicative of the presence of mC in a region can be derived from the combination of a feature indicative of the presence of mC or hmC (i.e. all modified Cs) in the region and one of a feature indicative of the presence of hmC and a feature indicative of the presence of mC in a region. Similarly, a feature indicative of the length of the region can in some embodiments be omitted and replaced by the use of features normalised by the length of the region. A feature indicative of the length of the region may simply be a length in bp. As the skilled person understands, when using regions of fixed length, a feature indicative of the length of the region may not be as informative and may be omitted. Thus, the plurality of features may include: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, and a feature indicative of the number of CpGs in the region. The plurality of features may not be the same for all regions, and for all gene expression metrics. For example, the inventors have found the presence of hydroxymethylated cytosines to be more informative in some regions than in others. Further, the inventors have found informative features to depend on whether the gene expression metric predicted is one that is indicative of steady-state transcript levels or nascent (new) transcript levels. For example, when predicting a gene expression metric that is indicative of steady state transcript levels, regions such as introns, 3’UTR and one or more regions downstream of the 3’UTR may include a feature indicative of the presence of hydroxymethylated cytosines in the region instead or in addition to a feature indicative of the presence of methylated cytosines in the region, whereas regions such as regions upstream of the TSS of the gene may not include a feature indicative of the presence of hydroxymethylated cytosines in the region, but may include a feature indicative of the presence of methylated cytosines in the region. As another examples, when predicting a gene expression metric that is indicative of nascent transcript levels (newly synthesised transcripts over a predetermined period of time), regions such as regions upstream of the TSS, and in particular those proximal to the TSS (e.g. within 1 kb of the TSS), regions around the TSS, promoter-like region, 5’URTR, first exon and optionally also first intron regions may include a feature indicative of the presence of hydroxymethylated cytosines in the region instead or in addition to a feature indicative of the presence of methylated cytosines in the region. Regions such as exons (combined), introns (combined), 3’UTR and regions downstream of the 3’UTR may not include a feature indicative of the presence of hydroxymethylated cytosines but would typically include a feature indicative fo the presence of methylated cytosines. The epigenetic data may comprise sequence reads. In such embodiments, a feature indicative of the presence of hydroxymethylated cytosines in the region can be selected from: (i) a summarised value (e.g. mean), over one or more CpGs in the region, of the fraction of reads indicative of the presence of a hydroxymethylated cytosine at the CpG; (ii) a density of hydroxymethylated cytosines in the region (this can be calculated as the mean hmC*count / length, where mean hmC is the mean fraction of reads indicative of the presence of a hydroxymethylated cytosine over the CpGs in the region, count is the number of CpGs in the region, and length is the length of the region (e.g. in base pairs); (iii) a standard deviation of the fraction of reads indicative of the presence of a hydroxymethylated cytosine at each CpG in the region; and (iv) an entropy of the hmC read counts over the region (this can be calculated as y k ; K-p_k (logi'“i(p_k))J where p_k is a summarised hydroxymethylation fraction calculated for a bin k of the region, where the region of divided into bins of equal length). Similarly, a feature indicative of the presence of methylated cytosines in the region can be selected from: (i) a summarised value (e.g. mean), over one or more CpGs in the region, of the fraction of reads indicative of the presence of a methylated cytosine at the CpG; (ii) a density of methylated cytosines in the region (this can be calculated as the mean mC*count / length, where mean mC is the mean fraction of reads indicative of the presence of a methylated cytosine over the CpGs in the region, count is the number of CpGs in the region, and length is the length of the region (e.g. in base pairs); (iii) a standard deviation of the fraction of reads indicative of the presence of a methylated cytosine at each CpG in the region; and (iv) an entropy of the mC read counts over the region (this can be calculated as -pfc(log (pfc)) where pk is a summarised methylation fraction calculated for a bin k of the region, where the region of divided into bins of equal length). Similarly, a feature indicative of the presence of modified cytosines in the region can be selected from: (i) a summarised value (e.g. mean), over one or more CpGs in the region, of the fraction of reads indicative of the presence of a modified cytosine at the CpG; (ii) a density of modified cytosines in the region (this can be calculated as the mean modC*count / length, where mean modC is the mean fraction of reads indicative of the presence of a modified cytosine over the CpGs in the region, count is the number of CpGs in the region, and length is the length of the region (e.g. in base pairs); (iii) a standard deviation of the fraction of reads indicative of the presence of a modified cytosine at each CpG in the region; and (iv) an entropy of the modC read counts over the region (this can be calculated as £fc-pfc(log (pfc)) where pk is a summarised modified C fraction calculated for a bin k of the region, where the region of divided into bins of equal length). The plurality of regions can include one or more or all of: one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene, a region corresponding to a promoter like sequence, a region of a predetermined length centred on the transcription start site of the gene, a region corresponding to the 5’UTR of the gene, a region corresponding to the first exon of the gene, a region corresponding to the first intron of the gene, a region combining regions corresponding to all exons of the gene, a region combining regions corresponding to all introns of the gene, a region corresponding to the 3’UTR of the gene, and one or more regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene. The one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene can comprise a plurality of nonoverlapping regions of equal length between 100bp and 500bp, such as about 200bp, together covering a region up to 2kb, 1.8kb, 1 .6kb, 1.5kb, 1.4kb, 1 .2kb or 1 kb from the transcription start site of the gene. The one or regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene can comprise a plurality of non-overlapping regions of equal length between 500bp and 1500bp, such as about 1 kb, together covering a region up to 5kb, 6kb, 7kb or 8kb from the 3’ end of 3’UTR of the gene. The one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene may be non-overlapping or partially overlapping. They may together cover the whole region between the TTS and a predetermined distance upstream from the TSS. The predetermined distance may be 2kb, 1 .8kb, 1 .6kb, 1 .5kb, 1 .4kb, 1 .2kb or 1 kb. The one or more regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene may be non-overlapping or partially overlapping. They may together cover the whole region between the 3’ end of 3’UTR of the gene and a predetermined distance downstream from the 3’ UTR. The predetermined distance may be 3kb, 4kb, 5kb, 6kb, 7kb or 8kb. Advantageously the predetermined distance may be at least 5kb or about 5kb. Overlapping regions may overlap by a stride of half the size of the region.

[0062] The regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene may be referred to as “promoter” regions. Promoters are typically considered to extend up to about 1 kb from the transcription start site. The present inventors have identified useful signal further than this distance and in particular up to about 1 .4 to 1 .6kb. The exact distance may depend on the length of the bins that is chosen when the bins are of equal length as the total distance is in such cases a multiple of the length of the non-overlap bins. The present inventors have found bins of approximately 200bp to be useful by examining the contribution of bins to the R2of the resulting models, showing and variable information content in such bins. In embodiments, determining, using the sequence data, values for each of said plurality of features comprises imputing a missing value for a feature using the values for other genes in a dataset (E.g. the training dataset). For example, where a gene does not have a region for which features are determined, the mean value of the respective features in that region in all other genes in the training data set may be used. A region corresponding to a promoter like sequence can be defined as a region annotated as a candidate cis-regulatory element (cCRE) that is a “promoter-like sequence” (PLS) in a database such as ENCODE. A region corresponding to a promoter like sequence can be defined as a region that has high DNase and H3K4me3 signal and that lies within 200 bp of an annotated transcription start site (TSS) (e.g. an annotated GENCODE TSS). In embodiments, the features include any one or more or all of the features listed in Table 1 . In embodiments, the machine learning model is a random forest model, a gradient boosted tree model, or a neural network model. Such architectures are particularly well suited for embodiments that use one or more features derived from the sequence data to make predictions. In other embodiments, the machine learning model is a deep learning model trained to take as input sequence data associated with a gene and produce as output a gene expression metric for the gene, wherein the sequence data further comprises genetic data associated with the gene. The deep learning model may be a transformer based model. The deep learning model may be or a model comprising one or more convolution layers, one or more transformer layers and a classification or regression layer. The deep learning model may take as input a sequence derived from 6-letter sequence data associated with a gene and the deep learning model may have been trained using training data comprising sequences derived from 6-letter sequence data associated with a plurality of genes and expression data associated with said plurality of genes. The input sequence may be a sequence derived from the 6-letter sequence data encoded using for each base of the sequence the following encoding scheme: A=[1 ,0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C=[0,0,0,f1 ,f2,f3], and N=[0, 0,0, 0,0,0], where f1 is the fraction of unmodified C reads in the 6-letter sequence data at the base, f2 is the fraction of methylated C reads in the 6-letter sequence data at the base, and f3 is the fraction of hydroxymethylated C in the 6-letter sequence data at the base. The deep learning model may take as input a sequence comprising data for a predetermined number of bases corresponding to the sequence of the gene, a sequence of a first predetermined length upstream of the gene and a sequence of a second predetermined length downstream of the gene.

[0063] At step 140, the results of any one or more of the preceding steps (including e.g. a predicted expression metric) can optionally be provided to a user (e.g. through a user interface) or data store.

[0064] The methods described herein find application in a variety of contexts. For example, the methods described herein can be used to predict gene expression in a genetic assay. A genetic assay is an assay that determines genetic information about a subject. Genetic assays are typically used to inform the likelihood of a subject developing a disease, to characterise a subject’s disease, to determine a prognostic or treatment response for a subject, etc. In embodiments, a genetic assay is an assay that selectively determines the presence of genetic variants at a one or more predetermined loci in a sample from a subject (e.g. mutations in one or more predetermined genes). This may also be referred to as a targeted assay, gene panel assay or simply panel assay. The genetic variants may be somatic variants or germline variants. The predetermined loci may be selected as loci that are relevant to the likelihood of developing a disease, to the prognostic associated with a disease or disorder, to the likelihood of response to a particular therapy, or to the characterisation of a disease (e.g. to characterise the subject as having a disease of a particular subtype or severity). Thus, the genetic assay may be used to obtain a polygenic risk score for a particular disease or disorder, based on the genotype of a subject at each of the predetermined loci in the genetic assay. A polygenic risk score is a score that indicates the risk of developing a particular disease or disorder (which can refer to a disease or a particular subtype or severity of a disease) based on the genotype of a subject at a plurality of loci. In other words, a polygenic risk score is a score that is calculated based on the determined genotype of a subject at a set of predetermined genetic loci, and which is indicative of the probability or odds ratio of a subject having a particular phenotype (such as e.g. developing a disease, acquiring resistance to a therapy, having a poor prognosis, etc.). The disease or disorder may be referred to as “complex disease”, by contrast with single-gene diseases which can be associated with a single gene (i.e. having a genetic variant at the particular gene causes the disease). Complex diseases are influenced by multiple genes and environmental factors such that the genotype at multiple genes, optionally in combination with one or more additional factors such as environmental factors, enables the prediction of a likelihood of developing the disease. Gene expression-based risk scores are also commonly obtained but typically require separate acquisition of gene expression information. This can also be performed in a targeted manner, using e.g. RT-PCR or Nanostring technologies. According to the present disclosure, it is possible to accurately predict gene expression from epigenetic information, it is possible to obtain insights obtainable by both a genetic assay (e.g. a polygenic risk score) and insights obtainable by an expression assay (or even new scores combining both genetic and expression related effects) using a single source of sequence data that contains genetic and epigenetic information. In other words, according to the present disclosure, a genetic assay can be used to provide information about gene expression.

[0065] Further, the methods of the present disclosure find use in the screening of one or more perturbations for an effect on one or more of: gene expression, and cytosine methylation / hydroxymethylation. Indeed, the methods described herein provide an opportunity to characterise perturbations in relation to their effect on the genome, the epigenome and the transcriptome using data from a single assay.

[0066] Indeed, the improved prediction accuracy provided by the methods described herein are particularly significant in the context of clinical assays and perturbation (e.g. drug) screening. In the context of clinical assay, any improvement in the accuracy of a prediction has meaningful implications in terms of the likelihood of patients being assigned to the correct clinical pathways and / or the likelihood of patients needing to undergo additional tests to confirm a clinical insight of the assay. In the context of drug screening, improved prediction accuracy (and especially in the context of the ability to obtain multiple types of information characterising the state of the perturbed system from a single assay) enables screening of larger sets of perturbations with lower risk of missing relevant candidates or spending additional resources on irrelevant candidates. Even small difference in accuracy can make large practical differences at when screening is performed at scale.

[0067] Additionally, the methods described herein are particularly useful in the context of analysing samples from which transcriptomic material cannot easily be obtained or where sample availability is limited, such as e.g. cell free DNA samples. Indeed, cfDNA samples such as those obtained from liquid biopsies are an increasingly common and convenient tool to diagnose, characterise and / or monitor diseases such as cancer. The present disclosure provide methods that augment the possibilities associated with analysis of such samples by enabling the quantification of expression metrics without the need to perform additional assays that may either not be possible to perform on this type of material or for which sufficient amounts of material is not available.

[0068] Figure 1 B is a flowchart illustrating methods of performing a genetic test and / or screening one or more perturbations according to general embodiments of the disclosure. At optional step 1100, sequence data is obtained from a sample, the sequence data comprising genetic and epigenetic information, including 5mC and optionally 5hmC information. Step 1100 may comprise analysing the sample using a sequencing protocol that identifies genetic and epigenetic bases (including mC and hmC) on the same reads. At step 1200, a method as described in relation to Fig. 1A is performed using the sequence data obtained at step 1 100. This may use sequence data that has been acquired using a genome-wide assay or a targeted assay that obtains sequence data about a selected set of genes. The selected set of genes may be disease associated genes. For example, in the context of a genetic test for cancer diagnosis or prognosis, genes (or specific genetic loci within genes) that are known to be associated with the risk of a subject developing cancer, the severity and / or drug response of a subject with cancer (e.g. mutations in tumour suppressor or tumour driver genes, mutations known to be associated with drug resistance or sensitivity, etc.) may be selected. At step 1300, the predictions obtained at step 1200, and optionally the sequence data obtained at step 1100 are analysed to obtain one or more metrics characterizing the sample. Obtaining sequence data may comprise sequencing material in a sample or receiving previously acquired sequence data from a memory, database or user interface.

[0069] In embodiments where the method is used to perform a genetic test, the sample may be a sample from a subject who has been diagnosed as having or being likely to have a disease or disorder. In such embodiments, the one or more metrics may include e.g. metrics indicative of a prognostic, treatment response, or diagnosis (including identification of a disease subtype such as a molecular subtype). For example, the genetic information (e.g. presence of one or more mutations at respective genetic loci) obtained at step 1100 may be used step 1300 to obtain a polygenic risk score. Defining a polygenic risk score for a particular phenotype of interest using a set of genetic loci has been identified as associated with this phenotype is within the capability of the skilled person. The epigenetic information (optionally in combination with the genetic information) is used at step 1200, to obtain one or more predicted gene expression metrics. The predicted gene expression may then be used at step z to obtain an expression based metric indicative of a prognostic, treatment response, or diagnosis (including identification of a disease subtype such as a molecular subtype). Alternatively, a multiomic metric may be obtained which combines information about the genotype of the subject as specific loci and the predicted level of expression of one or more genes.

[0070] In embodiments where the method is used to determine the effect of a perturbation (such as e.g. for drug screening, genetic KO screening, exposure to one or more physico-chemical or metabolic stresses etc.), the sample may be a sample of cells that have been exposed to the perturbation or a sample previously obtained from a subject that has been exposed to the perturbation. The one or more metrics may include metrics that compare the predictions obtained at step y to predictions obtained for a suitable control, in order to characterize the effect of the perturbation on gene expression and / or the state of epigenetic loci. Methods according to these embodiments may comprise performing the method for a plurality of samples exposed to respective perturbations, and selecting one or more of the plurality of perturbations for further screening based on the values of the one or more metrics obtained at step 1300. Further, methods according to this embodiment may further comprise exposing the one or more samples to the one or more perturbations. At optional step 1400, the results of any one or more of the preceding steps (including e.g. a predicted expression metric obtained at step 1200, metric characterising the sample obtained at step 1300, genetic and / or epigenetic information obtained at step 1100) are provided to a user (e.g. through a user interface) or data store.

[0071] Figure 2 is a flowchart illustrating a method of providing a trained model for predicting gene expression, according to embodiments of the disclosure.

[0072] In embodiments comprising training a machine learning model for predicting the gene expression of one or more genes in a sample, the method comprises at step 210 receiving a training data set comprising (i) sequence data associated with a plurality of genes, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine, and (ii) expression data comprising for each of the plurality of genes, one or more measured values of a gene expression metric. The training dataset may comprise data obtained from a plurality of samples. The training data set may comprise sequence data obtained from a first plurality of samples, and expression data obtained from a second plurality of samples. In other words, the samples from which sequence data and expression data was obtained may be different, although they are advantageously from the same organism, tissue and / or cell type or cell line (i.e. same type of sample). They may also have been subject to the same treatment (e.g. they may all be control samples that have not been perturbed, or they may all be samples that have been perturbed in the same way). When different types of samples or perturbations are used, it may be advantageous for samples of the same type and / or exposed to the same perturbation to be treated as forming pairs of ground truth data for the purpose of training. For example, expression values obtained in one cell type are used as ground truth for predictions from sequence data from the same cell type, whereas expression values obtained in another cell type are used as ground truth for predictions from sequence data from the same other cell type. The machine learning model may be trained using training data comprising sequence data and expression data from the same organism, tissue and / or cell type or cell line as that of a sample for which gene expression is to be predicted (i.e, intended use of the machine learning model). For example, the machine learning model may be trained using human data and deployed (used for predictions) on human data. The plurality of genes may comprise at last 1000 genes, at least 2000 genes, at least 3000 genes, at least 5000 genes, at least 6000 genes, at least 7000 genes, at least 8000 genes, at least 9000 genes or at least 10000 genes. Thus, the training data set may comprise sequence data associated with a plurality of genes, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine; and expression data comprising for each of the plurality of genes, one or more measured values of a gene expression metric.

[0073] The machine learning model may be trained using a training dataset comprising expression data comprising for each of the plurality of genes, a single value of a gene expression metric, the single value obtained by summarising one or more values corresponding to respective transcripts associated with the gene and / or respective samples. Thus, the model may comprise obtaining, for each gene of a plurality of genes represented in the training data, a single value for a gene expression metric, for example as a mean value across a plurality of transcripts and / or samples for the gene. Further, the machine learning model may be trained using a training dataset comprising sequence data comprising, for each of the plurality of genes, sequence data for a single transcript of one or more transcripts associated with the gene (e.g. the longest transcript). In other words, the model may be trained to predict gene expression at the gene level using sequence data at the gene level. The gene expression at the gene level may be obtained as a summary of multiple metrics for respective transcripts or samples, and the sequence data at the gene level may be obtained as the sequence data that maps to one of the transcripts of the gene. The definition of the chosen transcript can be used to specify regions to obtain values of features defined from the sequence data.

[0074] The machine learning model may be trained using a training dataset comprising expression data obtained by RNA sequencing, and the predicted gene expression metric may therefore be a gene expression metric obtainable by RNA sequencing. The RNA sequencing may be an RNA sequencing technology indicative of steady state transcript levels, such as bulk RNA sequencing. Alternatively, the RNA sequencing may be an RNA sequencing technology indicative of newly synthesised RNA over a predetermined period of time, such as TT-seq. Therefore, the gene expression metric can be a metric derived from read counts, such as RPM (reads per million), RPKM (reads per kilo base per million mapped reads), or TPM (transcripts per million), or a metric derived therefrom by log transformation. In embodiments, the gene expression metric is a RPKM, TPM, log(RPKM) or log(TPM). The method may comprise processing the training gene expression data to obtain gene expression metrics (for use as ground truth) in the chosen units, for example transforming RPM data to RPKM data. The machine learning model may be a regression model configured to output a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene, in which case the training gene expression data can be used directly as a ground truth for training. The machine learning model may be a regression model configured to output a value between predetermined bounds, such as between 0 and 1. In such embodiments, predicting the value of a gene expression metric associated with the gene further comprises obtaining a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene using a predetermined function derived from the observed dynamic range of the one or more transcripts in a plurality of samples. Thus, the training can further comprise obtaining, from the training gene expression data, values between predetermined bounds (e.g. between 0 and 1) for each gene or transcript by normalising the data between the bounds of the observed dynamic range of the gene or transcripts (in which case a linear mapping can be used) or by fitting any predetermined function that maps the observed values in the training gene expression data to values between the predetermined bounds. The reverse of this function can then be used when making predictions, to obtain absolute gene expression values for a gene. In embodiments, the machine learning model is a classification model, and predicting the value of a gene expression metric associated with the gene comprises classifying the gene, using the machine learning model, between a plurality of classes associated with respective ranges of values of the gene expression metric. In such embodiments, the methods can comprise defining a plurality of classes of gene expression (e.g. classes associated with respective ranges of gene expression metric values) and assigning ground truth class labels to the genes or transcripts in the training data (e.g. depending on which range the value for a gene or transcript falls in). For example, the plurality of classes and associated ranges can correspond to respective quartiles of a distribution of values of the gene expression metric obtained from a plurality of samples, such as quartiles of a distribution ofvalues ofthe gene expression metric over the training dataset.

[0075] At step 220, the training sequence data is encoded as described herein. The encoded training sequence data includes one or both of: (i) for each position of the sequence, encoded genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxymethylated cytosine, and (ii) values of one or more features derived from data including, for each position of the sequence, genetic data and encoded epigenetic data indicative of the presence of one or more epigenetic bases at the position, the one or more epigenetic bases including methylated cytosine and hydroxy methylated cytosine, and the one or more features including at least one feature indicative of the presence of hydroxymethylated cytosine. As explained above, the machine learning model may take as input a plurality of features derived from sequence data including at least one feature indicative of the presence of hydroxymethylated cytosine. In such embodiments, the method may further comprise determining at step 220, using the training sequence data, values for each of said plurality of features, the plurality of features may be as described above, for example in relation to Figure 1 A. As explained above, the machine learning model may be a random forest model, a gradient boosted tree model, or a neural network model. Each of these can be implemented as a regressor or a classifier, and is particularly well suited to making predictions from features rather than sequences (strings). In other embodiments, the machine learning model is a deep learning model trained to take as input sequence data associated with a gene and produce as output a gene expression metric for the gene, wherein the sequence data further comprises genetic data associated with the gene. The deep learning model may be a transformer based model or a model comprising one or more convolution layers, one or more transformer layers and a classification or regression layer. The deep learning model may take as input a sequence derived from 6-letter sequence data associated with a gene and training the deep learning model may comprise obtaining from the training data, sequences derived from 6-letter sequence data using a predetermined encoding scheme. For example, input sequences provided to the model (either at training time or at deployment time) may be a sequence derived from the 6-letter sequence data encoded using for each base of the sequence the following encoding scheme: A=[1 ,0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C=[0,0,0,f1 ,f2,f3], and N=[0, 0,0, 0,0,0], where f1 is the fraction of unmodified C reads in the 6-letter sequence data at the base, f2 is the fraction of methylated C reads in the 6-letter sequence data at the base, and f3 is the fraction of hydroxymethylated C in the 6-letter sequence data at the base. The deep learning model may take as input a sequence comprising data for a predetermined number of bases corresponding to the sequence of the gene, a sequence of a first predetermined length upstream of the gene and a sequence of a second predetermined length downstream of the gene. The first and second predetermined lengths can be optimised when training the machine learning model, for example by choosing lengths that provide the best prediction accuracy.

[0076] In embodiments comprising training a machine learning model for predicting the gene expression of one or more genes in a sample, the method further comprises at step 230 training a machine learning model, using said training data, to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a predicted gene expression metric for the gene. The machine learning model, its inputs and outputs and the training data may have any of the features described in relation to Figure 1A.

[0077] At step 240, the results of any preceding step may be provided to a user, computing device or data store. These may include one or more of: one or more trained models, one or more parameters of a trained model, one or more training and / or performance statistics that may be obtained as part of the train step, etc.

[0078] The training dataset may comprise data obtained from a plurality of samples. The training data set may comprise sequence data obtained from a first plurality of samples, and gene expression data obtained from a second plurality of samples. In other words, the samples from which sequence data and gene expression data was obtained may be different, although they are advantageously from the same organism, tissue and / or cell type or cell line (i.e. same type of sample). They may also have been subject to the same treatment (e.g. they may all be control samples that have not been perturbed, or they may all be samples that have been perturbed in the same way). When different types of samples or perturbations are used, it may be advantageous for samples of the same type and / or exposed to the same perturbation to be treated as forming pairs of ground truth data for the purpose of training. For example, values of gene expression metrics obtained in one cell type may be used as ground truth for predictions from sequence data from the same cell type, whereas values obtained in another cell type may be used as ground truth for predictions from sequence data from the same other cell type. The machine learning model may be trained using training data comprising sequence data and gene expression data from the same organism, tissue and / or cell type or cell line as that of a sample for which gene expression metrics are to be predicted (i.e, intended use of the machine learning model). For example, the machine learning model may be trained using human data and deployed (used for predictions) on human data. The plurality of genomic regions / genes may comprise at last 1000 regions, at least 2000 regions, at least 3000 regions, at least 5000 regions, at least 6000 regions, at least 7000 regions, at least 8000 regions, at least 9000 regions or at least 10000 regions. As explained above, the machine learning model may be a random forest model, a gradient boosted tree model, or a neural network model. Each of these can be implemented as a regressor or a classifier, and is particularly well suited to making predictions from features rather than sequences (strings). In other embodiments, the machine learning model is a deep learning model trained to take as input sequence data and produce as output position specific (also referred to herein as base specific) values or gene specific (also referred to herein as gene level) values. Position specific values may be e.g. a count of reads mapped to the position. Gene level values may be transcript level values or summarised values across a plurality of transcripts for the gene.

[0079] Figure 3 shows an embodiment of a system for implementing methods of the disclosure, such as e.g. providing a trained model, using a trained model to predict gene expression, performing a genetic test and / or screening one or more perturbations. The system comprises a computing device 1 , which comprises a processor 101 and computer readable memory 102. In the embodiment shown, the computing device 1 also comprises a user interface 103, which is illustrated as a screen but may include any other means of conveying information to a user such as e.g. through audible or visual signals. The computing device 1 is communicably connected, such as e.g. through a network, to one or more databases 2 storing sequence data from a sample or a cohort of samples, and / or to one or more sequence data acquisition means 4. In embodiments, the sequence data acquisition means 4 is configured to obtain sequence data from samples, in the form of DNA sequencing reads or RNA sequencing reads (where RNA sequencing reads may be obtained by sequencing of cDNA derived from RNA). Thus, the sequence data acquisition means may comprise a sequencing machine, such as e.g. an Illumina sequencer when using the duet multiomics solution from biomodal (such as e.g. the duet multiomics solution +modC which performs 5-letter sequencing or the duet multiomics solution evoC which performs 6-letter sequencing), the PromethlON sequencer from Oxford Nanopore Technologies, or the Revio or Sequel long-read sequencers from Pacific Biosciences. The samples may be samples that have been processed to enable the sequencing of reads in 5 or 6 letter code (i.e. A, C, T, G, modified C and / or methylated C and hydroxy methylated C). The sequence data acquisition means 4 may be connected to the database 2. Connection between the sequence data acquisition means 4 and the database 2 may be through a wired or wireless connection. The sequence data may be raw data (e.g. reads) and / or processed versions thereof, such as e.g. variant calling files, counts files, etc. The one or more databases 2 may further store one or more of: one or more reference sequences, one or more gene expression data sets, parameters (such as e.g. parameters of a trained, parameters of a sequence data preprocessing methods, etc.), etc. The computing device may be a smartphone, tablet, personal computer or other computing device. The computing device is configured to implement a method as described herein. In alternative embodiments, the computing device 1 is configured to communicate with a remote computing device (not shown), which is itself configured to implement a method as described herein. In such cases, the remote computing device may also be configured to send the result of the method to the computing device. Further, the various steps of the methods described herein may be split between the computing device 1 and the remote computing device. The remote computing device may be a cloud computing device, a server node, etc. Any processing device known in the art may be used for this purpose. Communication between the computing device 1 and the remote computing device may be through a wired or wireless connection, and may occur over a local or public network 3 such as e.g. over the public internet.

[0080] The features disclosed in the foregoing description, or in the following claims, or in the accompanying drawings, expressed in their specific forms or in terms of a means for performing the disclosed function, or a method or process for obtaining the disclosed results, as appropriate, may, separately, or in any combination of such features, be utilised for realising the invention in diverse forms thereof.

[0081] While the invention has been described in conjunction with the exemplary embodiments described above, many equivalent modifications and variations will be apparent to those skilled in the art when given this disclosure. Accordingly, the exemplary embodiments of the invention set forth above are considered to be illustrative and not limiting. Various changes to the described embodiments may be made without departing from the spirit and scope of the invention.

[0082] For the avoidance of any doubt, any theoretical explanations provided herein are provided for the purposes of improving the understanding of a reader. The inventors do not wish to be bound by any of these theoretical explanations.

[0083] Any section headings used herein are for organizational purposes only and are not to be construed as limiting the subject matter described. Throughout this specification, including the claims which follow, unless the context requires otherwise, the word “comprise” and “include”, and variations such as “comprises”, “comprising”, and “including” will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.

[0084] It must be noted that, as used in the specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the context clearly dictates otherwise. Ranges may be expressed herein as from “about” one particular value, and / or to “about” another particular value. When such a range is expressed, another embodiment includes from the one particular value and / or to the other particular value. Similarly, when values are expressed as approximations, by the use of the antecedent “about,” it will be understood that the particular value forms another embodiment. The term “about” in relation to a numerical value is optional and means for example + / - 10%.

[0085] Examples

[0086] EXAMPLE 1 - Prediction of RNA transcription

[0087] Methylation patterns have been shown to be associated with gene expression. At the highest level, methylation is thought to enhance gene expression. The relationship with hydroxymethylation is less well studied but hydroxymethylation has been postulated to suppress gene expression. However, the relation is more nuanced, and often depends on specific patterns of (hydroxy)methylation at various genomic regions. In this example, the inventors sought to establish a list of features at various genomic regions, based on methylated CpGs assessed using a 6-letter sequencing method, that are deemed relevant for gene expression, and use them to predict expression levels from RNA sequencing data. The inventors found that the use of 6-letters information did improve prediction of gene expression, uncovered major determinants of gene expression involving the complex interplay of methylation and hydroxymethylation in different regions in and around a gene, and demonstrated that the 6-letters information is able to predict dynamic transcription patterns with surprisingly high accuracy.

[0088] Methods

[0089] Data - Bulk RNAseq. Bulk RNAseq data was obtained from the Gene Expression Omnibus database, under accession GSE135509 (www.ncbi.nlm. nih.gov / geo / query / acc.cgi?acc=GSE135509). This contains bulk RNA-seq data for E14 and SAM mouse embryonic stem cells. SAM mESCs are a clonal cell line generated from E14s by lentiviral integration of dCas9-VP64 and MS2-p65-HSF1 . The data was originally published in Alda-Catalinas et al. 2021. The dataset contains data for 5 samples, 3 E14 replicates and 2 SAM replicates. Only data for the three replicates of an ES-E14 cell line were used here. The dataset consists of multiple rows (one per gene). For each gene, the columns are chromosome, start, end, and name of gene, as well as 3 additional columns: one expression level per replicate for each of the three replicates. The expression is given in log2(RPM) where RPM are Reads per million mapped reads (see definition here and discussion below). We create an additional column for the mean gene expression, as the mean of the 3 replicates. We further normalise this column by the gene length and add Iog2(103), effectively giving us an RPKM value. This is the column we will use to train our model. We favour RPKM over RPM as it removes a bias inherent of RPM towards longer genes having more transcripts simply by virtue of being longer, not necessarily because they are more expressed.

[0090] Data - TT-seq. TT-seq data was obtained from the Gene Expression Omnibus (GEO) database under accession GSE183278 (www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE183278). This contains TT-seq data for 4 replicates of E14 cells mouse embryonic stem cells. The TT-seq data is in the form of a bigwig file which we processed to produce a bedgraph file containing expression data. For each gene, we take the median expression of all the transcripts, giving us one single TT-seq expression level per gene. Similarly to the bulk RNA-seq dataset, we further normalise the expression value by the gene length and add Iog2(103) to create an RPKM value. Finally, we also normalise the expression by the mean gene length across the genes in the dataset.

[0091] Data - methylation. The 6-letters sequencing data used is available from the Gene Expression Omnibus (GEO) database under accession GSE251668

[0092] (www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE251668). The data was obtained using duet multiomics solution evoC - a new sequencing technology that simultaneously derives all four genetic bases without ambiguity in C or T calls alongside the modification status of Cytosines (5mC and 5hmC, i.e, 6-Letter data) in a single read from a single DNA molecule. The technology consists of pre-sequencing library prep and post-sequencing analysis pipeline, providing single-base resolution of genetics and epigenetics at high accuracy. The technology is described in Fullgrabe et al. 2023 and was implemented as described. We sequenced DNA from a mouse embryonic stem cell-line, ES-E14TG2A at high depth to obtain simultaneous reading of the genome and epigenome (both 5mC and 5hmC). The data was processed using a software developed in house as a python package to efficiently analyse the quant files and extract relevant information regarding the methylation status of each CpG site that we cover.

[0093] Feature engineering. Features were created for each gene analysed based on genomic and methylation information (e.g. mean methylation, number of CpG, length of the region) at different genomic annotations. Gene annotations from Gencode for GRCm38.p6 were used (GENCODE - Mouse Release M25, see www.gencodegenes.org / mouse / release_M25.html). A set of regions of interest were identified for each gene including promoters, TSS (region around the transcription start site), 5' UTR, exons, introns, 3' UTR, and 5kb downstream from the gene. Promoters and the downstream regions were split into sub-regions of equal length to enable proximal / distal regions to have a different effect in the model. For each region, the following features were calculated: (i) mean mC fraction; (ii) mean hmC fraction; (iii) mean modC fraction; (iv) length of the region; (v) CpG count in the region. Additional features were tested but not used in the results shown in the present examples, including density of mC, hmC or modC (e.g. mean mC*count / length; this captures information already captured by the combination of the corresponding mean, CpG count and region length and therefore could be used instead of one or more of these), entropy of mC, hmC or modC (which could be used to measure the degree of variability of the methylation across a region), and the standard deviation of the mC fraction, hmC fraction and / or modC fraction. Entropy can be calculated by calculating a probability density function 'p' for mC, hmC or modC over the region, which is essentially a histogram representing the methylation fraction (mC fraction), hydroxymethylation fraction (hmC fraction) or modified C fration (modC fraction) in each of a plurality of bins over the region. The entropy is then obtained as the sum over k of -p_k*log(p_k), where p_k would represent the probability of bin k (i.e. entropy=fc-pfc(log (pfc)) where pkis the probability for bin k). A complete list of features used to generate the results shown, and their description, is provided in Table 1.

[0094] The mean 5mC fraction is the mean over a defined region of the fraction of reads at each CpG in the region that have a 5mC relative to all reads mapping to the C. The mean 5hmC fraction is the mean over a defined region of the fraction of reads at each CpG in the region that have a 5hmC relative to all reads mapping to the C. The number of CpGs is the number of CpG dinucleotides in the region. The length of the region is the length of the region in number of bases (bp). This is fixed for regions of fixed length (e.g. bins of fixed length in promoter regions, TSS) but is variable for regions like exons, introns and UTRs. For each region, a mean modC fraction feature was also calculated as the mean over a defined region of the fraction of reads at each CpG in the region that have a 5mC or a 5hmC (i.e. conflating 5mC and 5hmC) relative to all reads mapping to the C. This feature is termed mean modC and can be used e.g. in comparative models that do not differentiate between mC and hmC (or as an additional feature in models that do include hmC information). It is not used in the results shown. Table 1. Features used and their description. Individual features of the model are the “features calculated in the region”. Regions are defined individually for each gene based on annotations from Gencode as explained above. The total model comprises 4 features calculated in each of 23 regions, i.e. 92 features per gene.

[0095] Annotations and missing features. The data was analysed for missing features to ensure that the chosen features could be calculated for sufficient numbers of genes. 85% of genes contained a 5’UTR, 96% contained a 3’UTR, and 99% contained an exon. 7.5% of genes had missing CpGs in the 5’UTR, 8.5% had missing CpGs in 3’ UTRs, and 0.2% have missing CpGs in exons (where “missing” CpGs are CpGs with not reads mapping to them). 3% of genes had a mean coverage lower than 5 and only those were excluded. The GENCODE annotation is made by merging the manual annotation produced by the Ensembl-Havana team (levels 1 and 2 annotations, level 1 being validated annotations and level 2 being manual annotations) and the Ensembl-genebuild automated gene annotation (level 3). Regions that were annotated as a gene by manual annotation but where the corresponding annotation was from automated annotation were discarded, resulting in the loss of a small number of genes. Soft-masked genes were not included (i.e. genes in soft-masked regions of the genome which are regions that are difficult to sequence accurately due to repetitive sequences). These represent very small numbers of genes. Only manual annotations from Gencode (Havana) were kept (level 1 and level 2 annotations in Gencode), with priority for level 1 annotations when available. If a region has 0 CpGs, the mean of the region is set to 0 for all of 5mC, 5hmC and modC. In the rare case where a gene has a missing region (e.g. a 5' UTR), feature values (other than length) were imputed as the mean of the respective feature in all other genes at that region. This strategy has very little effect on the overall R2, with variations of the order of I O3. A total of 11449 genes for which expression data was available remained after filtering. This number represents the overlap between the genes that we selected above from our annotation database, and the genes for which we simultaneously had gene expression.

[0096] Multiple transcripts. The GENCODE annotation datasets often gives multiple transcripts per genes, and each transcript will have its own set of exons, introns, 5’ and 3' UTRs. The bulk RNA seq dataset that was used for training does not give any information about the transcript (e.g., a transcript ID). Instead it gives the start and end position of the gene for which expression is measured, as well as gene name. In order to map our annotations to a gene, we therefore decided to compute features (mean methylation fraction, etc) across the longest transcript of a gene, when multiple transcripts are present. The TT-seq data did provide transcript specific information but for consistency (and to avoid having to recalculate all of the features), the data was summarised across transcripts for the same gene. It is however possible to train a model using transcript specific data, resulting in transcript specific predictions. Note that the model could have equally been trained at a transcript level had the information been available in the specific data used.

[0097] Units. The bulk RNA-seq dataset used to train the model (GSE135509) gives the expression in log2(RPM) (reads per million). The TT-seq data used was also provided in RPM. The RPM is defined as the number of reads that map to a given gene divided by the number of reads matched to the genome sequence, and then multiplied by 106. Another commonly used unit is the RPKM (Reads per kilo base per million mapped reads). It’s essentially the RPM normalised by the gene length. In practice, the gene length itself is normalised by 1 kbp, so that RPKM = RPM * 1000 / gene length. Other studies use TPM (transcripts per million). TPM can be derived from RPKM as TPM = RPKM * 10*6 / sum (RPKM). Hence, because they only differ by a normalisation factor, TPM and RPKM are linearly correlated, and in principle TPM and RPKM could be used with similar accuracy. This is illustrated on Figure 4, which shows the relation between RPM, RPKM, and TPM in the bulk RNAseq data used (top panels) as well as the distribution of RPM, RPKM, and TPM in this data (bottom panels).

[0098] The data shows an expected relationship between these units. The excess of genes with values of 0 in RPM occurs because RPM simply counts reads that map to a transcript, it penalises small genes, because not a lot of reads will map to them (just by virtue of them being small). Hence a gene could have a low RPM because it’s not expressed, but also because it’s small. The data shown for the TPM and RPKM applies a cutoff at log2(RPM) of about -2.5. When converting to RPKM, in Iog2 space, this substracts Iog2(length) to the RPM value, Because different genes will have different lengths, this transformation smoothes out the data. The present inventors found the use of a unit normalised by length (such as RPKM or TPM) to improve performance of the model, and therefore the models were trained using RPKM values.

[0099] Dynamic range. Genes show a high level of variability in dynamic range. The model demonstrated in the results below was trained to predict an absolute RPKM value for each gene. The model can be improved by correcting for gene specific dynamical range. In particular, the model can be configured to predict an output between predetermined bounds (e.g between 0 and 1 would be a natural choice), and this can be mapped to a gene-specific dynamic range using a predetermined function (e.g. linearly). The gene-specific dynamic range can be obtained by determining the minimum and maximum observed expression values across the training data, or across an expanded dataset including the training data and one or more of the data in GEO: GSE90277 (www.ncbi.nlm. nih.gov / geo / query / acc. cgi?acc=GSE90277, ENCODE Project Consortium 2012), GEO: GSE36025 (www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE36025, Lin et al. 2014) and GEO: GSE136199 (www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE136199, Kraushar et al. 2021) (after transforming the data to be in the same units as explained above).

[0100] Model training. Once a list of features is obtained, the data comprising values for each of the epigenetic features for each of the genes observed in the RNAseq data and corresponding ground truth RNA expression values was passed to a machine learning algorithm to predict gene expression. Two types of architectures were used: a random forest (RF) regressor from scikit learn (scikit- learn.org / stable / modules / generated / sklearn. ensemble. RandomForestRegressor.html) and a gradient boosted tree model (XGBoost) model from PyPI (xgboost 2.0.3, pypi.org / project / xgboost / ). All results shown here are using the gradient boosted regressor from XGBoost. Best hyperparameters were identified through a grid search using the following parameters: (i) max_depth (max depth of tree): searched between 3 and 7, selected value of 5; (ii) n_estimators (number of trees): searched between 100 and 600, selected 500; (iii) Eta (learning rate, the step size at each iteration while moving toward a minimum of the loss function): searched between 0.01 and 0.05, selected 0.03; (iv) Subsample (fraction of the training data to be randomly sampled for building each tree) searched between 0.2 and 0.6, selected 0.5; (v) colsample_bytree (fraction of features to be randomly sampled for building each tree), searched between 0.8 and 1 , selected 0.9. Results are shown for models with the “selected” values of these parameters. All genes on chromosome 8 were used as a test set (468 genes), and the rest were used as a training set (10981 genes). Note that classifiers could have been used instead of regressors. Both of the above types of models have classifier equivalents (see e.g. scikit- learn.org / stable / modules / generated / sklearn. ensemble. RandomForestClassifier.html). Further, any machine learning model that is suitable for regression (linear or non-linear, although models that are able to capture non-linear relationships, such as RF models, XGboost and neural networks, are preferred) or classification could have been used. When classification is used, genes can be classified between a plurality of categories (e.g. 4), corresponding to non-overlapping ranges of the RNA expression values observed in the training data (or an expanded version thereof as discussed above in relation to dynamic range). The following categories are suitable choices: not expressed, poorly expressed, moderately expressed and highly expressed, where the range for each category can be chosen to match the range of the first, second, third and fourth quartiles of the distribution of expression values in the training data set or an expanded version thereof as discussed above. Alternative approaches forselecting ranges can be used, such as e.g. by clustering expression values to identify groups of genes and boundaries between them. Performance of the models was assessed in terms of the RA2 (variance explained) and Spearman R between predicted and observed expression in the test set. Further, in order to assess the impact of including the 5hmC information, a comparative model was trained using exactly the same procedure but using only mC (i.e. using a mean mC fraction feature alone instead of the mean mC fraction and mean hmC fraction features).

[0101] Feature contribution. The contribution of each region to the performance of the model was evaluated by rerunning the model one region at a time (e.g. exons), and recording the R2value of that simulation (in practice, we run the model 10 times per region, and took an average R2over those 10 instances). We did that for mC only (in which case the features per region are mean mC, CpG count, and region length), and for the usual mC+hmC model (in which case the features per region are mean mC, mean hmC, CpG count, and region length). This allows us to pull out the contribution of hmC to each region by seeing how many RA2 points we gained by going from an mC-only to an mC+hmC model.

[0102] Results

[0103] Figure 5 shows the results obtained by training the XGBoost regression model for bulk RNAseq expression prediction from the complete set of features in Table 1 . The model obtained a R2=0.75 and Spearman rank correlation^.86. A comparative model with only the mC information (no 5hmC) obtained a R2=0.73 and Spearman rank correlation^.85. This shows that the addition of 5hmC information enabled to explain an additional 3% of the unexplained variability between predicted and observed expression (an improvement of about 18% in explaining the 17% unexplained variability). This is a significant gain particularly when considering tat not all of the 17% unexplained variability can be explained by a predictive model as inherent biological variability imposes a theoretical maximum variability that can be explained by a predictive model (likely in the region of 90 to 95%). Therefore, the true gain in explainable variability is likely to be even higher than 17%.

[0104] Figure 6 shows the results of a feature importance analysis for the model used to generate the results in Figure 5. The total length of the bar is the contribution of the features of the region as provided in Table 1 , and the lighter part of the bar is the contribution of the features of the region where mean mC fraction is used instead of mean mC fraction and mean hmC fraction. This analysis shows that regions upstream and downstream of the gene become less and less important as their distance to the gene increases. In particular, the data shows that promoter regions between 0 and 1400bp from the TSS were more informative than further upstream regions, and a model that excludes the regions further upstream of these would likely have a similar accuracy. Further, the main contribution to the RA2 appears to come from introns. Finally, in the 3' UTRs and part of the downstream regions (e.g. 1-2 kb downstream of the 3’UTR and 4- 5kb downstream of the 3’UTR), the contribution from hmC dominates that of mC. Even though these regions are overall of small influence on the RA2, collectively all regions showing significant 5hmC influence result in the above prediction performance improvement.

[0105] The same approach as illustrated above for bulk RNA-seq data was also applied to train a model to predict TT-seq (transient transcriptome sequencing) expression data. TT-seq (Schwalb et al. 2016) is a protocol that measures transcription for a limited period of time by labelling newly synthesized RNA through metabolic incorporation of 4-thiouridine (4sU) in live cells. Analysis of this RNA can provide information on the velocity of transcription, whereas bulk RNA-seq relates to a steady-state between transcription and decay. The inventors postulated that methylation is more likely to give an indication of the rate of transcription rather than the steady state (since the decay rate is independent of the epigenome whereas the rate of transcription is not). Therefore, the inventors set out to test whether methylation can be used to predict TT-seq expression better than it predicts bulk RNA-seq expression.

[0106] The results of this are shown on Figure 7, which shows the results obtained by training the XGBoost regression model for TT-seq expression prediction from the complete set of features in Table 1 . The model obtained a R2=0.85 and Spearman rank correlation^.91 . A comparative model with only the mC information (no 5hmC) obtained a R2=0.83 and Spearman rank correlation^.90. To the best of the inventors knowledge, this is the first time that prediction of TT-seq expression from methylation data has been attempted. The results show an extremely high prediction accuracy, achieving levels of explained variability likely close to the theoretical maximum, with mC alone, and even more so with mC and hmC.

[0107] Figure 8 shows the results of a feature importance analysis for the model used to generate the results in Figure 7. Similar to the bulk RNA-seq analysis, the inventors evaluated the importance of each region by assessing their independent contribution to the RA2. Interestingly, hmC makes a significant contribution in many regions. However, the main gain in RA2 is driven by the contribution of introns, and is largely independent of hmC levels. TT-seq data contains many reads in genomic regions that give rise to shortlived, non-coding RNAs, such as introns (e.g. this paper), which could explain the spike in the introns contribution, although it remains to be understood why hmC is so prevalent in other regions. This could indicate a prominent role of 5hmC in the onset of gene expression.

[0108] EXAMPLE 2 - Prediction of RNA transcription usinq a sequence model

[0109] In the examples above, the inventors used feature based classification and regression algorithms to predict orthogonal types of omics data related to gene expression, using epigenetic features derived from 6-letters sequencing data (illustrated on Figure 9, right hand side). However, the data used (6-letters sequencing data) inherently contains more information than just epigenetic data, since it contains genetic and epigenetic information on a single read. Genetic information (i.e. the DNA sequence in traditional 4 bases code - A, C, T, G) has been shown to be predictive of gene expression in Avsec et al. 2021. Building on the findings in Example 1 above, the present inventors set out to investigate whether 6-letters sequencing data could be used to train improved models for predicting gene expression.

[0110] The designed an approach illustrated on Figure 9, left hand side, where instead of using the 6 letters sequencing data from the indicated regions to calculate features, 6 letters sequence data in the entire region around a TSS (e.g. from 2kb upstream of the TSS to 5kb downstream of the 3’UTR, or in 2000bp windows centred at TSSs (Transcription Start Sites)) is used to obtain an input sequence provided as an input to a sequence based deep learning model trained to predict a gene expression metric from bulk RNA seq or TT-seq. Note that other / longer regions of the genome could be used instead, but the above regions were found to be informative in Example 1 , and information content was shown to decrease with larger distances from the genes such that there is likely to be diminishing returns in including longer sequences as inputs (considering the increased model complexity required when increasing the number of input tokens).

[0111] Model architectures. A transformer based deep learning model architecture can be used due to its computational efficiency in handling large input strings. Specifically, an architecture similar to that used in Avsec et al. 2021 (see Avsec et al. 2021 , Methods) can be built, including a one or more convolutional layers and one or more transformer layers, as well as a regression head that is organism, cell type or tissue specific depending on the training data available. The convolution and transformer layers would together produce embeddings which would be used by the regression head to predict the expression metrics (either individually or using a single head with multiple outputs to predict a plurality of expression metrics, such as e.g. bulk RNA expression, and / or TT-seq expression). Alternatively, a model using convolutional layers only can be used. For example, CNN models can be built which consist of three primary components: (i) an initial 1 -layer convolutional block, (ii) a 9-layer residual convolutional stack, and (iii) a final convolution and cropping output layer. Each convolutional layer can consist of a 1 D convolution followed by a ReLU activation function. The initial ‘stem’ block is used to upscale the initial channel dimensions (5 initial channels for 4-base encoding, 8 initial channels for 6-base encoding) to 64 channels (kernel_size=17, stride=1 , padding =“same”). Following this initial upscaling, the data flows through a 9-layer 1 D convolutional stack; all layers within this stack have 64 channels (kernel_size=17, stride=1 , padding - ‘same”). This stack has residual connections throughout and also uses exponential dilation (exponent of 2), such that the span of the dilation increases from 0 to 512 throughout the stack. The final layer first uses a 1 D convolution to downscale the channels from 64 to 1 (kernel_size=3, stride=1), followed by a central cropping layer to restrict the output dimensions prior to flattening for a final output target tensor with dimensions (1 , 2000) (when using 2kbp focal regions, see below). The model may not use any pooling operations up to this point, in order to maintain the single-base-resolution of the models Then, a classification or regression layer (e.g. a fully connected layer with one output node for regression, or one or more output nodes for classification depending on the number of classes) can be used to produce a gene expression output.

[0112] Alternative architectures were also designed by the inventors, including a transformer and ConvFormer architecture. The Transformer models can include three primary components: (i) a layer normalization block, (ii) a configurable-depth transformer decoder stack, and (iii) a final prediction layer with GELU activation. Each transformer decoder layer can consist of multi-head self-attention followed by a feedforward network with ReLU activation. The initial layer normalization block processes the input tensor of shape (batch_size, sequencejength, num_channels) to standardize the features. Following this normalization, learnable positional encodings are added to the input to provide position information to the transformer (shape: 1 , sequencejength, num_channels). The data then flows through the transformer decoder stack; all layers within this stack maintain the input dimension (sequencejength * num_channels) and use a causal attention mask to ensure each position only attends to previous positions. The number of attention heads is dynamically determined based on the number of channels (this can default to matching the number of attention heads to the number of channels). The feedforward networks within each transformer layer can have 32 hidden units, use ReLU activation, and apply dropout (probability=0.5) . The final prediction layer first applies mean pooling across the tracks dimension, followed by a linear transformation to predict the target at each position, producing an output tensor with dimensions matching the input sequence length. A GELU activation can be applied as the final non-linearity. This can be used as input to a final classification or regression layer as explained above. Alternatively, mean pooling can be applied at the final prediction layer across the positions to obtain a single prediction.

[0113] ConvFormer models can simply take a combination of convolutional and transformer blocks as described above.

[0114] Generating feature and target data. The present example focusses on the task of predicting functional genomic readouts in windows constructed around at TSSs (Transcription Start Sites), referred to below as “focal regions”. The inventors generated focal regions using GENCODE gene annotations for the M25 release (www.gencodegenes.org / mouse / release_M25). Specifically, they downloaded the “basic annotation” and parsed strands separately to identify the transcription start sites as the starts of protein coding genes, while matching the strand orientation. After identifying the start sites they then got the region coordinates for regions around the TSS (e.g. 2000bp regions centred around the start sites (+ / - 1000bp) or regions comprising 2kb upstream of the TSS to 5kb downstream of the 3’UTR, as in Example 1 ). In total, they generated focal region coordinates for 21 ,674 TSSs. The inventors next extracted the GRCm38 reference genome sequence at each of these focal genomic regions and then performed a one-hot encoding of these DNA sequences to generate inputs compatible with our model architecture (A=[1 ,0,0,0]; C=[0,1 ,0,0]; G=[0,0,1 ,0]; T=[0, 0,0,1]; N=[0, 0,0,0]). Alongside the genomic one-hot encoding, the inventors added some additional information. As a standard data augmentation approach the inventors appended the reverse complement for each one-hot sequence to the input in order to capture the alternative strand state (such that we represent sense and anti-sense sequences) and added an additional strand element to the one-hot arrays, to indicate from which strand transcription initiates. This fully captures information that would be represented in a genomic-only model (4-base encoding). This is similar to the encoding illustrated on Fig. 10A, except with an additional dimension to encode the strand, i.e. 5 bits per base, 2000 bases per input sequence on each strand (assuming a focal region of 2kbp), i.e. 4000 bases, leading to an encoded dimension of (5, 4000) for each focal genomic region. At the simplest, the model can take as input a modified one-hot encoded DNA sequence in the form of A=[1 , 0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C= [0,0,0 ,f1 ,f2,f3], and N=[0, 0,0, 0,0,0], where f1 is the fraction of unmodified C reads at the base, f2 is the fraction of methylated C at the base, and f3 is the fraction of hydroxymethylated C at the base. Alternatively, a one-hot encoded DNA sequence in the form of A=[1 ,0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C=[0, 0,0,1 ,0,0], mC=[0,0,0,0,1 ,0], hmC=[0, 0,0, 0,0,1], and N=[0, 0,0, 0,0,0] could be used with the base with the highest fraction of c, mC and hmC called at the location. The former encoding is expected to be more informative as it preserves the information about the fractions of mC and hmC.

[0115] 6-base data for the prediction of gene expression can be obtained using the using duet multiomics solution evoC as described above. The inventors obtained one such dataset to investigate encoding of 6-base data. The inventors generated four technical replicate samples from the mouse embryonic stem cell-line (ESEI 4TG2A), extracted DNA (80ng input), prepared libraries using a duet multiomics solution evoC kit, and sequenced with an S4 2x 150bp kit on an Illumina Novaseq 6000. They performed initial post-sequencing analysis (read resolution, trimming, alignment to GRCm38, quantification and QC) using the duet multiomics analysis pipeline to generate methylation quantification files for each of the 4 samples. Mean coverage of samples ranged from 35 to 50X. For quantification they pooled these technical replicates to realise a single biological sample with approximately 100x CpG coverage. Pooling was performed only because the samples are technical replicates of the same cell pool, and therefore treating them as independent samples would have amounted to pseudo-replication. Such high levels of coverage are not necessary to train the models. Indeed, the inventors have done some downsampling experiments showing that the performance of the models is robust to lower coverage of methylation data at least as low as 15- 20x coverage.

[0116] For genomic + evoC methylation models (6-base encoding) the inventors designed a scheme that uses the 4-base encoding described above and additionally overlaid: (i) the number of C calls, (ii) the number of 5mC calls and (iii) the number of 5hmC calls, which were quantified in the 6-base data, at all CpG locations. They used strand-specific methylation quantification in order to specifically overlay these counts to each of the sense and anti-sense components of the input sequences. This is similar to the encoding illustrated on Fig. 10B (Methylated C example 1), except with an additional dimension to encode the strand, i.e. 8 bits per base, 2000 bases per input sequence on each strand (when 2kbp focal regions are used), i.e. 4000 bases, leading to an encoded dimension of (8, 4000) for each focal genomic region. Alternative encodings are possible, such as e.g. using fractional counts instead of counts, as illustrated on Fig. 10B, Methylated C example 2. Further, alternative encodings may use statistics from model based estimates of fractions of modified and unmodified bases, to account for non-biological noise owing to sampling variation and methylation miscalls that may influence counts / fractional counts data. For example, a model based on a Dirichlet process can be used, which can capture uncertainty in multinomial count data. For example, one possible encoding would be to use the 95% confidence interval for the fraction of each possible state (C, 5mC, 5hmC), or the mean and variance of fractional estimates, as illustrated on Fig. 10C. The mean (also known as expectation) and variance of the Dirichlet can be calculated exactly using statistical formulae, while confidence intervals can be be estimated empirically using a sampling approach. For a Dirichlet distribution with X=(Xi,..., Xk)-Dir(a), where Xi,..., Xk are K fractions that sum to 1 and a =(ai,..., ak) are concentration parameters for each of the k categories, where ai>0, the mean for each category can be calculated as and the variance can be calculated as where a0= at. This is similar to the encoding illustrated on Fig. 10C (Methylated C example 1), except with an additional dimension to encode the strand, i.e. 11 bits per base (4 for A, C, G, T, and 2 for each of C, mC, hC fractions, and 1 for the strand), 2000 bases per input sequence on each strand (when 2kbp focal regions are used), i.e. 4000 bases, leading to an encoded dimension of (11 , 4000) for each focal genomic region. In methylated C example 3 on Fig. 10C, the order of the encoded information after the 4 bits encoding the genetic information (N, A, C, G, T) is set out as: [C fraction left boundary of Cl, C fraction right boundary of Cl, mC fraction left boundary of Cl, mC fraction right boundary of Cl, hmC fraction left boundary of Cl, hmC fraction right boundary of Cl], where “Cl” stands or confidence interval. In methylated C example 4 on Fig. 10C, the order of the encoded information after the 4 bits encoding the genetic information (N, A, C, G, T) is set out as: [C fraction mean, C fraction sd, mC fraction mean, mC fraction sd, hmC fraction mean, hmC fraction sd], where “sd” stands for standard deviation. However, the exact order of each of these bits of information is not essential and any order can be used provided that the order is consistently used. The same applies to methylated C examples 1 and 2 on Fig. 10B. It is also irrelevant whether the epigenetic information is provided in the first few bits or the last few bits of the encoding, as long as the location of each bit of information is consistent. Similarly, the specific choice of order for the one-hot encoding of N, A, C, G, T is also arbitrary, and any other order could be used (e.g. where A is [0,1 ,0,0] and C is [1 , 0, 0, 0]). The orders presented are intuitive and therefore advantageously easier to interpret and to encode with minimum risk of confusion.

[0117] The input length can be reduced to the minimum input size that enables to include the regions in Example 1 (i.e. gene +2kb upstream and 5kb downstream), using padding characters for genes of size such that the gene+upstream and downstream sequences is smaller than the predetermined input size (e.g. left or right padding). Alternatively, the input length can be set to the minimum input size that enables to include the regions in Example 1 (i.e. gene +2kb upstream and 5kb downstream) and sequence including the gene and upstream and downstream sequence to fill the input length can be used. Human transcript lengths vary between a few dozen bases to 50k bp, with an average around 200bp. Thus, such an input length would still likely be significantly smaller than the 200kb used in Avsec et al. 2021 , in which case the number of convolution layers can be reduced, for example including one, two or three convolution layers instead of 7). Alternatively, a 200kb length of sequence can be used (preserving the input dimensions in Avsec et al. 2021) that includes the 7kb of sequence defined in Example 1 and additional sequence upstream and downstream. The first alternative above (gene + predetermined upstream and downstream sequence length, with padding characters) may be advantageous in order to maintain a gene centric approach, as that is what the ground truth training data (expression data from bulk RNA-seq or TT-seq) contains - whereas 200kb may well contain multiple genes.

[0118] The regression head (note that a classification head could be used instead, as explained in Example 1) will be substantially simpler than that of the Enformer as it would only output expression metrics for a particular gene (i.e. a few single scalar values rather than multiple genomic tracks at 200bp resolution). For example, a single pooling pooling layer, a ReLU or GeLU layer and a softmax layer may be included.

[0119] The model will be trained using 6 letters sequence data and bulk RNA sequence data (or TT-seq data) for a plurality of samples from the same species, ideally for a plurality of cell lines and / or tissues. Note that Enformer did not use a gene centric approach and did not predict gene expression of single genes. Instead, Enformer takes as input bins of 200kb of genomic sequence and produces as output transcription factor chromatin immunoprecipitation and sequencing (TF ChlP-seq), DNAse-seq, ATAC-seq and CAGE genomic tracks for the same region. By contrast, the proposed model will use gene centric stretches of 6- letters sequence data as input and produce as output single scalar expression (or accessibility) values for the gene. Thus, the proposed model is fundamentally different in structure and training data, even though it may use a similar architecture.

[0120] References

[0121] A number of publications are cited above in order to more fully describe and disclose the invention and the state of the art to which the invention pertains. Full citations for these references are provided below. The entirety of each of these references is incorporated herein.

[0122] The ENCODE Project Consortium., Moore, J.E., Purcaro, M.J. et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature 583, 699-710 (2020).

[0123] Alda-Catalinas C, Bredikhin D, Hernando-Herraez I, Santos F et al. A Single-Cell Transcriptomics CRISPR-Activation Screen Identifies Epigenetic Regulators of the Zygotic Genome Activation Program. Cell Syst 2020 Jul 22;11 (1 ) :25-41 ,e9. PMID: 32634384

[0124] ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012 Sep 6;489(7414):57-74. doi: 10.1038 / nature11247. PMID: 22955616; PMCID: PMC3439153.

[0125] Lin S, Lin Y, Nery JR, Urich MA et al. Comparison of the transcriptional landscapes between human and mouse tissues. Proc Natl Acad Sci U S A 2014 Dec 2;1 11 (48):17224-9. PMID: 25413365

[0126] Kraushar ML, Krupp F, Harnett D, Turko P et al. Protein Synthesis in the Developing Neocortex at Near- Atomic Resolution Reveals Ebp1 -Mediated Neuronal Proteostasis at the 60S Tunnel Exit. Mol Cell 2021 Jan 21 ;81 (2):304-322.e16. PMID: 33357414

[0127] Schwalb B, Michel M, Zacher B, Fruhauf K, Demel C, Tresch A, Gagneur J, Cramer P. TT-seq maps the human transient transcriptome. Science. 2016 Jun 3;352(6290):1225-8. doi: 10.1126 / science.aad9841 . PMID: 27257258.

[0128] Cruz-Molina S, Respuela P, Tebartz C, Kolovos P, Nikolic M, Fueyo R, van Ijcken WFJ, Grosveld F, Frommolt P, Bazzi H, Rada-Iglesias A. PRC2 Facilitates the Regulatory Topology Required for Poised Enhancer Function during Pluripotent Stem Cell Differentiation. Cell Stem Cell. 2017 May 4;20(5):689- 705. e9. doi: 10.1016 / j.stem.2017.02.004. Epub 2017 Mar 9. PMID: 28285903.

[0129] Avsec Z, Agarwal V, Visentin D, Ledsam JR, Grabska-Barwinska A, Taylor KR, Assael Y, Jumper J, Kohli P, Kelley DR. Effective gene expression prediction from sequence by integrating long-range interactions. Nat Methods. 2021 Oct;18(10):1196-1203. doi: 10.1038 / s41592-021 -01252-x. Epub 2021 Oct 4. PMID: 34608324; PMCID: PMC8490152. Lou, S., Lee, HM., Qin, H. et al. Whole-genome bisulfite sequencing of multiple individuals reveals complementary roles of promoter and gene body methylation in transcriptional regulation. Genome Biol 15, 408 (2014). doi.org / 10.1186 / s13059-014-0408-0

[0130] Pongor LS, et al. Integrative epigenomic analyses of small cell lung cancer cells demonstrates the clinical translational relevance of gene body methylation. iScience. 2022 Oct 12;25(11 ): 105338. doi: 10.1016 / j.isci.2022.105338. PMID: 36325065; PMCID: PMC9619308.

[0131] Villicana and Bell. Genetic impacts on DNA methylation: research findings and future perspectives. Genome Biology (2021) 22:127

[0132] Frommer, M. et al. A genomic sequencing protocol that yields a positive display of 5-methylcytosine residues in individual DNA strands. PNAS 89, 1827-1831 (1992).

[0133] Vaisvila, R. et al. Enzymatic methyl sequencing detects DNA methylation at single-base resolution from picograms of DNA. Genome Res. 31 : 1280-1289 (2021).

[0134] WO 2022 / 023753 A1 - COMPOSITIONS AND METHODS FOR NUCLEIC ACID ANALYSIS, assignee: CAMBRIDGE EPIGENETIX LIMITED.

[0135] Simpson JT, Workman RE, Zuzarte PC, David M, Dursi LJ, Timp W. Detecting DNA cytosine methylation using nanopore sequencing. Nat Methods. 2017 Apr;14(4):407-410. doi: 10.1038 / nmeth.4184. Epub 2017 Feb 20. PMID: 28218898.

[0136] Pacific Biosciences 2022, “Measuring DNA methylation with 5-base HFi sequencing”, available at www.pacb.com / wp-content / uploads / application-brief-measuring-dna-methylation-with-5-base-hifi- sequencing.pdf

[0137] Fullgrabe, J., Gosal, W.S., Creed, P. et al. Simultaneous sequencing of genetic and epigenetic bases in DNA. Nat Biotechnol 41 , 1457-1464 (2023). doi.org / 10.1038 / s41587-022-01652-0

[0138] Takahashi H, Kato S, Murata M, Carninci P. CAGE (cap analysis of gene expression): a protocol for the detection of promoter and transcriptional networks. Methods Mol Biol. 2012;786:181-200. doi:

[0139] 10.1007 / 978-1 -61779-292-2_11. PMID: 21938627; PMCID: PMC4094367.

[0140] Vaswani et al. 2017. Attention Is All You Need. arXiv:1706.03762

[0141] Yu & Koltun 2016. Multi-scale context aggregation by dilated convolutions. arXiv:1511 .07122v3

[0142] For standard molecular biology techniques, see Sambrook, J., Russel, D.W. Molecular Cloning, A Laboratory Manual. 3 ed. 2001 , Cold Spring Harbor, New York: Cold Spring Harbor Laboratory Press

Claims

Claims:

1. A computer-implemented method of predicting the gene expression of one or more genes in a sample, the method comprising: receiving (110) sequence data associated with the one or more genes obtained from the sample, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine; and predicting (130), for each gene of the one or more genes and using the sequence data associated with the gene, the value of a gene expression metric associated with the gene, wherein the predicting is performed using a machine learning model trained to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a gene expression metric for the gene.

2. The method of any preceding claim, wherein the machine learning model has been trained using a training dataset comprising expression data obtained by RNA sequencing, and the predicted gene expression metric is a gene expression metric obtainable by RNA sequencing.

3. The method of claim 2, wherein the RNA sequencing is an RNA sequencing technology indicative of steady state transcript levels, optionally bulk RNA sequencing.

4. The method of claim 2, wherein the RNA sequencing is an RNA sequencing technology indicative of newly synthesised RNA over a predetermined period of time, optionally TT-seq.

5. The method of any preceding claim, wherein the gene expression metric is a metric derived from read counts, optionally RPM (reads per million), RPKM (reads per kilo base per million mapped reads), or TPM (transcripts per million), or a metric derived therefrom by log transformation, optionally wherein the gene expression metric is a log(RPKM) or log(TPM).

6. The method of any preceding claim, wherein the machine learning model is a regression model, and the machine learning model is configured to output a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene.

7. The method of any of claims 1 to 6, wherein the machine learning model is a regression model, and the machine learning model is configured to output a value between predetermined bounds, optionally between 0 and 1 , and predicting the value of a gene expression metric associated with the gene further comprises obtaining a gene expression metric indicative of the absolute level of one or more transcripts associated with the gene using a predetermined function derived from the observed dynamic range of the one or more transcripts in a plurality of samples.

8. The method of any of claims 1 to 5, wherein the machine learning model is a classification model, and predicting the value of a gene expression metric associated with the gene comprises classifying the gene, using the machine learning model, between a plurality of classes associated with respective ranges of values of the gene expression metric, optionally wherein the plurality ofclasses and associated ranges correspond to respective quartiles of a distribution of values of the gene expression metric obtained from a plurality of samples.

9. The method of any preceding claim, wherein the machine learning model takes as input a plurality of features derived from sequence data including at least one feature indicative of the presence of hydroxymethylated cytosine, and wherein the method further comprises determining (120), using the sequence data, values for each of said plurality of features.

10. The method of claim 9, wherein the plurality of features include, for each of a plurality of genomic regions associated with the genes, one or more of: a feature indicative of the presence of hydroxy methylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the presence of cytosines that are either methylated or hydroxymethylated in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region.

11. The method of claim 10, wherein the plurality of features include: a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, a feature indicative of the length of the region, and a feature indicative of the number of CpGs in the region.

12. The method of claim 10 or claim 11 , wherein the plurality of features include the same features for all regions or a different set of features depending on the region, provided that at least one of the regions includes a feature indicative of the presence of hydroxymethylated cytosines in the region and all regions include at least one of a feature indicative of the presence of hydroxymethylated cytosines in the region, a feature indicative of the presence of methylated cytosines in the region, and a feature indicative of the presence of methylated or hydroxymethylated cytosines in the region.

13. The method of any of claims 10 to 12, wherein the epigenetic data comprises sequence reads, and wherein: the feature indicative of the presence of hydroxymethylated cytosines in the region is a summarised value, over one or more CpGs in the region, of the fraction of reads indicative of the presence of a hydroxymethylated cytosine at the CpG.

14. The method of any of claims 10 to 13, wherein the feature indicative of the presence of methylated cytosines in the region is a summarised value, over one or more CpGs in the region, of the fraction of reads indicative of the presence of a methylated cytosine at the CpG; optionally wherein a summarised value is a mean value.

15. The method of any of claims 10 to 14, wherein the plurality of regions include one or more or all of: one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene, a region corresponding to a promoter like sequence, a region of a predetermined length centred on the transcription start site of the gene, a region corresponding to the 5’UTR of the gene, a region corresponding to the first exon of the gene, a region corresponding to the first intron of the gene, a region combining regions corresponding to all exons of the gene, a region combining regions corresponding to all introns of the gene, a regioncorresponding to the 3’UTR of the gene, and one or more regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene.

16. The method of claim 15, wherein the one or more regions of predetermined lengths at a predetermined distance upstream of the transcription start site of the gene comprise a plurality of non-overlapping regions of equal length between 10Obp and 500bp, such as about 200bp, together covering a region up to 2kb, 1.8kb, 1.6kb, 1 .5kb, 1.4kb, 1.2kb or 1 kb from the transcription start site of the gene.

17. The method of claim 15 or claim 16, wherein the one or regions of predetermined lengths at a predetermined distance downstream of the 3’ end of 3’UTR of the gene comprise a plurality of nonoverlapping regions of equal length between 500bp and 1500bp, such as about 1 kb, together covering a region up to 5kb, 6kb, 7kb or 8kb from the 3’ end of 3’UTR of the gene.

18. The method of any preceding claim, wherein the machine learning model is a random forest model, a gradient boosted tree model, or a neural network model.

19. The method according to any of claims 1 to 8, wherein the machine learning model is a deep learning model trained to take as input sequence data associated with a gene and produce as output a gene expression metric forthe gene, wherein the sequence data further comprises genetic data associated with the gene.

20. The method according to claim 19, wherein the deep learning model is a transformer based model, a model comprising one or more convolution layers, a model comprising one or more transformer layers, or a model comprising one or more convolution layers and one or more transformer layers, and / or wherein the deep learning model comprises a classification or regression layer.

21. The method of claim 19 or claim 20, wherein the deep learning model takes as input a sequence derived from 6-letter sequence data associated with a gene and the deep learning model has been trained using training data comprising sequences derived from 6-letter sequence data associated with a plurality of genes and expression data associated with said plurality of genes, optionally wherein the input sequence is a sequence derived from the 6-letter sequence data encoded using for each base of the sequence the following encoding scheme: A=[1 ,0,0, 0,0,0], G=[0,1 ,0,0, 0,0], T=[0, 0,1 ,0,0,0], C=[0,0,0,f1 ,f2,f3], and N=[0, 0,0, 0,0,0], where f1 is the fraction of unmodified C reads in the 6-letter sequence data at the base, f2 is the fraction of methylated C reads in the 6- letter sequence data at the base, and f3 is the fraction of hydroxymethylated C in the 6-letter sequence data at the base, and / or wherein the deep learning model takes as input a sequence comprising data for a predetermined number of bases corresponding to the sequence of the gene, a sequence of a first predetermined length upstream of the gene and a sequence of a second predetermined length downstream of the gene.

22. The method of any preceding claim, wherein the machine learning model has been trained using a training dataset comprising expression data comprising for each of the plurality of genes, a single value of a gene expression metric, the single value obtained by summarising one or more valuescorresponding to respective transcripts associated with the gene and / or respective samples, and wherein the summarised value is a mean value across a plurality of transcripts and / or samples.

23. A computer-implemented method of training a machine learning model for predicting the gene expression of one or more genes in a sample, the method comprising: receiving (210) a training data set comprising (i) sequence data associated with a plurality of genes, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine, and (ii) expression data comprising for each of the plurality of genes, one or more measured values of a gene expression metric; and training (230) a machine learning model, using said training data, to take as input sequence data for a gene or one or more features derived therefrom including at least one feature indicative of the presence of hydroxymethylated cytosine, and produce as output a predicted gene expression metric for the gene.

24. The method of any preceding claim, wherein the sequence data has been obtained using a sequencing technology from which epigenetic and genetic bases can be called on the same read, and / or wherein the sequence data comprises or consists of reads in 6-letters code.

25. A method of performing a genetic test for a subject, the method comprising: obtaining (1100) sequence data associated with a sample previously obtained from a subject, the sequence data comprising epigenetic data indicative of the presence of one or more epigenetic bases at one or more positions in the genome, the one or more genetic bases including methylated cytosine and hydroxymethylated cytosine; and performing (1200) the method of any preceding aspect to obtain a predicted gene expression metric for one or more genes, optionally wherein the sequence data further comprises genetic data indicative of the presence of one or more genetic bases at one or more positions in the genome, and / or wherein the method further comprises analysing (1300) the predicted gene expression metrics, and / or genetic data by obtaining the value of one or more metrics characterising the sample, optionally wherein the one or more metrics are selected from a diagnostic or prognostic metric.

27. A system comprising at least one processor and a non-transitory computer readable medium comprising instructions that, when executed by at least one processor, cause the at least one processor to perform the method of any of claims 1 to 25.

28. A non-transitory computer readable medium comprising instructions that, when executed by at least one processor, cause the at least one processor to perform the method of any of claims 1 to 25.

Images