Scientist Stories: Kasper Hansen, Computational analysis of DNA methylation in Disease

My name is Kasper Daniel Hansen, and I'm welcoming you to this talk on computational analysis of DNA methylation data. First, we're going to talk about how to measure DNA methylation. In mammals, DNA methylation occurs exclusively at CpG dinucleotides, with the exception of neurons and stem cells that have some amount of non-CpG methylation. CpGs occur in specific patterns along the genome. We sometimes find them clustered together in regions of high of high degree of CpG occurrences called CpG islands. Humans have about 28 million CpGs, and they're depleted genome-wide because DNA methylation in the germline increases the mutation rate of CpGs. The most popular way these days to measure DNA methylation is by sulfite conversion coupled with some kind of readout, either microarray or sequencing. Historically, there was an abundance of assays proposed to measure DNA methylation, and there was substantial disagreement about which assays were best. There are two important exceptions to bisulfite conversion that are sort of important use today. First, third-generation sequencing, specifically nanopore sequencing but also pipe biosequencing, allows the detection of DNA modifications without bisulfite conversion. This typically requires a custom base caller. In the future, it's pretty clear that all Nanopore sequencing, all Nanopore DNA sequencing, will result in DNA methylation calls being provided to the end user. And then there's an emerging and important field of cell-free DNA, very important in cancer diagnostics and surveillance, where they use a different set of assays including MeTip, which is an immunoprecipitation-based assay, And they do this because they have a requirement of being able to handle very low amounts of input material. Now we're mostly going to focus on bisulfide— we're almost exclusively going to focus on bisulfide-converted assays in this talk here. And I have to mention that bisulfide conversion cannot distinguish methylation from hydroxymethylation. Most tissues, but not all, have very low levels of hydroxymethylation. That doesn't mean it doesn't matter, but, it means that it's relatively hard to detect. The main bisulfite platforms available today are the Illumina methylation microarrays. Historically, we had the 450K array, and currently there's the EPIC array. And then there's bisulfite sequencing, which can be either whole genome bisulfite sequencing, which is sort of the most expensive and most comprehensive way of measuring DNA methylation, at least in terms of the number of CPGs. And then there's two methods for sort of basically sequencing a smaller part of the genome, either capture sequencing where you use some kind of method to enrich for specific regions along the genome, or, and also a method called reduced representation by sulfide sequencing where you use kind of a chemical trick to enrich for a specific fixed part of the genome. So one disadvantage with reduced representation by sulfide sequencing is you can't just say I want to measure these genes you're basically given by the assay what you're seeing. It's important when you analyze bisulfide data to be aware that technical aspects of the assay have a large impact on data quality. For example, the type of polymerase you use, which can really dramatically change the data. There's actually a very nice paper in Frontiers in Genetics on this, and there was a recent review of library preparation protocols for whole genome bisulfite sequencing. In bisulfite treatment, it's a chemical treatment that converts unmethylated Cs to Gs. So along the genome, we have Cs here. Some are unmethylated and some have a methyl group attached. And after bisulfite treatment, Cs that are unmethylated get converted into Gs, and unmethylated— methylated Cs are untouched. Now remember that most, in mammals, most Cs are not in a CpG context and therefore unmethylated. So this essentially changes the genome into a genome with 3 bases. This converts the chemical information carried by the chemical modification into a DNA difference, which can now be assayed by sequencing or other methods. Here's an example with sequencing. I've run it through my sequencer, I've generated 8 reads. I've aligned them to the genome, and then I compare the sequence of the read over the C, and I count. In this case, I have 5 reads where they have a T, I have 3 reads that has a C, and I conclude that this position is 3 divided by 8% methylated. Note here that we have 2 measurements from the same sample. We can call them M and U, and here you can think of it as coverage. This is the number 8, and M is the number 3, and we end up with sort of percent methylation by forming a ratio of these two numbers. This is the same situation with methylation arrays, and in my opinion, it's one of the reasons for why DNA methylation data is relatively— by sulfite converter data— is relatively well behaved, because we can think of these two measurements as kind of an internal control, control internal to the sample. So let's look at some sequence data and discuss smoothing. Here I've displayed data from the first sample of the first paper where they assayed a human sample with whole genome bisulfite sequencing. At the top, we have individual— we have the percent methylation for individual CpGs. Each CpG is a dot, and this is the number 3 divided by 8 that we saw in previous slides. At the bottom, I've depicted the coverage, which is the number 8 depicted earlier, and tells us like how deeply we've sequenced each single CpG. What we see here immediately is that there's a high degree of correlation in the plot. There's a high degree of spatial structure. We have regions of high methylation and regions of low methylation and sort of some transition between those two states. When I see data like this, I always want to smooth it, and I would say, I would argue that at least for this particular genomic region, the black line is a quite faithful representation of the methylation level of individual CpGs, perhaps with the exception of the region out here where there's a little higher sort of amount of variability around the line. This is important because in sequencing we pay money to increase the coverage, increase the height of the blue points below, and in some way what we're doing is we're trying to reduce technical variation or sequencing variation by sequencing deeper. In this paper, they actually generated data from 6 individual library preparations of the same sample. And here I've taken one of the 6 library preparations, and I've plotted the data here. We can see it's much lower coverage than previously, and we can see there's sort of some spatial structure that looks a little bit similar to before, but there's a lot more sort of variation, or it doesn't look quite as convincing as before. We can smooth the data using a method that I developed in 2012. We get this pink line here, and we can compare the high coverage and the low coverage data and see that there's basically, in this particular region, no difference between the two different approaches, suggesting that we can make do with perhaps lower coverage and use that money we've saved on sequencing the sample lower to generating more replicates, which is something I am a big believer in. Yeah, so this is a technique that I like to use for sequencing data. Yeah, so let's talk a little bit about the structure of DNA methylation because it's pretty important to know what do things look like when we study them. Here I've taken samples from liver, spleen, and brain regions at the same individuals. Who contribute different samples here. And I have smoothed my methylation data with a very large window size. So the lines that are depicted on this plot here do not represent the methylation level in individual CpGs. Rather, they represent the methylation level of sort of bigger regions of methylation. What we have depicted here is a 9-megabase region across the genome.. And we see this sort of very distinct structure where we have some areas here where there's a big difference between the three tissues. And you can see that the three tissues are not actually three tissues because it's different brain regions, but in these plots here, they are basically indistinguishable. We have liver, we have spleen, and we have brain regions, and we have sort of the same ordering. And that's actually true genome-wide that you're going to see this ordering everywhere. Then we have focal areas of low methylation. And the reason why this methylation dip doesn't go all the way down to 0% is really because we're smoothing over a big window. So there's always going to be some adjacent regions that's going to drive the average up. And then we have areas of very high methylation, substantially higher than the intermediate level. As soon as our methylation hits this intermediate level of methylation, we see this separation into the different tissues. These here were examples of— I would call them examples of partially methylated domains that was defined in this 2009 paper, where they defined partially methylated domains using a 70% threshold. So they basically looked at the methylation level inside these domains here. I will argue that this 70% threshold doesn't really make sense. What happens is that each tissue has its own baseline partially methylated domain level and gives rise to these sort of rather big differences between regions. But there's nothing special. Brain doesn't really have very much partially methylated domains at a 70% threshold, but in reality, they have these domains here. They're just high up because brain has a higher methylation level. We know that these dips with very low methylation level are likely regulatory regions like enhancers, promoters, regions of like potentially functional importance. Here's a different part of the genome. It's, I think, 12 megabases. Here we just see that these partially methylated domains are not everywhere. We can have like long stretches of areas where we don't see these partially methylated domains. And note here it's the same ordering we see again. When we zoom in and we look at a plot like we saw before, here we have 10 and 12 kilobases. We we see sort of this somewhat slightly, gives a little different impression. We have these sort of high levels of methylation, and then we have low levels of methylation, and then we have, we can see when we compare across tissues, like areas that differentially methylate regions where some tissue is different from the other tissues. This can happen either here like this, where we suddenly see a domain of low or a region of low methylation appearing, in one of the tissues, but it can also happen because the boundaries of these low methylation levels change. Here the different tissues disagree a little bit or are different in respect to how fast and when they go from sort of 0% methylation to high methylation. We can see that these— in these pictures here, we can see that this sort of classic pattern of low methylation on top of a CpG island. A final thing that occurs that is perhaps not appreciated so much is something I would call 'variable methylated regions,' and different people use this term 'variable methylated regions' in different ways, but I really— I think of it as— I intend for it to mean sort of a high degree of variation between replicates. That variation might be driven genetically or maybe, you know, some environmental exposure or some confounder or whatever, But basically, we take a lot of samples, we look at areas, and we see here very clearly that there's areas of the genome where there's just high variation. In this case here, we see kind of a tissue-specific, or in this case, we have 3 different brain regions, a brain region-specific level of variability. And so it's true over here to the left. Let's talk a little bit about what kind of differences or signal we expect to find when we do a DNA methylation study. Different studies have very different expected differences. The biggest type of differences are between cancer-normal or perhaps between tissues or between cell types. In those cases, we expect to see a lot of large differences, even more than you would expect to see with expression sometimes, in the sense that a large part of the methylome may be changing. In contrast, some comparisons with DNA methylation have minor expected differences. Here I'm typically thinking about sort of case-control comparisons within a cell type, although it depends a little bit on Things like, is the cell type you're profiling relevant to the exposure you may be studying or the intervention? And also, for example, does the exposure intervention result in some clear differences at the cellular level? Having some expectations on this really is important when you do quality control. Of your data, and it helps you guide several analysis choices. But the gamut, or the sort of the, the range of expected differences, is quite vast within the DNA methylation area. There's also different effect sizes. First, let me put in sort of an editorial comment and say that methylation in a single cell is really binary, or really has three different values, because we have 2 copies of each chromosome, but we either methylate or unmethylate. So every time we see methylation values that are different from 0 or 100%, we are really looking at something that has to be cell-to-cell variability in methylation. Furthermore, we almost never see changes in methylation going from 0 to 100 or from 100 to 0.. It's often perhaps, you know, a methylation change of 50% or so would be considered quite large, I would say. But it's really rare that it goes from like fully methylated to fully unmethylated or vice versa. Broadly speaking, I would say tissue comparisons results in many differentially methylated regions with large effect sizes, where many of them have large effect sizes, which I think of as like say 10 to 50 percentage points on the methylation scale, and we would expect these differentially methylated regions to overlap putative, at least, regulatory regions. If we also differentially express data, differentially expressed genes, and we have a differentially methylated region that's close to a differentially expressed gene, we often observe quite good correspondence between the magnitude of the change with methylation and the magnitude of the expressing change. When that is said, there are always many differentially expressed genes that don't have differentially methylated regions nearby and vice versa. So it's absolutely not, you know, a strong relationship, but there is a relationship. Cancer is similar in the sense that we find many DMRs often with large effect sizes, but perhaps I would, you know, somewhat generally say that we have poor expression correspondence in cancer. Then finally, sort of case-control studies within a cell type often report highly significant differences, at least if they have like thousands of samples, but sometimes these have tiny effect sizes down to 1% or 2% on the methylation scale. It's really an open question how functional they are, but irrespective of that, sometimes these small changes replicate. Let's talk a little bit about quality control and preprocessing. First, for sequencing. For sequencing, I found these unbiased plots that basically shows the fraction of methylation evidence along the read as being particularly revealing. They are produced by the Bismarck aligner, which is a popular aligner for bisulfite data, and I have many times, or several times at least, found unexpected issues with data that has been revealed with these M-bias plots. We really want the plots to be flat. The level of the plot can be a little sample-dependent, but they really ought to be flat assuming that everything is random, or the steering and the— basically that everything works well. In this case here, we are showing some data, some captured data, and for technical reasons, the start of the read basically is biased. There's going to be a mixture of sort of artificially methylated CBGs and sort of true signal. And this plot here suggests that we should trim off, say, the first 20 bases or so of the reads. PCA plots, which are widely used in genomics, is something that I tend to find quite useful for for DNA methylation data, as I do for most other expression data, for most other omics data. It can be a great way to discover sample mix-ups or issues with your data. It's of course helped a lot by having a good annotation you can use to color the PCA plots with, and it also helps a lot to have a decent amount of samples. In this dataset here, we have some PCA plots. We have a fair amount of whole genome bisulfite data, We average them into 1-kb bins, we discard bins that have very few CpGs in them, and then we do sort of a PCA plot here. And this plot actually revealed a couple of sample mix-ups where we are quite confident that we had mislabeled the brain samples that we were assaying. This is not just based on the PCA plot, we have sort of additional evidence for that. For methylation arrays, I tend to check the marginal beta distributions. And if you're familiar with methylation arrays, you know that they come in— or the arrays themselves have room for multiple samples, and they're typically processed with 96 samples at a time in a plate setup. You should really, when you, when you design an experiment, you should really try to randomize samples to the plates. But in practice, this is not always done and can therefore, if the samples are not randomized, plate effects can be confounded with comparisons of interest, and that can be a real issue. I also recommend using genotype probe to track samples on the methylation arrays. This is only really relevant if you have multiple samples, multiple samples from the same person. And finally, a cute little thing is that it's very easy to predict the sex of the individuals using DNA methylation data, and we have sometimes used that to track sample mix-ups by comparing predicted sex to reported sex. Now for preprocessing. Early on, I mentioned that methylation arrays or methylation data—bisulfite data in general—has kind of an internal control because you're missing a methylated and unmethylated channel. There's a lot of good normalization methods for Illumina methylation arrays, and— but broadly speaking, this existence of the internal control makes it less important to do preprocessing, and the choice of preprocessing has a smaller impact compared to expression microarrays. That doesn't mean it doesn't matter. It just means that the impact is smaller, and doing nothing is often surprisingly competitive in comparisons. If I have to cut to the chase and recommend something of all of these different methods, I tend to like to start with NOOB, partly because it has excellent performance, but also partly because it's a single sample method, which means that the normalization of one sample is independent of the other samples you have assayed. Almost all other normalization methods do some things where your results are going to change ever so slightly depending on which samples you normalized it together with. So NOOB has great performance and it's really convenient, and that's where I like to start. And I'm not an author of NOOB. Um, other alternatives that I always have in the back of my mind are SISME, sort of a new pipeline, uh, for that is particularly relevant. I would really look into that a lot if I was assessing cancer data. I have a method on functional normalization that is also relevant for cancer sometimes. With functional normalization, you can do a little bit better than NOOB, but not always. And then we have— we and other people have come up with sort of an adaptation of quantile normalization to methylation arrays. I'm mentioning this because I use this method if I really want to potentially even over-normalize the data. I may be in a situation where I really want to make sure that the signal I found is real, and I would rather normalize too much than too little. Let's talk a little bit about cell type composition. It's well understood that changes in cell type composition between samples, especially systematic changes in cell type composition between samples, may result in the appearance of differential methylation or false positives. This was illustrated quite well in an important paper by Hausmann from 2012. Because of this, a large number of algorithms have been proposed to estimate cell type composition from DNA methylation data. We can roughly divide them into reference and non-reference based. A reference-based algorithm means that we are using a reference panel of samples or sorted purified cell types that have been assayed, so we know what the sort of reference level of what the methylation landscape should look at in the different cell types. And non-reference-based methods basically don't use a reference method. So the reference-based algorithms really require the existence of these reference data. Hausmann has a method which he has expanded on that there's some later extensions on it, and there's a method called abidis from Zheng that I also recommend taking a look at. But this is a relatively crowded field, although there are sort of subtleties. A lot of the reference data is on methylation arrays. You may be interested in, like, using those in sequencing on sequencing data, And, yeah, another question you have is that after you've estimated, you may have, after you estimated a cell type composition, you may, and you see some differences, you may be interested in like asking, well, which of the cell types actually had the difference? And there's an extension to Ibitis that tries to answer that question. Most of the reference data for human is in blood and brain samples. And in brain, it's really neurons versus non-neurons. Cell type composition can also result in false negatives. That's in some way intuitively obvious, but we— I was recently involved in a study that really drove this home. And I would say eliminating false negatives might actually require sorting to overcome.. As an example, we studied methylation differences between brain regions, and we did it both in bulk data, unsorted data, and then we sorted on the NeuN marker, which basically means we are sorting into two populations: neurons and non-neurons. And doing this increased our number of differentially methylated regions by more than 1,000-fold. We basically went from finding nothing, like 100 DMRs or something like that genome-wide, to finding 100,000. This is, in this particular case, caused by the large technical variation in the spatial density of neurons in the brain. Here we have taken different samples from the same slice of the hippocampus. These two arrows indicate basically cutting out a very small piece of the slice and assaying And by basically moving the scalpel over a couple of millimeters and slicing a new part of the brain, the proportion of neurons in the same section changes from roughly 10% to 40%. This basically introduces a massive amount of technical variation in the amount of neurons we have in each sample and basically obscures the signal. And sorting the cells on this new end marker really dramatically increased our power, like dramatically. This might be, you know, this might be relevant to most tissue samples, but perhaps this spatial confounding is particularly relevant in something like that. Let's talk about genetics and DNA methylation. It's well recognized that a genetic variant can control or influence the DNA methylation in a nearby or in a distal site. This is called a methylation QTL, quantitative trait loci. And is the end goal of quite a few methylation studies. There are two important additional scenarios describing the interplay between DNA methylation and genetics which are relevant. Situation 1: a genetic variant may affect the ability to quantify DNA methylation and thereby wrongly appears at methylation QTL. In this example below here, we have a variant, we have a CpG we're interested in measuring, We have a genetic variant amongst people and we have a probe that is designed against the reference genome. So individuals with this variant present, the probe will not match the reference genome and this will impede the ability of the probe to measure DNA methylation correctly. This is a case not of a— this will appear as a methylation QTL, but it's not really. It's basically a problem with quantification. Situation number 2 is the fact that CpG to TPG variants are the most common SNPs in the genome due to the mutational pressure of DNA methylation in the germline. These variants are invisible to bisulfite conversion and they register as 0% methylation because a TPG is evidence for a lack of methylation. How you want to think of these variants is a bit philosophical, and that's because we really have like these 3 situations here. We can have a CpG, whereas a C is ready to take methylation and may be more or less methylated. Then we can have a CpG that is unmethylated, and then finally we can have a TpG. And with bisulfite conversion, the last two situations are basically impossible to differentiate, or at least naively impossible to differentiate. And whether or not you think of these as different depends a little bit on your point of view. I mean, clearly there is no CpG methylation when you have a T instead of a C, but you can argue that you also don't have the ability to be methylated. This is important to be aware of because these SNPs are quite prevalent in the human genome or any mammalian genome. And I'm going to note that you can think of it as a change in methylation, right? In terms of evolution and selection, the easiest way to remove methylation from a part of the genome is to basically you know, a tolerant variant where the CpG gets mutated away. So the question you have in analysis is, should these situations be excluded or should they be studied for their interest? A common strategy, especially for analysis of methylation arrays, is exclusion. And some people or some work exclude as much as more than 100,000 CpGs from consideration because they overlap SNPs. This is particularly relevant when we compare mice strains between each other. Two different mice strains are typically, typically separated by millions of years in evolution, and we therefore have sometimes millions of CPGs that are gained and lost between the strains. We started this extensively in a paper from 2019 where we proposed to use smoothing as imputation. The region we have here is an example of the kind of pattern you may see when you don't take this into account. Down here, the top tick marks are common CpGs, and the bottom tick marks are CpGs that are only present in one of the two reference genomes. In this, in this picture, we have taken all our samples from the two different mice strains and aligned them to BLAC6. And what happens is that the Cas genome has lost CpGs, and that appears as lower methylation. If we take the exact same data and align the exact same region to the Cas genome, we now see that the BLAST6 genome is less methylated than the Cas genome, and basically completely reversing the effect. And we have a proposal for dealing with this. We can think of smoothing as a form for imputation, and we can think of smoothing as computing kind of a local percent methylation in a regional fashion. And doing this, aligning to strain-specific genomes and thinking about moving the data between the two different coordinate systems, results in basically no difference between the two genomes in this particular region. This is an interesting situation. The situation is much stronger for mice than it is for humans, simply because there's much more difference between— genetic differences between two mouse strains and between two humans we've sampled. But it's worth thinking about. Finally, let's talk a little bit about differential methylation analysis. Broadly speaking, you can think of analysis of methylation data as being focused on finding differentially methylated positions or CpGs or single CpGs. You do an analysis that's unaffected by the nearby CpGs, and then there's another set of methods where you arrive at differentially methylated regions where you want to find many CpGs that are jointly differentially methylated between your conditions. You, when you do this analysis, you really need to consider confounders, especially all type of genomics data is affected by batch effects, and we have a lot of knowledge about this, and it goes way beyond this lecture to talk about this. And then it's also common to think about cell type composition. There's of course other confounders, possible confounders such as age or other aspects. So looking for differentially methylated positions is a type of analysis that's quite well understood. That doesn't mean that it's easy to do, but it's well understood. We basically treat each CBG individually and have a number of CBGs, and that basically means this type of analysis has many similarities to expression analysis, or for that matter, to genome-wide or GWAS. Key points in your analysis pipeline is thinking about filtering and removing CBGs from considerations. We typically apply various forms of linear models to include covariates. If you have arrays, oh, I guess for sequencing as well, you could consider transforming from the beta scale or the percent scale into M values, or basically take the logit of the percent methylation. And if you're doing sequencing, you really need a model that considers coverage, for example, through a negative binomial model or by modeling the mean-variance relationship. And there's a huge variety of literature on these aspects, and some. There's also a lot of useful work in expression analysis and other omics types of analysis. It is pretty well understood how to account for batch effects in this situation using, for example, surrogate value analysis or remove unwanted variation. And if you're using arrays, there's an easy-to-use adaptation of IUV in the missmethl package. I'm going to highlight— this is the missmethl package and The Ms. Methyl authors also has an extended analysis workflow for methylation arrays in F1000 Research that I really would recommend as a great starting point for learning how to do analysis of methylation arrays. Differentially methylated regions can roughly be obtained by kind of first doing an analysis of individual CpGs and then aggregating them into regions, or through some kind of statistical analysis where you directly discover differentially methylated regions. For sequencing, there's also a wide variety of methods far beyond what's listed here. I've developed a method called B-SMOOTH. My advisor later on extended that into something called DMR-Seq, and there's DMR-CATE that's also a great method. Most of the work on statistical analysis of whole genome bisulfite sequencing has been done and been assessed using very small sample sizes. There's a lot of focus on methods that work well with 1 or 2 or 3 samples, and I think that's a little bit problematic. We know that DNA methylation varies between people, but it's really mostly a cost issue. These days sequencing, whole genome bisulfite sequencing, is starting to get cheap enough that we are beginning to see studies with tens or hundreds of samples, and it's going to be interesting to see what kind of methods emerges here. It's also possible to go beyond sort of the CpG-level analysis and analyze things at sort of the domain level. We call this block finding. It's more unusual to do. It's a little bit similar, if you think back on the talk I had on structure of DNA methylation, to find basically large domains of difference of different methylation between the samples, which usually corresponds to kind of PMTs. Finally, after differential analysis, many people are interested in doing some type of gene set enrichment analysis or gene ontology analysis to kind of potentially provide some insight into the result. It's well understood that it is critical that when you do such analysis, you have to account for the unusual spatial distribution of CpGs. This is usually done by modifying the background distribution utilized in these gene set enrichment methods, and this is typically not done for tools— if you have access to tools that were developed for analysis of RNA expression data. There's a— one of the first papers I read on this was from 2013, and there was a new publication that does this at the differentially methylated region level instead of at the individual CpGs levels in Genome Biology that I would recommend having a look at. This is really critical to think about. So this concludes the lecture on computational analysis of DNA methylation data. I've only like really scratched the surface. There's many specific analysis choices and methods that we don't have time to go into, but I hope this will provide a good starting point for reading the literature.