Methylation quantitative trait loci identification

Abstract

Computer-implemented methods of identifying a methylation quantitative trait locus are provided The method comprise receiving sequence data comprising sequence reads for a plurality of samples; and identifying, using said sequence data, a statistically significant relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus, thereby identifying a methylation quantitative trait locus, wherein the relationship is allele-specific and the presence of said statistically significant allele-specific relationship indicates that the fraction of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy. Related methods, systems and products are described.

Description

[0001] Methylation Quantitative Trait Loci Identification

[0002] Field of the Disclosure

[0003] The present invention relates to methods of identifying relationships between epigenetic modifications and genetic variants and particularly, although not exclusively, to methods of identifying methylation quantitative trait loci (mQTLs), methods of determining the expected level of base modification at an epigenetic locus, and methods of providing 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 studies epigenetic mark. DNA methylation is typically measured using whole genome bisulfite sequencing (WGBS, see Frommeret al. 1992), enzymatic methyl-sequencing (EM-seq, see Vaisvila et al. 2021) or DNA methylation microarrays. Both 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. (2022) 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”).

[0006] The establishment of DNA methylation marks is influenced by environmental exposures but also by genetic variation. For example, germline genetic alterations have been shown to cause changes in DNA methylation that ultimate impact the risks of developing diseases like cancer (see e.g. Ahmed et al. 2021). One way to study this influence is through mQTL (methylation quantitative trait loci, also referred to as “meQTL”). An mQTL is a genetic locus at which genetic variation is associated with variation in DNA methylation at a particular CpG site. mQTL identification (also referred to as “calling”) relies on finding statistically significant associations between single nucleotide polymorphisms (SNPs) and the level of methylation at candidate CpG sites. Once identified, variation as such loci can be used as a proxi to methylation, thereby enabling the probing of e.g. disease relevant methylation patterns using genetic assays. However, mQTLs can have proximal (cis) or distal (trans) effect, and the numbers of SNPs and CpG in a genome such as the human genome are extremely high, resulting in the need to analyse large populations of individuals in order to be able to confidently identify mQTLs, and imposing practical limits on the sensitivity of the detection (i.e. limiting detection to strong and / or common associations).

[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 genetic and epigenetic information about the same DNA molecule (i.e. combining genetic and epigenetic information on a single read) could significantly improve on the prior art mQTL calling methods, provided that a new method was developed to make use of this information. They devised a method that makes use of such information to identify mQTLs with significantly increased power, by identifying relationships between genetic and epigenetic loci at the level of the genetic locus allele rather than at the level of the genotype at the locus using phased genetic and epigenetic locus data from a plurality of samples. The inventors demonstrated that the approach enabled the identification of mQTLs with higher power than was previously possible, drastically reducing the sample requirement for mQTL identification as well as the range of strength of mQTLs that can be identified.

[0009] Thus, according to a first aspect, the disclosure provides a computer-implemented method of identifying a methylation quantitative trait locus, the method comprising: receiving sequence data comprising sequence reads and / or information derived therefrom including information about the presence and location of one or more genetic variants and modified bases in the sequence reads, for a plurality of samples; and identifying, using said sequence data, a statistically significant relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus, thereby identifying a methylation quantitative trait locus, wherein the relationship is allele-specific and the presence of said statistically significant allelespecific relationship indicates that the fraction of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy. mQTL relationships are typically genotype specific, rather than allele-specific. The presence or absence of a genetic variant at the genetic locus on the same chromosome copy refers to the presence of a genetic variant in phase. In other words, according to the present method a relationship between base modification at an epigenetic locus and in phase genotype at a genetic locus is determined.

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

[0011] The genetic variant may be a single nucleotide variant or polymorphism. The modified base may be a modified cytosine. The epigenetic locus may be a CpG locus.

[0012] The sequence data may comprise a plurality of sequence reads for each sample and genetic variants identified in said sequence reads. The sequence data may comprise a plurality of sequence reads for each sample and the method may comprise identifying one or more genetic variants present in each of said plurality of samples using said sequence reads. For example, the sequence data may include one or more reads file (e.g. a SAM or BAM file, typically one per sample) and a variant calling file (e.g. a VCF file, which may be a joint VCF file listing all variants identified across the plurality of samples). Identifying genetic variants is also referred to as variant calling. The sequence data may further comprise, for each sample, information identifying a plurality of sets of heterozygous genetic loci that are in phase with each other. The method may comprise phasing the sequence data, optionally by applying a read-based phasing method to the sequence reads for each sample. The step of identifying a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus may comprise, for each sample individually: identifying a plurality of haploblocks, wherein each haploblock corresponds to a particular set of alleles for respective genetic loci that are in phase with each other in the particular sample; and obtaining, for each haploblock, a modified read count and an unmodified read count, wherein the modified read count is a count of the number of sequence reads that align to and are consistent with the haploblock and that show a modified base at an epigenetic locus, and the unmodified read count is the count of the number of sequence reads that align to and are consistent with the haploblock and that show an unmodified base at the epigenetic locus, optionally wherein said counts are obtained individually for each epigenetic locus associated with a haploblock. The step of identifying a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus may further comprise: for each of one or more variable genetic loci identified in the plurality of samples, fitting a relationship between the modified read counts, the unmodified read counts and the allele at the variable genetic locus for any haploblock comprising the variable genetic locus.

[0013] Identifying a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus may comprise fitting a regression model where the response variable is the count of reads with a modified base at the epigenetic locus and the predictive variables include the total number of reads overlapping the epigenetic locus and the allele at the genetic locus on a single chromosome copy, wherein only reads that map to the same chromosome copy are counted. The total number of sequence reads overlapping the epigenetic locus and the allele at the genetic locus on a single chromosome copy may be obtained as the sum of the modified read count and the unmodified read count. The allele at the genetic locus on a single chromosome copy may be the allele present at the genetic locus on a haploblock encompassing the genetic locus, and the reads that map to the same chromosome copy may be reads associated with the haploblock. The regression model may be a negative binomial regression model. Reads that map to the same chromosome copy as the allele may be reads that overlap the epigenetic locus and at least one heterozygous genetic locus comprised in the haploblock. Reads that are associated with a haploblock may be reads that overlap the epigenetic locus and at least one heterozygous genetic locus comprised in the haploblock. The regression model may be fitted using data comprising, for each haploblock encompassing the genetic locus in each sample: counts of reads associated with the haploblock in the sample that show the presence of a modified base at the epigenetic locus, counts of reads associated with the haploblock in the sample that do not show the presence of a modified base at the epigenetic locus, and the allele present at the genetic locus in the haploblock in the sample. Reads associated with a haploblock may be reads that at least partially map to the haploblock. The reads may include a portion at the 5’ end of the read that extends beyond the 5’ end of the haploblock, and / or a portion at the 3’ end that extends beyond the 3’ end of the haploblock. These portions that extend beyond the end of a haploblock typically do not comprise any heterozygous genetic loci. Indeed, if they did then the loci are either in phase with the loci in the haploblock and the haploblock would extent to encompass them. The regression model may be a model of the form modC~totalC+SNPallele+x where mode the count of reads with a modified base at the epigenetic locus associated with a single chromosome copy in each sample, total C is the count of reads at the epigenetic locus associated with the chromosome copy, and SNPalllele is the allele at the genetic locus on the chromosome copy, and x represents one or more optional additional variables. In other words, the model may be modC~totalC+SNPallele or modC~totalC+SNPallele+x where x may be absent, or may be a single variable x or a plurality of variables x1 , x2, etc. The additional variables may include one or more of: an interaction term between the totalC and SNPallele variables, and a condition variable (e.g. phenotype or exposure variable such as presence absence or identity of exposure to one or more perturbations), or an interaction term between the totalC, the SNPallele and / or a phenotype or exposure variable (e.g. an interaction term between the SNPallele and a condition variable). The variable SNPallele may be a variable that takes the value of 0 or 1 depending on the allele present at the genetic locus on the chromosome copy. It is irrelevant whether the reference or variable allele are associated with the values of 0 and 1 , but conventionally the reference allele would be associated with the value of 0. A reference allele is the allele that is present in a reference genome and / or the allele that is most common in a particular population I cohort.

[0014] Fitting the regression model may comprise identifying a parameter associated with each predictive variable and a metric of statistical significance associated with each parameter, and a statistically significant allelespecific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus may be identified when the metric of statistical significance associated with the variable that represents the allele at the genetic locus satisfies one or more predetermined criteria. For example, the metric of statistical significance may be a p-value and a statistically significant relationship may be identified when the p-value is below a predetermined threshold (optionally after multiple testing correction). Any multiple testing correction known in the art may be used, such as e.g. using the Benjamini-Hochberg procedure or Bonferroni method.

[0015] The method may comprise fitting the regression model for each of a plurality of candidate pairs of epigenetic and genetic loci, and identifying one or more mQTLs by selecting those pairs of epigenetic loci and genetic loci for which the regression model indicates a statistically significant allele-specific relationship between the fraction of a modified base at the epigenetic locus and the genotype at the genetic locus. The plurality of candidate pairs of epigenetic and genetic loci may comprise pairs of epigenetic and genetic loci that are associated with a haploblock identified using the sequence data. The plurality of candidate pairs of epigenetic and genetic loci may comprise all pairs of epigenetic and genetic loci that are associated with a haploblock identified using the sequence data and where the plurality of samples include samples that differ in their genotype at the genetic locus. Subsets of these pairs may also be used, which may be selected based on one or more predetermined criteria, such as e.g. distance between genetic and epigenetic loci in a pair, prevalence of a genetic variant at the genetic locus, pre-existing knowledge of relevance of the genetic locus, etc.

[0016] The sequence data may have been obtained using a sequencing technology from which epigenetic and genetic bases can be called on the same read. The sequence data may comprise or consist of reads in 5 or 6-letters code. The samples of the plurality of samples may each be associated with a respective subject. The / each subject may be a diploid organism. Each subject may be a eukaryote organism, such as a vertebrate, mammalian and / or a human subject.

[0017] According to a second aspect, there is provided a method of identifying a methylation quantitative trait locus, the method comprising: obtaining sequence data comprising sequence reads for a plurality of samples by sequencing genetic material in said samples; and analysing said sequence data using the computer-implemented method of any embodiment of the first aspect. Obtaining sequence data may comprise using a sequencing technology that provides an output from which epigenetic and genetic bases can be called on the same read.

[0018] According to a third aspect, there is provided a method of designing a genetic test, the method comprising: identifying one or more methylation quantitative trait loci using the method of any embodiment of the first or second aspects, and including the one or more genetic loci of the identified mQTLs in a panel of genetic loci that are measured in the genetic test. The method may be computer implemented, e.g. when using the method of any embodiment of the first aspect. The panel of genetic loci may comprise the one or more genetic loci of the identified mQTLs and one or more disease-associated genetic loci. The one or more epigenetic loci of the identified mQTLs may be disease-associated epigenetic loci. The disease-associated genetic loci and / or disease-associated epigenetic loci may be loci that are associated with the risk of a subject having a disease or disorder, and / or loci that are predictive of a prognosis of diagnosis in a subject. The method may further comprise defining a polygenic risk score associated with one or more of the genetic loci that are measured in the genetic test. The method may further comprise designing reagents for the targeted measurement of the panel of genetic loci. Also described according to the present aspect is a method of providing a genetic test, the method comprising: designing a genetic test using a method according to the present aspect, and manufacturing a genetic test comprising reagents for the targeted measurement of the panel of genetic loci.

[0019] According to a fourth aspect, there is provided a method of determining the expected level of base modification at an epigenetic locus in a subject (e.g. likely methylated fraction at one or more CpG contexts in the subject), the method comprising: receiving a previously identified statistically significant allele-specific relationship between the fraction of a modified base at the epigenetic locus and the genotype at a genetic locus, wherein the presence of said statistically significant allele-specific relationship indicates that the fraction of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy; receiving genotype data about the subject at the genetic locus; and determining the expected level of base modification at the epigenetic locus on each chromosome copy of the subject using the genotype data and the identified allele-specific relationship. The level of base modification may be a methylated fraction (e.g. fraction of methylated and / or hydroxy methylated cytosine).

[0020] The previously identified statistically significant allele-specific relationship may have been identified using the method of any embodiment of the first aspect. Note that an allele-specific relationship that has been identified using embodiments of the first aspect is inherently different from an mQTL relationship of the prior art because it is at the level of single alleles rather than genotype. In other words, an mQTL relationship of the prior art is a relationship between the global methylated fraction at an epigenetic locus and the global genotype at a genetic locus (specified as e.g. 0 / 0, 0 / 1 , 1 / 1). By contrast, an mQTL relationship according to the present disclosure is a relationship between the methylated fraction at an epigenetic locus on a single chromosome copy and the allele on the same chromosome copy (specified as e.g. 0 or 1). The method may comprise identifying the statistically significant allele-specific relationship using the method of any embodiment of the first aspect. In other words, determining the expected methylated fraction at the epigenetic locus is done in an allele specific manner. Thus, determining the global expected methylated fraction at the epigenetic locus, taking into account all (e.g. both) chromosomal copies), may comprise averaging the expected methylated fraction at the epigenetic locus when the allele in the subject’s genotype differ between chromosome copies (e.g. the subject is heterozygous at the genetic locus). Alternatively, the global expected methylated fraction at the epigenetic locus may be identical to the expected methylated fraction at the epigenetic locus any one chromosome copy when all chromosome copies have the same allele at the genetic locus (e.g. the subject is homozygous at the genetic locus). Thus, the step of determining the expected methylated fraction at the epigenetic locus on each chromosome copy of the subject using the genotype data and the identified allele-specific relationship may comprising determining a single expected methylated fraction at the epigenetic locus (the single expected methylated fraction being identical for all chromosome copies as they all contain the same allele at the genetic locus). A global methylation fraction is a methylation fraction that is determined by combining information across all chromosome copies at the epigenetic locus.

[0021] The method may further comprise receiving an observed methylated fraction at the epigenetic locus, and comparing the observed methylated fraction with the expected methylated fraction. The observed and expected methylated fractions may be global methylated fractions at the epigenetic locus. Therefore, a global expected methylated fraction may be obtained from allele-specific expected methylated fractions determined as described herein, and compared to a measured global methylated fraction. Alternatively, allele specific methylation fractions may be observed and compared to expected allele specific methylation fractions. This is possible when using a methylation detection technology that also includes genetic information such that reads can be phased with genetic loci. This is however not necessary in order to determine a difference between observed and expected methylation fractions, as global methylation fractions can easily be compared, which can be obtained using any methylation detection technology known in the art. The difference between the observed and expected methylated fractions may be indicative of the effect of an exposure factor / perturbation on the subject.

[0022] According to a fifth aspect, there is provided a method of determining the effect of one or more perturbations on the level of modified base at an epigenetic locus, the method comprising: receiving observed levels of modified base at the epigenetic locus in one or more samples that have been exposed to the one or more perturbations; receiving a previously identified statistically significant allele-specific relationship between the level of a modified base at the epigenetic locus and the genotype at a genetic locus, wherein the presence of said statistically significant allele-specific relationship indicates that the level of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy; receiving genotype data about the samples at the genetic locus; determining, for each sample, the expected level of modified base at the epigenetic locus on each chromosome copy of the subject using the genotype data and the identified allele-specific relationship; and comparing the observed and expected levels of modified base at the epigenetic locus in each sample, wherein the presence of a difference between the observed and expected level of modified base for a sample is indicative of an effect of the one or more perturbations to which the sample has been exposed on the level of modified base at the epigenetic locus.

[0023] The methods according to the present aspect may have any of the features described in relation to any preceding aspect. The level of modified base may be a methylated fraction. The method may further comprise obtaining the methylation fractions by determining the proportion of methylated bases at the epigenetic locus in the one or more samples. 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.

[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 first, second, third, fourth or fifth 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, second, third, fourth or fifth 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, second, third, fourth or fifth 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.

[0028] Summary of the Figures

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

[0030] Figure 1 is a flowchart illustrating a method of identifying a methylation quantitative trait locus (mQTL) according to a general embodiment of the disclosure.

[0031] Figure 2 is a flowchart illustrating a method of designing a genetic assay, providing a genetic assay and / or determining the likely modified base level (e.g. methylation fraction) at a locus in a sample, according to embodiments of the disclosure.

[0032] Figure 3 shows schematically a system for implementing methods of the disclosure. Figure 4 shows simulated methylation fraction data at a particular locus, processed using a method according to the disclosure where methylation fraction is determined on a single allele basis using phased genetic and epigenetic information (A), and using a method of the prior art where methylation fraction is determined on a genotype basis (B).

[0033] Figure 5 illustrates schematically a step of phasing heterozygote genotypes into haploblocks using sequencing data for samples in a cohort (A) and a step of counting methylated reads at CpG contexts associated with a haploblock allele (B).

[0034] Figure 6 shows an example of data that can be obtained as a result of applying the steps of Figure 5 on sequencing data for a plurality of samples (samples 1 to n) in a cohort.

[0035] Figure 7 illustrates schematically a step of aggregating data obtained using the steps of Figure 5 over samples and haploblocks compatible with each of the two alleles (referenced, variant=1) at a locus.

[0036] Figure 8 shows methylation fraction data for each of the two alleles at a locus quantified using a method of the disclosure, for 3 different loci. Each plot shows data (methylation fraction, i.e. haploblock modified C counts scaled by coverage at the CpG context), for a SNP. A. Known disruptive SNP at a CpG context (strong positive association between genotype and methylation fraction). B. New mQTL site that could only be identified with the methods of the disclosure, with high confidence negative association between genotype and methylation fraction. C. Site with high variance methylation fraction.

[0037] Figure 9 shows results of a statistical power analysis comparing the power for identifying mQTLs using methods of the disclosure and methods of the prior art, for strong effect sizes. The plots illustrate, using simulated data, power as a function of minor allele frequency, effect size, and number of samples, showing that the mQTL calling method using phased data consistently outperforms the one using unphased data. Thus, for a rare allele, a small effect size, a small number of samples, or a combination thereof, phasing allows for a larger power of detection of mQTLs. A. Statistical power of identifying the mQTLs as a function of the number of samples in the cohort, using the methods described herein based on individual alleles (blue data series at the top) and the conventional method based on genotypes (orange data series at the bottom), assuming a minor allele frequency of 0.5 in the study cohort. B. Statistical power of identifying the mQTLs as a function of the minor allele frequency in the cohort of samples analysed, assuming a cohort of 25 samples, using the methods described herein based on individual alleles (blue data series at the top) and the conventional method based on genotypes (orange data series at the bottom). C. Probability density function of methylation fraction for the effect sizes used in A and B.

[0038] Figure 10 shows results of a statistical power analysis comparing the power for identifying mQTLs using methods of the disclosure and methods of the prior art, for weaker effect sizes. A. Statistical power of identifying the mQTLs as a function of the number of samples in the cohort, using the methods described herein based on individual alleles (blue data series at the top) and the conventional method based on genotypes (orange data series at the bottom), assuming a minor allele frequency of 0.5 in the study cohort. B. Probability density function of methylation fraction for the effect sizes used in A.

[0039] Figure 11 shows the distribution of p-values obtained using a method of the disclosure applied to data from the Genome In A Bottle (GIAB) project, separated between 4 categories of SNP-CpG context pairs: disruptive mQTLs (SNPs overlapping with a highly methylated CpG and the disrupted CpG context, labelled “dis”), SNPs overlapping with a highly methylated CpG and a nearby CpG context (not the one that the SNP overlaps with, labelled “dis_snp”), disrupted CpG and SNPs linked to the disrupted CpG but not overlapping with it (labelled “dis_cpg”), and all other pairs (labelled “null”).

[0040] 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.

[0041] Detailed Description

[0042] 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.

[0043] 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 model fitting (requiring the optimisation of model parameters for a regression model between epigenetic modifications and allele level genotype, typically by maximum likelihood estimation using data from multiple samples and a parameter search strategy), is such that it cannot be performed in the human mind.

[0044] 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. 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.

[0045] A methylation quantitative trait locus (mQTL) is a relationship between a genetic locus and an epigenetic locus, where the genotype at the genetic locus is indicative of the level of modification at the epigenetic locus. The term is also used to refer to the genetic locus itself, i.e. a mQTL may refer to a genetic locus, the genotype of which is indicative of (i.e. associated with) the level of DNA modification at a particular locus. In the context of the present disclosure, a genetic locus refers to a genomic location at which there is sequence variation between individuals (i.e. different genotypes). The variation typically refers to a single nucleotide variant (SNV, in the somatic context) or single nucleotide polymorphism (SNP, in the germline context), although insertions / deletions (indels) can also be used in the context of the present disclosure (particularly small indels, such as indels that are less than 24 bases, or less than 12 bases in length). Throughout the present disclosure, the terms “locus”, “genetic locus”, “variant locus” and “SNP” will be used interchangeably to refer to loci at which there is variability (whether in the germline or somatic context). In the context of the present disclosure, an epigenetic locus typically refers to a genomic location at which there is a base modification in at least some individuals. An epigenetic may be a cytosine methylation (presence of a methylated or hydroxy methylated cytosine, which are the most common forms of base modification) or an adenosine methylation (e.g. presence of a N6-Methyladenosine, m6A). Thus, an epigenetic locus may refer to a site of cytosine modification. As used herein, cytosine modification encompasses one or both of methylated cytosines (also referred to as 5mC as the cystosine is 5- methylcytosine) and hydroxy methylated cytosines (also referred to as 5hmC as the cystosine is 5- hydroxymethylcytosine). Cytosine methylation typically occurs at CpG dinucleotides, also referred to as “CpG islands” and “CpG contexts”. Thus, the wording “epigenetic locus”, “CpG context” and “context” are used interchangeably herein to refer to sites of cytosine modification, although an epigenetic locus may also be a site of modification of another base such as adenosine. Thus, the term “mQTL” refers to an association between any genetic locus and any genetic locus. However, it is most commonly applied to the context of relationships between single nucleotide variants / polymorphisms and cytosine modification. Therefore, in the examples below the terms SNPs and modified I methylated cytosine may be used to illustrate the method, but the skilled person would understand that the same principles are applicable to other types of genetic and epigenetic variants.

[0046] Figure 1 is a flowchart illustrating a method of identifying an mQTL according to a general embodiment of the disclosure.

[0047] At step 1 10, sequence data is obtained for a plurality of samples. Each sample comprises genomic and / or mitochondrial DNA for a respective subject (also referred to interchangeably as “individual”). Thus, 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. For example, the sample may be a blood sample or sample derived therefrom (e.g. purified peripheral blood mononuclear cells), 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). Without wishing to be bound by theory, it is believed that many mQTLs are tissue specific and therefore samples of genomic and / or mitochondrial DNA from a single tissue or cell population is likely to be more informative than samples comprising DNA from multiple tissue types such as e.g. circulating DNA. The sample may be a sample of cells including a cell culture, and hence the terms “subject” and “individual” refer to samples comprising genomic material from a single identity rather than to organisms from which such cells may have been originally obtained. The plurality of samples or individuals from which the samples have been previously obtained may be referred to as “cohort” or “study cohort”. The data comprises sequence data obtained using any sequencing technology from which epigenetic and genetic 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 I modified or unmodified base). 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 , 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).

[0048] At step 112, the sequence data may be processed to identify genetic variants present in each of the plurality of samples. This process is commonly referred to as “variant calling”. This step is optional because the method may start from sequence data in which variants have already been called. Many variant calling methods exist in the art and any of these methods may be used. These include for example the GATK haplotype caller (DePristo et al. 201 1), Freebayes (Garrison and Marth, 2012), Samtools in combination with BCFtools (Li 2011), etc (reviewed in Kobolt, 2020). The output of the variant calling process may be recorded in a standard file format such as VCF (variant calling file) files. Individual files may be produced for each sample, and / or a combined file listing all variants identified may be obtained (e.g. joint VCF file, for example by merging individual VCF files using packages such as BCFtools, Li 201 1 , or by performing variant calling jointly). Advantageously, a joint VCF file may be obtained, describing the variants present in the cohort of samples to be analysed. A joint VCF file advantageously lists all genotypes of interest that will be analysed. A single sample variant call file generally only includes non-reference alleles. This means that a single sample VCF does not include any information at every SNP of potential interest to be analysed for which the particular sample does not have a non-reference allele confidently identified. By contrast, a joint VCF includes all genotype information at all SNPs of potential interest (i.e. genotype information for all samples at any SNP where at least one of the samples in the cohort comprises one or more non-reference alleles). This means in particular that the joint VCF includes information that distinguishes situations where a sample has a confident 0 / 0 genotype call (i.e. homozygous for the reference allele) for a SNP and situations where a sample has no evidence for the SNP (represented as - / - genotype, i.e. genotype missing). This information is not present in a single VCF and taking this information into account rather than assuming a 0 / 0 genotype results in higher accuracy of mQTL calls. Thus, step 112 may comprise obtaining, by analysing the sequence data for the plurality of samples in the cohort, a file comprising genotype information for each of a plurality of genomic locations to be analysed (such as e.g. each genomic location at which at least one of the samples has a genotype comprising a non-reference allele, i.e. a variant allele) and for each sample. The genotype information for a genomic location and sample may indicate whether the sample is: homozygous for the reference allele (0 / 0), homozygous for a non-reference (variant) allele (1 / 1), heterozygous (0 / 1) or genotype missing (- / -, i.e. a genotype at the location could not be inferred with confidence based on the sequence data obtained from the sample).

[0049] At step 114, the sequence data for each sample is phased. Phasing refers to separating the two alleles of a heterozygote into haplotypes. A haplotype is a set of genetic variants (e.g. SNPs) adjacent to one another on a chromosome which are likely to be inherited together. This step is optional because the method may instead start from sequence data that has already been phased. This may comprise a tag identifying variants that are in phase with each other. Phasing may use read-based phasing or population-based phasing. Read-based phasing uses mapped reads spanning at least two heterozygous variants to infer the phase. In other words, read-based phasing allows to reconstruct haplotypes of a sample purely from sequence reads. Multiple read-based phasing approaches are known in the art, and any of these approaches can be used. For example, commonly used read-based phasing approaches include WhatsHap (Martin et al. 2016), GATK’s HaplotypeCaller (Polin et al. 2017), HapCUT (Bansal et al. 2008) and phASER (Castel et al. 2016). Population based phasing methods (also referred to as “statistical phasing”) use large datasets to predict haplotypes based on statistical likelihood. These methods are based on the idea that shared ancestry and limited recombination give rise to shared haplotype blocks, such that it is possible to infer a maximum-likelihood model of haplotypes given a data set (either phasing a large cohort such as the 1000 genomes cohort, or phasing a single additional individual using a haplotype reference panel created previously from a large cohort). Multiple population-based phasing approaches are known in the art (as reviewed in Browning and Browning, 2011), and any of these approaches can be used. For example, commonly based approaches include MaCH (Li et al. 2010), IMPUTE2 (Howie et al. 2009), PHASE (Stephens and Sheet 2005), Eagle2 (Loh et al. 2017), SHAPEIT2 (O’Connell et al. 2014), SHAPEIT5 (Hofmeister et al. 2023) and BEAGLE (Browning and Browning, 2007). The use of populationbased phasing methods advantageously enables to identify longer haploblocks (thereby increasing the range of the associations between genetic and epigenetic loci that can be identified), but increases the probability of errors in the identification of haploblocks as those are not supported by direct physical evidence of the joint inheritance of two loci for each individual. Thus, conversely, the use of read-based phasing advantageously results in higher confidence mQTL calls, even though the range of relationships that can be discovered can be more limited. Note that this limitation also depends on the length of the reads available, and is less and less prominent as longer reads are achievable.

[0050] At step 116, one or more mQTLs (each corresponding to a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus) are identified using the sequence data, genetic variants and phasing information. This may comprise a step of 116A of identifying a plurality of haploblocks for each haplotype identified in each sample. Haploblocks are runs of heterozygotes that we believe to be in phase with each other in a sample. Each block is considered to represent one allele. Haploblocks are identified on an individual by individual (i.e. sample by sample) basis. Thus, the coordinates of the blocks and the loci that they encompass may differ between individuals. The fact that haploblocks may differ between individuals means that it is possible to differentiate between the effect of two SNPs that are on the same block in some individuals but not others. However, if two SPNs are always inherited together in the study cohort then it is not possible (with any computational method) to tell which SNP is associated with a change in methylation at any CpG. In other words, the two SNPs will be associated with the same effect even though one or both of them may be in fact associated with the methylation change. This is rare in human as the human genome has a deficit of rare alleles compared to other species. At such, it is rare for a whole cohort of human genomes to have identical haploblocks.

[0051] At step 116B, for each sample, the counts of (i) reads with a methylated C and (ii) reads with an unmethylated C at the epigenetic locus associated with each haploblock (i.e. each haploblock allele) are obtained. At this step, a count of modified and unmodified C for each haploblock for each sample (i.e. block x - allele H1 , block x - allele H2) is obtained. Modified C may encompass methylated C and / or hydroxymethylated C. Importantly, the counts are not at the level of SNPs but at the level of haploblocks. A count can be obtained for any epigenetic locus (such as e.g. CpG context) associated with an allele of the haploblock. This can be any epigenetic locus that overlaps with any read that also overlaps with a heterozygous position in the haploblock. Indeed, any such epigenetic locus can also be phased with the heterozygotes in the haploblock (i.e. it is possible to determine which allele of the haploblock a particular modified or non-modified read belongs to). Note that the terms “haploblock” and “haploblock allele” are used interchangeably. Thus, a haploblock can refer to a particular set of alleles for respective genetic loci that are in phase with each other in the particular sample or to a genomic block (range of genomic coordinates) in which all heterozygotes are phased (producing two haploblock alleles, where the genotype at each genetic locus in the genomic block is known, i.e. phased genotype information is available throughout the block). The alleles for respective genetic loci that are in phase with each other in the particular sample are necessarily heterozygotes in the sample, because phase is irrelevant for homozygous loci (since the genotype is the same on both phases).

[0052] At step 116C, the haploblock level data obtained for all samples is aggregated at the level of SNPs. For each SNP and each allele of the SNP (labelled as 0 or A for the reference allele and 1 or B for variant alleles), this step aggregates data from all samples using all counts for haploblocks that are consistent with the respective alleles of the SNP. Counts for haploblocks that are consistent with the respective alleles of the SNP can be: (i) when the SNP is heterozygous in a sample: counts associated with the haploblock allele that comprises the respective allele (i.e. either H1 or H2 counts), or (ii) when the SNP is homozygous in the sample: counts associated with either of the haploblock alleles (from reads that overlap the homozygous SNP position but could be phased with a particular allele of the haploblock because they also encompass a heterozygous position) and / or counts associated with degenerate reads (reads that overlap the homozygous SNP position but could not be phased with a particular allele of the haploblock because they do not encompass a heterozygous position). Each set of counts per haploblock and the corresponding allele at the SNP now forms a set of data that can be used to fit a relationship between the counts and the allele at the SNP (rather than the diploid genotype at the SNP). This aggregates the data into a two-state model (allele specific model, states 0 and 1) rather than the conventional 3 state model (genotype specific model, states 0 / 0, 0 / 1 and 1 / 1). As explained above, mQTLs are identifiable using the methods described herein when there are haploblocks in the cohort of samples that encompass both the genetic locus (e.g. SNP) and the epigenetic locus (e.g. CpG context), since the methods rely on phasing between genetic and epigenetic loci (i.e. counts at epigenetic loci that are on the same haploblock as a genetic locus). Thus, the power of identification of mQTLs may decrease with the distance between the genetic and epigenetic loci as fewer individuals in the discovery cohort will have haploblocks that encompass both loci. In other words, as the distance increase, the effective size of the population of samples from which a call can be made may decrease.

[0053] Aggregating the data may comprise, for each SNP, fitting a generalised linear model (e.g. a negative binomial regression model) to the number of modified C as a function of (at least) the total number of reads (sum of the number of modified Cs and the number of unmodified Cs) and the SNP allele. The model may further comprise one or more additional terms for other factors that are believed to be possible predictors of the methylation fraction at the epigenetic locus. This may include e.g. a term for a perturbation, treatment or phenotype (including e.g. exposure to a drug, stress, etc.), a term for an interaction between the total number of reads (also referred to as totalC or total count) and the allele and / or a term for an interaction between the allele and a perturbation / treatment / phenotype (i.e. condition). Negative binomial regression is a known method in the art and multiple implementations exist. For example, a negative binomial regression model can be fitted to counts data using the function glm.nb in the R package MASS (see cran.r- project.org / web / packages / MASS / index.html). Alternatively, the model used at this step may be a Poisson regression model, which is commonly used for count data. As another alternative, weighted least square regression may be used, where each methylation fraction is weighted based on the underlying coverage at the location in the sample. For example, weights for each observation may be determined using an empirical formula, such as e.g. w( ) = + Xqor a formula that penalises sharply coverages lower than a threshold and does not penalise coverages above another threshold may be used, such as e.g. w( ) = where x is the coverage (with xo=10, observations with coverage lower than 5x are sharply penalised and observations with higher than 20x are not penalised). The use of the negative binomial model is believed to be preferable as it accounts for the extra variability often observed in count data. Fitting a model may comprise determining the value of one or more parameters of the model including a parameter associated with the SNP allele, and determining a p-value associated with said parameters (i.e. a metric of statistical significance associated with each parameter). The p-value quantifies the statistical confidence that the value of the parameter is not 0, i.e. the statistical confidence in there being an association between the SNP allele and the counts of modified C. Step 116 may be performed for a plurality of genetic loci (e.g. SNPs) for any epigenetic locus, result in a plurality of p-values. These p-values may be corrected for multiple testing using any method known in the art, such as e.g. Benjamini-Hochberg or Bonferroni correction.

[0054] At step 116D, a statistical significance criterion (e.g. cutoff applied to the p-value associated with the model parameter associated with the SNP allele) is applied to the results of step 116C, where each genetic locus- epigenetic locus pair which satisfies the statistical significance criterion represents a mQTL. At step 1 18, the results of any one or more of the preceding steps (including e.g. the identity of a genetic locus identified at step 1 16D and / or the parameters of a corresponding model identified at step 116C) can optionally be provided to a user (e.g. through a user interface) or data store.

[0055] The methods described herein are applicable to any non-monoploid species in which epigenetic DNA modifications occur. The example implementations described refer to diploid genomes, but the same concept could be extended to genomes with higher ploidy. In embodiments, the subject is a subject with a diploid genome. Thus, a subject may be a eukaryote organism, including but not limited to vertebrates, and in particular mammalians such as a human, dog, cat, horse, or a model animal (such as a mouse, rat, etc.). In embodiments, a subject is a human subject, a model animal (e.g. mouse or rat), or a pet animal. In embodiments, the subject is a mouse. In embodiments, the subject is a human.

[0056] The methods described herein find application in a variety of contexts. For example, the methods described herein can be used to identify loci to be included in a genetic assay. 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. 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 has 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. 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. Variants at loci identified using the methods described herein may contribute to such risks by modifying the methylation fraction at disease relevant CpG loci.

[0057] Epigenetic variants have been found to be associated with a variety of diseases (reviewed in e.g. Villicana and Bell, 2021) including cancer, metabolic diseases such as type 2 diabetes, neu rod egene rative disease such as Alzheimer’s disease, Parkinson’s disease and multiple sclerosis, as well as complex phenotypes such as platelet function (linked to cardiovascular disease such as the risk of adverse events in patients with acute coronary syndrome) and fatty acid levels (linked to cardiovascular and metabolic diseases). Thus, the loci identified using methods described herein may be used to determine the risk of a subject developing a disease or disorder selected from: cancer, a metabolic disease, a neurodegenerative disease, or a cardiovascular disease.

[0058] Figure 2 is a flowchart illustrating a method of providing or designing a genetic assay according to embodiments of the disclosure, and a method of determining the expected base modification level at an epigenetic locus (e.g. methylation fraction at a CpG locus) in a subject and / or determining the effect of a perturbation on base modification level at an epigenetic locus (e.g. methylation fraction at a CpG locus.As explained above, CpG methylation is the most common and extensively studied base modification, and as such the description below will refer to methylation, CpG loci and methylation fraction for ease of reading. However, the methods are applicable to any base modification. At step 210, one or more mQTLs are identified as described by reference to Figure 1 . At optional step 212, one or more of these mQTLs may be selected, for example for disease relevance. For example, methylation at specific epigenetic loci has been shown to be associated with risks of developing certain diseases as explained above. The mQTLs identified (and optionally selected) can be used to design and optionally provide a genetic test, and to determine the expected level of base modification (e.g. methylated fraction) at epigenetic loci in these mQTLs (instead or in addition to obtaining an observed level of base modification - for example for the purpose of perturbation screening).

[0059] In relation to the former, at optional step 214, a plurality of genetic loci are selected. These include the genetic loci that are part of the selected mQTLs, and optionally also additional disease relevant genetic loci (disease-associated epigenetic and / or genetic loci). For example, in the context of a genetic test for cancer diagnosis or prognosis, genetic loci 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 216, a genetic testing panel that targets the selected genetic loci (including the genetic loci in the selected mQTLs as well as any additional optional genetic loci) is designed. This may comprise designing reagents for the targeted measurement of the panel of genetic loci (e.g. capture and / or detection probes, primers etc.). This may also comprise defining a polygenic risk score associated with one or more of the genetic loci that are measured in the genetic test. 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.). Defining a polygenic risk score can be done using any method known in the art, as the present disclosure is primarily concerned with the identification of genetic loci that are potentially informative (because they influence base modification at one or more epigenetic loci). Fitting a polygenic risk score for a particular phenotype of interest once a set of genetic loci has been identified is within the capability of the skilled person. The method may further comprise manufacturing a genetic test panel at step 218 (e.g. by producing reagents for the targeted measurement of the selected genetic loci and / or software for the calculation of a polygenic risk score based on measured genotypes at the selected genetic loci).

[0060] Instead or in addition to steps 214-218, mQTLs identified and optionally selected at steps 210-212 can be used to determine the expected level of base modification (e.g. methylated fraction) at an epigenetic locus in a subject, and to screen perturbations (e.g. drug candidates) for an effect on the level of modified base at an epigenetic locus. For example, at step 220, the genotype of one or more subjects at the genetic loci of the mQTLs identified (or selected) is obtained. This can be performed using any sequencing technology known in the art, or even targeted genotyping assays. At step 222, the genotype obtained at step 220 and the corresponding mQTL relationship from step 210 can be used to determine the expected level of methylation at any epigenetic locus for which an mQTL is known involving the obtained genotype. Thus, the expected methylation fraction of the subjects can be determined using genotype information without the need to perform epigenetic measurements. The expected methylation fraction may be allele specific or global. A global methylation fraction is a methylation fraction that is determined by combining information across all chromosome copies at the epigenetic locus. For example, in the case of a subject that is heterozygous, a different methylation fraction will be estimated for each of the chromosome copies of the epigenetic locus. A global expected methylation fraction can be derived from this by averaging. In the case of a subject that is homozygous, the global expected methylation fraction is the same as the allele specific one.

[0061] In embodiments where the method is used to determine the effect of a perturbation on the level of modified base at an epigenetic locus (such as e.g. for drug screening, genetic KO screening, expsorue to one or more physico-chemical or metabolic stresses etc.), an observed methylation fraction may also be obtained at step 224, this can be compared with the expected methylated fraction obtained at step 222. The presence of a difference between the observed and expected level of modified base for the sample is indicative of an effect of the one or more perturbations to which the sample has been exposed on the level of modified base at the epigenetic locus, and therefore the comparison can be used at step 224 to determine the effect of a perturbation applied to the sample on the level of modified base at the locus. Methods according to this embodiment may therefore comprise determining the proportion of methylated bases at the epigenetic locus in one or more samples that have been exposed to one or more perturbations to be screened. Further, methods according to this embodiment may further comprise exposing the one or more samples to the one or more perturbations.

[0062] Figure 3 shows an embodiment of a system for implementing methods of the disclosure, such as e.g. identifying one or more mQTLs, designing a genetic assay, etc. 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. The sequence data acquisition means 4 is configured to obtain sequence data from samples, in the form of DNA sequencing reads. 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), the PromethlON sequencerfrom Oxford Nanopore Technologies, orthe 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 fitted models between a genetic locus and an epigenetic locus, parameters (such as e.g. parameters of a fitted model between a genetic locus and an epigenetic locus, 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.

[0063] 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.

[0064] 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.

[0065] 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.

[0066] Any section headings used herein are for organizational purposes only and are not to be construed as limiting the subject matter described.

[0067] 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.

[0068] 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%.

[0069] Examples

[0070] EXAMPLE 1

[0071] Introduction

[0072] In the landscape of genetics and epigenetics, the discovery of associations between genetic variation and epigenetic modifications plays an important role in our understanding of the link between genetics and gene regulation. A substantial volume of research focuses on the detection of Methylation Quantitative Trait Loci (mQTLs), which are genetic variants linked to variations in DNA methylation patterns at specific CpG sites. In essence, mQTL calling relies on finding statistically significant associations between single nucleotide polymorphisms (SNPs) and methylated cytosines.

[0073] The examples below demonstrate the use of sequencing technology combining genetic and epigenetic information on a single read, to exploit phasing information and enhance mQTL detection. mQTL calls are most commonly identified by fitting a regression model of the methylation fraction vs genotype (see, e.g., review by Villicana & Bell 2021), as illustrated on Figure 4B. In heterozygotes, the methylation fraction should be determined by an even balance between alleles. However, in practice binomial sampling creates an imbalance of reads and a large spread of methylation fraction. Using the methods described herein, by phasing the contexts, it is possible to separate the two alleles of a heterozygote (i.e. identify which haplotype a particular methylated site is associated with) and gain statistical power to detect mQTLs by removing an important source of noise. This gives the approach a significant advantage compared to existing (unphased) mQTL discovery methods.

[0074] Methods

[0075] Data. The present examples use 2 technical replicates of 7 samples from the GIAB cohort (Zook et al. 2016). These are all extensively characterised human samples, available as cells or extracted DNA from the Coriell Institute for Medical Research. See www.nist.gov / programs-projects / faqs-genome-bottle. The data was aligned to GRCh38 using bwa_mem (Li, 2013), SNPs were called using GATK haplotype caller and phased using Whatshap.

[0076] Input. The methods used in these examples use as input: a joint VCF (variant call format) file describing the variants present in the cohort of samples to be analysed.

[0077] Obtaining haploblocks. The joint VCF file is then phased using a phasing algorithm. In practice, phasing means a tag is applied to heterozygotes in the VCF file that are within phase of each other. This is done per individual using read-backed phasing, which uses reads that overlap 2 heterozygous positions in order to phase heterozygotes (see below). A number of tools are available for both variant calling and phasing. In the present examples, the GATK haplotype caller (Polin et al. 2016) and Whatshap (Martin et al. 2016) were used for variant calling and phasing, respectively. Phasing refers to separating the two alleles of a heterozygote into haplotypes. Read-based phasing was used in these examples. Read-based phasing uses mapped reads spanning at least two heterozygous variants to infer the phase. Alternative phasing methods can be used such as population level phasing. Phasing can apply from between just two heterozygotes to the whole chromosome. An allele is an alternative base at a polymorphic site. The word is also used to refer to the set of alternative bases associated with a haplotype or haploblock (see below).

[0078] Based on the phasing information, all SNPs identified in a cohort of samples to be analysed are associated with a sample-specific haploblock. A haplotype is a set of DNA variants (e.g. SNPs) adjacent to one another on a chromosome which are likely to be inherited together. A haploblock is similar to a haplotype in that it refers to a set of variants that are adjacent to each other, but emphasising that we do not know the complete haplotype but a series of “haploblocks” made up of phased heterozygotes, but blocks are not necessarily in phase with one another. In other words, haploblocks are sample-specific sets of heterozygotes that we know to be in phase in a sample, but where we do not know whether the heterozygotes in one haploblock are in phase with the variants in another haploblock.

[0079] Identifying mQTLs. When haploblocks have been identified, alleles are defined for these, then the counts of CpG contexts associated with each allele of the haploblock are separately recorded, for each sample (i.e. for each sample a count per allele per haploblock is obtained). A CpG context associated with an allele of the haploblock is any CpG that overlaps with any read that also overlaps with a heterozygous position in the haploblock. Indeed, any such CpG can also be phased with the heterozygotes in the haploblock (i.e. it is possible to determine which allele of the haploblock a particular modified or non-modified read belongs to). Then, the counts per sample are aggregated using the following process: for each SNP of interest from the joint VCF file: for each sample: identify haploblock alleles: determine which haploblock alleles are consistent with the observed SNP allele(s) for each consistent haploblock: extract the total (totalC) of modified (modC) and unmodified Cs (totalC-modC) for each CpG context to obtain a table where each row is (modC, totalC, SNPallele) aggregate data across samples by performing a negative binomial regression on these data for each CpG context: modC~totalC+SNPallele

[0080] Negative binomial regression was performed using statsmodels in python: 0.14.0 (www.statsmodels.org / stable / index.html; Seabold and Perktold, 2010).

[0081] A SNP of interest was defined here as any SNP for which at least one variant allele is present in the cohort (and which can be associated with a haploblock). In other words, any SNP for which the cohort comprises at least two different genotypes in the cohort (i.e. there are at least two samples that have a different genotype selected from 0 / 0, 0 / 1 and 1 / 1) was included.

[0082] A model that uses modC as response variable rather than modC / totalC (i.e. the absolute number of reads with a modified C rather than the proportion of reads with a modified C) was used in order to avoid overdispersion. This is because the uncertainty associated with a proportion (a number between 0 and 1) is much larger than that associated with an absolute number (i.e. the uncertainty is higher when evaluating 0.5~SNPallele than when evaluating 10~20+SNPallele, even though both are based on the same data).

[0083] As explained further below, the regression model is calculated at the level of alleles (haplotype) not genotype, i.e. instead of prior art models where the methylation fraction is estimated as a function of states 0, 1 , 2 (homozygous ref allele, heterozygous, homozygous variant allele), here the model is fitted using a SNP variable that represents an individual allele (i.e. 0 or 1 , reference or variant) depending on the haploblock allele. A consistent haploblock is a haploblock allele that is consistent with a particular allele of the SNP.

[0084] The data used for this step is summarized at the sample level by 2 files: (i) a “block” file that contains information about the SNPs: specifically which haplotype block a particular allele of a SNP can be found in; and (ii) a “counts” file which contains information at the haploblock level, specifically detailing the counts of C's and mode's per haploblock allele.

[0085] Thus, through this process, CpG reads are assigned to a haplotype using phasing information. A read assignment can be uncertain (degen) if the read does not contain any heterozygous location.

[0086] An example of the content of a block file is provided in Table 1 .

[0087] Table 1. Example content of a block file. Note only a few rows are shown, in practice many more rows would be present.

[0088] In the example in Table 1 , two SNPs are observed at positions 51478 and 61986 on chromosome 1. For that specific sample, the first SNP at 51478 is homozygous. Only the T allele is present, designated with allele code 0. Given a T allele here, all of the haploblock alleles H1 , H2, or degen are consistent (i.e. H1 and H2 are both T at this position in this sample). H1 and H2 correspond to haplotype 1 and haplotype 2, respectively, with degen representing reads that are ambiguous. A read is ambiguous if it is not possible to assign it to a hapolotype, i.e. the read does not overlap with any heterozygous locations in the haploblock. While these reads cannot be assigned to a haploblock allele, they can still be useful to identify mQTLs. For example, in a cohort of 100 samples with 5 heterozygous samples (0 / 1) and 95 homozygous samples (0 / 0) at the location, the methylation on the 95 homozygous (0 / 0) samples is strong evidence for the background rate of methylation at 0 alleles. Thus, this information is useful even if for these 95 we are not doing any phasing. The second SNP at 61986 is an A / G heterozygote in this sample. For that SNP there is no degeneracy, so no need for the degen entry. When aggregating over allele A, it is necessary to use the haploblock allele H1 , and for allele G, H2.

[0089] An example of the content of a counts file is provided in Table 2.

[0090] Table 2. Example of a counts file. Note only a few rows are shown, in practice many more rows would be present.

[0091] The counts file provides the positions of C's within CpG contexts, along with their respective haploblocks (i.e. the haploblock they belong to). In the example counts file in Table 2 (associated with the block file in Table 1), the initial CpG at position 55328 is affiliated with the block chr1 :51478. Forthis particular sample, we know from the blocks file that all reads exhibit degeneracy. Consequently, the counts for mode and C's are exclusively attributed to the "degen" category, with no counts allocated to "H1" and "H2". These counts cannot be used for a SNP that is heterozygous in this sample, but are still useful for any SPN that is homozygous in the sample (as explained above). By contrast, the CpG at position 61953 is associated with the block chr1 :61986, which is heterozygous (het) in this sample. There is no “degen” count forthis block because there were no ambiguous reads in that sample at that location. In this scenario, we can resolve the CpG counts across the two alleles of the haplotype blocks (i.e. we can determine which haploblock allele each read is associated with). In other words, for any SNP, the counts used will be: for any sample in which the SNP is heterozygous: the counts for the respective allele (H1 or H2), and for any sample in which the SNP is homozygous: the counts for degen reads and H1 , H2 reads (if any).

[0092] Using the haplotype block ID we can then match each CpG to a SNP, and aggregate this data for each sample. This leads to a table like Table 3.

[0093] Table 3. Example of aggregated data from counts file for mQTL estimation.

[0094] In Table 3, we see that some samples appear only once, that’s because they are homozygous (either reference or alternative allele), and some samples appear twice, once per allele, because they are heterozygous. Using this data, we can split our CpG counts by their alleles, and ask the question: do we see different methylation levels per allele ? If we can answer “yes” with some degree of confidence, then we have a mQTL. Answering this question requires a statistical test, which in the present examples is a negative binomial regression. A negative binomial regression is a type of generalized linear model (GLM) used to analyse count data that exhibit overdispersion, meaning the variance is greater than the mean. It is an extension of the Poisson regression model, which is commonly used for count data, but the negative binomial model accounts for the extra variability often observed in count data. The counts data analysed in these examples tends to be overdispersed, which motivates the choice of this model. In the present examples, the model fits mode as a function of totC and the allele. Additional terms can be added to the model above. For example, terms that capture an interaction between the allele and a sample characteristic (e.g. sex or other phenotype, condition, etc.) can be included. In particular, conditions (e.g. drugs or KO of genes) known or postulated to influence methylation may be represented through an additional term and / or interaction term. An interaction term can capture the differential effect of the condition in the reference and variant genetic background at the SNP of interest.

[0095] Negative binomial regression and Poisson regression are used in prior art mQTL calling but applied to unphased data. Negative binomial regression is a GLM where the dependent variable Y is a count of the number of times an event (here modC) occurs, where P(Y=y) is given by the negative binomial distribution: where f is the gamma distribution, p>0 is the mean of Y and a>0 is a heterogeneity parameter to be estimated, and Inn = f>o + Pixi + / ?2 2- ^PPXPwhere xltx2, ...xpare predictor variables and

[0096] Po’Pi’Pf -Pparepopulation regression coefficients to be estimated, parameter estimation is typically performed using maximum likelihood estimation, i.e. identifying the a and p values that maximise the likelihood of the observed data. This also enables estimation of a variance-covariance matrix of the estimators (equal to -FT1where H is the Hessian matrix of second derivatives of the log likelihood function), which in turn can be used to estimate p-values for the coefficient estimates.

[0097] In the present context, three p parameters are estimated (intercept, coefficient for the total C, coefficient for the SNP allele), and a significant association exists between the number of modC at a CpG site when the parameter for the SNP allele is significant at a chosen level of confidence (e.g. p<0.05).

[0098] Power analysis. As previously discussed, the advantage of using reads including both genetic and epigenetic information (e.g. +modC) to call mQTLs is that we are able to split heterozygotes into their constituent alleles and remove the noise associated with binomial sampling. As a proof of concept to demonstrate the gain in detection power induced by phasing the data, the inventors generated some simulated data of CpG methylation fraction as follows. Consider N samples, where each sample has a given genotype: 0 / 0, 0 / 1 , or 1 / 1. Balance between the three genotypes is given by the Hardy-Weinberg equilibrium. We assume each sample has 15X mean coverage following a Poisson distribution. Coverage is split into two alleles with probability following a binomial distribution. For each allele, modC are generated from a beta-binomial of the coverage. For each allele, we now have modC and coverage, so we can compute a methylation fraction. Figure 4 illustrates, using this simulated data, what an mQTL candidate would look like when being evaluated using phased data (split into two alleles, Fig. 4A) and unphased data (split into three genotypes, Fig. 4B). We can then perform a power analysis, where power is interpreted as sensitivity. For the unphased model, we use a negative binomial regression to fit mode as a function of totalC and the three genotypes 0 / 0, 0 / 1 and 1 / 1 . For the phased data, we have split the heterozygotes into their alleles, and we fit mode as a function of totalC and the allele code (0 or 1). The results are shown on Figure 9 and demonstrate that a model using phased data performs better than a model using unphased data.

[0099] Results

[0100] Current processes for identifying mQTLs (SNP loci where the allele present at the locus affects the methylation status at a CpG context nearby) typically rely on large cohorts of samples, where for each SNP, subjects are classified into 3 groups based on their genotype at the SNP (e.g. AA, AB, BB, where A is the reference allele and B the variant allele), and the methylation fraction at candidate CpG contexts are obtained for all individuals in each category, then the a linear model is fitted to determine whether there is a statistically significant associated between the genotype groups at the SNP and the methylation fraction at a candidate CpG context. The present inventors identified that a problem with this approach is that for individuals with an AB genotype, as CpG methylation is probed separately from the genotyping, it is not possible to determine which fraction of the methylated I non-methylated reads came from the chromosome copy with allele A and which came from the chromosome copy with allele B. This would not be a problem if the sequencing process resulted in equal proportions of reads from the two copies. However, sequencing samples the DNA molecules present in the sample and therefore by chance alone it is possible to get significantly more reads from one of the two alleles than from the other. This introduces significant amounts of noise in the data available (see Figure 4B, where the heterozygous genotype is shown to be associated with a greater spread of methylation fraction than the two homozygous genotypes). Discarding data from all heterozygotes at a particular SNP removes the noise associated with heterozygotes but also discards the information associated with these individuals, essentially reducing the size of the cohort, increasing the variance of the regression and reducing the statistical power of identification of mQTLs.

[0101] The present inventors recognised that another solution would be preferable, which would also remove this source of noise but would not reduce the statistical power of identification of mQTLs. In other words, the proposed solution increases the statistical power of identification of mQTLs given a cohort of a given size, or conversely enables identification of mQTLs with a given statistical power using fewer samples than was previously possible. The proposed solution phases the CpG contexts and SNPs to be analysed. This enables separation of the two alleles of each heterozygote, and fitting of a model of methylation fraction vs genotype at the allele level (i.e. considering individual allele genotypes rather than diploid genotypes, resulting in a 2 states model rather than a 3 states model), as illustrated on Figure 4A.

[0102] The method works by phasing heterozygote genotypes into haploblocks using short read information for the samples in the cohort. This is illustrated on Figure 5A, where the left column illustrates schematically a set of reads at heterozygous loci in a sample, and the right column illustrates two haploblocks (one at the top, one at the bottom), each associated with a haploblock allele, and comprising a respective set of SNPs that are inherited together. The haploblocks are therefore each associated with a haploblock allele. For example, for haploblock 1 the alleles are: [SNP1=1 , SNP2=1 , SNP3=0, SNP4=0] and [SNP1 =0, SNP2=0, SNP3=1 , SNP4=1], Having obtained haploblocks, the count of methylated reads at CpG contexts associated with each haploblock allele are obtained individually for each sample. This is illustrated on Figure 5B, where the reads aligning to allele 2 of haploblock 2 are shown, with CpG contexts associated with the haploblock highlighted as dots and shown in a different colour depending on whether the particular CpG in the read mapping to the particular haploblock allele was methylated or not. This can be repeated for each sample in the cohort, leading to data as illustrated on Figure 6, where for each haploblock in each sample, and each CpG context associated with each haploblock, the number of reads with modified and unmodified C are counted that are consistent with haploblock allele 1 or haploblock allele 2. Note that haplotype labels (haploblock allele 1 - H1 , or haploblock allele 2 - H2) have no meaning across blocks (i.e. H1 and H2 are simply a first and second haplotype for a particular block, with no particular relationship to either H1 or H2 of another haploblock), and haploblocks have no meaning across samples (i.e. blocks are identified individually for each sample - this is also illustrated on Figure 7, where the haploblocks that are associated with alleles 0 and 1 of a SNP all overlap the position of the SNP but may otherwise differ between the samples in the cohort).

[0103] Finally, the data is aggregated over samples and haploblocks to obtain all information to make a call. This is illustrated on Figure 7, where for a given SNP to be analysed, the counts data from all haploblocks alleles across samples that are consistent with the reference allele (left) are compared with the counts data for all haploblock alleles across samples that are consistent with the alternative (variant) allele (right). Note that as illustrated on Figure 7, there may be imbalanced numbers of haploblocks alleles associated with each SNP allele. This is because in the case of homozygote samples two haploblock alleles from a sample appear on the same side (i.e. providing data for the reference allele of the SNP), whereas with heterozygotes the two haploblock alleles appear on opposite sides (providing data for the reference and alternative alleles of the SNP). A call per SNP and CpG is then made using negative binomial regression as explained above. This is illustrated for real data on Figure 8 for selected SNP-CpG pairs, in the Genome In a Bottle (GIAB) data. Figure 8A shows data (methylation fraction, i.e. haploblock mode counts scaled by coverage at the CpG context, i.e. [C+modC]=[totalC]) for a known disruptive SNP at a CpG context (identified with -login p-value=8.0). Figure 8B shows data for a potential mQTL site that has never been previously identified, but could be identified using the methods described herein with high confidence (- logio p-value=7.4). Figure 8C shows data for a site where the methylation fraction at the CpG has high variance, and it is therefore not possible to identify a relationship between the SNP and the CpG context with confidence (-logio p-value=5.8). This shows that the method was able to identify with confidence known as well as novel mQTLs supported by allele level data, and also to discriminate sites where the association between alleles and CpG count is not significant.

[0104] The increased power of the proposed method was verified using a power analysis. The power of a method of identifying mQTLs is the probability that an effect is detected using the method given that the effect exists (i.e. P(effect detected | there is an effect). For any method this depends on the strength of the effect and the number of samples available. Higher number of samples are required to be able to identify a first mQTL having a smaller effect size than a second mQTL at a chosen level of significance. The present inventors hypothesised that by removing an important source of noise in conventional mQTL identification, the new methods would have a higher power at any given number of samples at which the power is not 1 for conventional methods. Simulated data was generated (10,000 simulations per data point) to investigate the effect of using the phased epigenetic and genetic information model on the power of detection of associations between SNPs and CpG contexts. This is done by specifying beta parameters from which the true fractions are sampled for each simulation. Figure 9 shows the results of this analysis for a hypothetical mQTL with strong effect size (i.e. a strong association between genotype and methylation fraction - in particular, with true fraction of methylation for allele A of 50 / 60 (10 unmethylated, 50 methylated) or 10 / 60 (50 unmethylated, 10 methylated)). Figure 9A shows the statistical power (proportion of all simulations that have a p-value below alpha=0.05, over all 10,000 simulations with the indicated parameters) of identifying the mQTLs with the stated effect size (beta parameters specifying the probability density function of methylation fraction illustrated on Figure 9C) as a function of the number of samples in the cohort, using the methods described herein based on individual alleles (blue data series at the top) and the conventional method based on genotypes (orange data series at the bottom), assuming a minor allele frequency of 0.5 (i.e. the two alleles are equally represented in the cohort). Figure 9B shows the statistical power (as above) of identifying the mQTLs with the stated effect size as a function of the minor allele frequency in the cohort of samples analysed, assuming a cohort of 25 samples, using the methods described herein based on individual alleles (blue data series at the top) and the conventional method based on genotypes (orange data series at the bottom). These data show that the methods proposed here have a higher statistical power for identifying mQTLs compared to the prior art over a range of number of samples (the exact number depending on the strength of effect size of the mQTL), and over all minor allele frequencies for cohorts within this range of samples. Figures 10A and 10B show data corresponding to Figures 9A and 9C, for a hypothetical mQTL with weaker effect size (i.e. a weaker association between genotype and methylation fraction), assuming a minor allele fraction of 0.25. This shows that the prior art methods do not even reach the levels of power achievable with the new methods at numbers of samples as high as 175. In other words, the weaker the QTLs the larger the range of number of samples at which the present methods have higher statistical power. In other words, the higher the mQTL that larger the difference between the number of samples that would be required to identify the mQTL using methods of the prior art vs using the methods described here.

[0105] This analysis demonstrates that the samples requirements for identifying mQTLs with a given strength of effect are reduced compared to the prior art, and that given a number of samples available, mQTLs are identified with greater sensitivity since smaller effects can be detected. In other words, given a number of samples available, the methods described herein are able to identify mQTLs that previously could not be identified. Note that in the implementation used in these examples the analysis is limited to CpG and SNP pairs that are within a phase-able distance of each other (i.e. where the data contains enough information to determine the phase of the CpG and SNP). This is a few hundred base pairs, the exact number depending on the location, as variable length stretches of the genome can be phased, in the present examples. However, much larger distances can be analysed using statistical phasing using reference panels, for example using the SHAPEIT5 algorithm as described in Hofmeister et al. 2023.

[0106] The methods proposed were additionally validated by investigating the results obtained for disruptive mQTLs (disQTLs). Disruptive mQTLs are the simplest example of mQTL. They occur when a SNP overlaps with a highly methylated CpG, disrupting the context. A method that correctly identified mQTLs would be expected to produce significant p-values at these sites (labelled “dis”). However, where these SNPs are located close to another CpG context, we would not necessarily expected significant p-values (labelled “dis_snp”). Further, the disrupted CpG context may be linked to other SNPs (labelled “dis_cpg”). We expect these to be significant but with lower power than 'dis'. Anything that does not fall into these 3 categories was labelled as “null”. The distribution of p-values for each of these 4 categories in data from the Genome In a Bottle (GIAB) project (using data from all 7 samples, phased using Whatshap as explained above) is shown on Figure 11 . This shows that the distribution of p-values for the dis_snp and null categories were identical and included mostly non-significant p-values, whereas the distributions of p-values for the dis and dis_cpg categories mostly included significant p-values (the former comprising more density towards the lower / more significant p-values than the latter). This further validates the methods described on real, widely used benchmark data.

[0107] In addition to the extension already mentioned above which considers interaction terms, the model above (modc~totalC+genotype, where genotype is an allele level genotype) could be extended to support other predictor variables as well as their interaction with genotype. For example, a model could consider the presence of a particular phenotype, such as sex, presence of an epigenetic disorder, treatment or manipulation that disrupts an epigenetic process such as an epigenetic enzyme knockout. This can be included in the model as an independent term such as modC~totalC+genotype+phenotype or as an independent term and an interaction term modC~totalC+genotype+phenotype+(genotype*phenotype)). Phenotypes are associated with samples and can easily be included as a predictor variable. This can allow the method to be used to answer questions like: are there any CpG contexts that behave differently given genotype and phenotype, and do mutations in methylation genes (e.g. TET1 or DNMT3A) lead to a different QTL relationship?

[0108] References

[0109] 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.

[0110] Ahmed, M. et al. CRISPRi screens reveal a DNA methylation-mediated 3D genome dependent causal mechanism in prostate cancer. Nat. Commun. 12, 1781 (2021).

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

[0112] 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).

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

[0114] Bansal, V., Bafna, V.: HapCUT: an efficient and accurate algorithm for the haplotype assembly problem. Bioinformatics 24(16), 153-159 (2008). DePristo, M.A., Banks, E., Poplin, R., Garimella, K.V., Maguire, J.R., Hartl, C., Philippakis, A.A., Angel, G.D., Rivas, M.A., Hanna, M., McKenna, A., Fennell, T.J., Kernytsky, A.M., Sivachenko, A.Y., Cibulskis, K., Gabriel, S.B., Altshuler, D., Daly, M.J.: A framework for variation discovery and genotyping using nextgeneration DNA sequencing data. Nature Genetics 43(5), 491-498 (2011).

[0115] Castel, S.E., Mohammadi, P., Chung, W.K., Shen, Y., Lappalainen, T.: Rare variant phasing and haplotypic expression from RNA sequencing with phASER. Nature Communications 7, 12817 (2016).

[0116] Marcel Martin, Murray Patterson, Shilpa Garg, Sarah O. Fischer, Nadia Pisanti, Gunnar W. Klau, Alexander Schoenhuth, Tobias Marschall. WhatsHap: fast and accurate read-based phasing bioRxiv 085050 doi: 10.1101 / 085050

[0117] Li Y, Wilier CJ, Ding J, Scheet P, Abecasis GR. MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes. Genet Epidemiol. 2010;34:816-34.

[0118] Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am J Hum Genet. 2005;76:449-62.

[0119] Howie BN, Donnelly P, Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet. 2009;5:e1000529.

[0120] Browning SR, Browning BL. Rapid and accurate haplotype phasing and missing-data inference for wholegenome association studies by use of localized haplotype clustering. Am J Hum Genet. 2007;81 :1084- 97.

[0121] Browning SR, Browning BL. Haplotype phasing: existing methods and new developments. Nat Rev Genet. 201 1 Sep 16;12(10):703-14.

[0122] Simpson, J.T., Workman, R.E., Zuzarte, P.C., David, M., Dursi, L.J. and Timp, W. (2017) Detecting DNA cytosine methylation using nanopore sequencing. Nat. Methods, 14, 407-410.

[0123] Koboldt, D.C. Best practices for variant calling in clinical sequencing. Genome Med 12, 91 (2020).

[0124] Erik Garrison, Gabor Marth. Haplotype-based variant detection from short-read sequencing. 20 Jul 2012 arXiv:1207.3907 Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 ;27(21):2987- 93.

[0125] Ryan Poplin, Valentin Ruano-Rubio, Mark A. DePristo, Tim J. Fennell, Mauricio O. Carneiro, Geraldine A. Van der Auwera, David E. Kling, Laura D. Gauthier, Ami Levy-Moonshine, David Roazen, Khalid Shakir, Joel Thibault, Sheila Chandran, Chris Whelan, Monkol Lek, Stacey Gabriel, Mark J. Daly, Ben Neale, Daniel G. MacArthur, Eric Banks. Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv 201178; doi: https: / / doi.org / 10.1101 / 201178

[0126] O'Connell J, Gurdasani D, Delaneau O, Pirastu N, Ulivi S, Cocca M, Traglia M, Huang J, Huffman JE, Rudan I, McQuillan R, Fraser RM, Campbell H, Polasek O, Asiki G, Ekoru K, Hayward C, Wright AF, Vitart V, Navarro P, Zagury JF, Wilson JF, Toniolo D, Gasparini P, Soranzo N, Sandhu MS, Marchini J. A general approach for haplotype phasing across the full spectrum of relatedness. PLoS Genet. 2014 Apr 17;10(4):e1004234.

[0127] Loh PR, Danecek P, Palamara PF, Fuchsberger C, A Reshef Y, K Finucane H, Schoenherr S, Forer L, McCarthy S, Abecasis GR, Durbin R, L Price A. Reference-based phasing using the Haplotype Reference Consortium panel. Nat Genet. 2016 Nov;48(11):1443-1448.

[0128] Zook, J., Catoe, D., McDaniel, J. et al. Extensive sequencing of seven human genomes to characterize benchmark reference materials. Sci Data 3, 160025 (2016). Seabold, Skipper, and Josef Perktold. “statsmodels: Econometric and statistical modeling with python.” Proceedings of the 9th Python in Science Conference. 2010.

[0129] Li H. (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2 [q-bio.GN], Hofmeister RJ, Ribeiro DM, Rubinacci S, Delaneau O. Accurate rare variant phasing of whole-genome and whole-exome sequencing data in the UK Biobank. Nat Genet. 2023 Jul;55(7):1243-1249.

[0130] 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 identifying a methylation quantitative trait locus, the method comprising: receiving sequence data comprising sequence reads and / or information derived therefrom including information about the presence and location of one or more genetic variants and modified bases in the sequence reads, for a plurality of samples; and identifying, using said sequence data, a statistically significant relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus, thereby identifying a methylation quantitative trait locus, wherein the relationship is allele-specific and the presence of said statistically significant allele-specific relationship indicates that the fraction of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy.

2. The method of claim 1 , wherein the genetic variant is a single nucleotide variant or polymorphism, the modified base is a modified cytosine and / or the epigenetic locus is a CpG locus.

3. The method of any preceding claim, wherein the step of identifying a statistically significant allelespecific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus comprises, for each sample individually: identifying a plurality of haploblocks, wherein each haploblock corresponds to a particular set of alleles for respective genetic loci that are in phase with each other in the particular sample; and obtaining, for each haploblock, a modified read count and an unmodified read count, wherein the modified read count is a count of the number of the sequence reads that align to and are consistent with the haploblock and that show a modified base at an epigenetic locus, and the unmodified read count is the count of the number of the sequence reads that align to and are consistent with the haploblock and that show an unmodified base at the epigenetic locus, optionally wherein said counts are obtained individually for each epigenetic locus associated with a haploblock.

4. The method of claim 3, wherein the step of identifying a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus further comprises: for each of one or more variable genetic loci identified in the plurality of samples, fitting a relationship between the modified read counts, the unmodified read counts and the allele at the variable genetic locus for any haploblock comprising the variable genetic locus.

5. The method of any preceding claim, wherein identifying a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus comprises fitting a regression model where the response variable is the count of reads with a modified base at the epigenetic locus and the predictive variables include the total number of sequence reads overlapping the epigenetic locus and the allele at the genetic locus on a single chromosome copy, wherein only sequence reads that map to the same chromosome copy are counted.

6. The method of claim 5, wherein the allele at the genetic locus on a single chromosome copy is the allele present at the genetic locus on a haploblock encompassing the genetic locus, and the reads that map to the same chromosome copy are reads associated with the haploblock.

7. The method of claim 5 or claim 6, wherein the regression model is a negative binomial regression model, and / or wherein reads that map to the same chromosome copy as the allele and / or reads that are associated with the haploblock are reads that overlap the epigenetic locus and at least one heterozygous genetic locus comprised in the haploblock.

8. The method of any of claims 5 to 7, wherein the regression model is fitted using data comprising, for each haploblock encompassing the genetic locus in each sample: counts of reads associated with the haploblock in the sample that show the presence of a modified base at the epigenetic locus, counts of reads associated with the haploblock in the sample that do not show the presence of a modified base at the epigenetic locus, and the allele present at the genetic locus in the haploblock in the sample.

9. The method of any of claims 5 to 8, wherein the regression model is a model of the form modC~totalC+SNPallele+x where mode the count of reads with a modified base at the epigenetic locus associated with a single chromosome copy in each sample, total C is the count of reads at the epigenetic locus associated with the chromosome copy, and SNPalllele is the allele at the genetic locus on the chromosome copy, and x represents one or more optional additional variables.

10. The method of any of claims 5 to 9, wherein fitting the regression model comprises identifying a parameter associated with each predictive variable and a metric of statistical significance associated with each parameter, and wherein a statistically significant allele-specific relationship between the fraction of a modified base at an epigenetic locus and the genotype at a genetic locus is identified when the metric of statistical significance associated with the variable that represents the allele at the genetic locus satisfies one or more predetermined criteria.11 . The method of any of claims 5 to 9, wherein the method comprises fitting the regression model for each of a plurality of candidate pairs of epigenetic and genetic loci, and identifying one or more mQTLs by selecting those pairs of epigenetic loci and genetic loci for which the regression model indicates a statistically significant allele-specific relationship between the fraction of a modified base at the epigenetic locus and the genotype at the genetic locus.

12. 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 5 or 6-letters code.

13. A method of identifying a methylation quantitative trait locus (mQTL), the method comprising: obtaining sequence data comprising sequence reads for a plurality of samples by sequencing genetic material in said samples; andanalysing said sequence data using the computer-implemented method of any of claims 1 to 12, optionally wherein the step of obtaining sequence data comprises using a sequencing technology that provides an output from which epigenetic and genetic bases can be called on the same read.

14. A computer-implemented method of designing a genetic test, the method comprising: identifying one or more methylation quantitative trait loci using the method of any of claims 1 to 12, and including the one or more genetic loci of the identified methylation quantitative trait loci in a panel of genetic loci that are measured in the genetic test.

15. The method of claim 14, wherein the panel of genetic loci comprises the one or more genetic loci of the identified mQTLs and one or more disease-associated genetic loci, and / or wherein the one or more epigenetic loci of the identified mQTLs are disease-associated epigenetic loci, optionally wherein the disease-associated genetic loci and / or disease-associated epigenetic loci are loci that are associated with the risk of a subject having a disease or disorder, and / or are predictive of a prognosis of diagnosis in a subject.

16. The method of any of claims 14 or 15, further comprising defining a polygenic risk score associated with one or more of the genetic loci that are measured in the genetic test.

17. The method of any of claims 14 to 16, further comprising designing reagents for targeted measurement of the panel of genetic loci.

18. A method of determining the expected level of base modification at an epigenetic locus in a subject, the method comprising: receiving a previously identified statistically significant allele-specific relationship between the fraction of a modified base at the epigenetic locus and the genotype at a genetic locus, wherein the presence of said statistically significant allele-specific relationship indicates that the fraction of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy; receiving genotype data about the subject at the genetic locus; and determining the expected level of base modification at the epigenetic locus on each chromosome copy of the subject using the genotype data and the identified allele-specific relationship.

19. The method of claim 18, wherein the previously identified statistically significant allele-specific relationship has been identified using the method of any of claims 1 to 12, or wherein the method comprises identifying the statistically significant allele-specific relationship using the method of any of claims 1 to 12.

20. The method of claim 18 or claim 19, further comprising receiving an observed methylated fraction at the epigenetic locus, and comparing the observed methylated fraction with the expected methylated fraction, optionally wherein the observed and expected methylated fractions are global methylated fractions at the epigenetic locus.21 . A method of determining the effect of one or more perturbations on the level of modified base at an epigenetic locus, the method comprising: receiving observed levels of modified base at the epigenetic locus in one or more samples that have been exposed to the one or more perturbations; receiving a previously identified statistically significant allele-specific relationship between the level of a modified base at the epigenetic locus and the genotype at a genetic locus, wherein the presence of said statistically significant allele-specific relationship indicates that the level of the modified base expected at the epigenetic locus depends on the presence or absence of a genetic variant at the genetic locus on the same chromosome copy; receiving genotype data about the samples at the genetic locus; determining, for each sample, the expected level of modified base at the epigenetic locus on each chromosome copy of the subject using the genotype data and the identified allele-specific relationship; and comparing the observed and expected levels of modified base at the epigenetic locus in each sample, wherein the presence of a difference between the observed and expected level of modified base for a sample is indicative of an effect of the one or more perturbations to which the sample has been exposed on the level of modified base at the epigenetic locus.

22. The method of claim 21 , wherein the level of modified base is a methylated fraction and / or wherein the method further comprises obtaining the methylation fractions by determining the proportion of methylated bases at the epigenetic locus in the one or more samples and / or wherein the method further comprises exposing one or more samples to one or more perturbations.

23. The method of claim 21 or claim 22, wherein the one or more perturbations are selected from: exposure to a drug, exposure to a physico-chemical stress, and exposure to a metabolic stress.

24. 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 23.

25. 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 23.

Images