--- title: "Data Tables and Quality Control" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Data Tables and Quality Control} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r chunk, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, out.width = "100%", dev = "ragg_png" ) ``` ```{r setup} library(mesa) library(dplyr) library(ggplot2) library(tidyr) library(stringr) ``` # Generating tables of values from a qseaSet `qsea` has several different types of values that can be returned per window. You can get a table of raw counts per window with `getCountTable`: ```{r generating-tables-of-values-fr} exampleTumourNormal %>% getCountTable() ``` but this is not particularly useful, as different samples have generally been sequenced to different read depths. Note that the annotation of the `qseaSet` regions is included in this output, notably the `CpG_density` that is described in `vignette("generation", package = "mesa")`. The next option is normalised reads per million (NRPM). This normalises the number of reads by the depth of sequencing, as well as by copy number variation. ```{r chunk-2} exampleTumourNormal %>% getNRPMTable() ``` Finally, we can request beta values. Beta values are measures of methylation, scaled to lie between 0 and 1, where 0 means that the window is unmethylated and 1 is methylated. These are scaled with respect to a number of further metrics, including the region's `CpG_density`, background noise and the sample enrichment profile, and are therefore the most useful value to consider when comparing methylation between samples. ```{r chunk-3} exampleTumourNormal %>% getBetaTable() ``` Note that for some windows an NA value will be returned. This occurs when fewer than three reads would be required for the window to be fully methylated, at the read depth and enrichment profile calculated per sample. These occur particularly at windows with a low `CpG_density`. All three of these functions take an option `useGroupMeans` that determines whether to average over the `group` column of the qseaSet `sampleTable`, which can be useful when you have biological or technical replicates for example. In the example qseaSet, `group` is the same as `sample_name`, so we need to change that to see the effect: ```{r chunk-4} exampleTumourNormal %>% getBetaTable(useGroupMeans = FALSE) exampleTumourNormal %>% mutate(group = str_remove(sample_name,"\\d")) %>% getBetaTable(useGroupMeans = TRUE) ``` Now the table has one column per group (four here) rather than one per sample. Windows where no sample in a group has a beta value are returned as `NaN`. We explicitly note that qsea incorporates the uncertainty of the measurements in its beta value estimations. This means that the values returned as `beta` values are actually the median of this estimate, with an associated confidence interval (which can be queried, see `qsea::makeTable`). In practice, this means that samples with high coverage can return values very close to 0 and 1, but the same sample with less coverage might return 0.1 or 0.9. This can have an effect on clustering and PCAs for example. # Quality Control One extremely important aspect of any sequencing experiment is quality control (QC), and ME-seq is no exception. Alongside the typical QC steps on the bam files themselves (although a raised CG content should obviously be expected following enrichment), there are several methylation-specific steps to consider. ## Relative Enrichment The first metric used for QC is a relative enrichment metric, known as `relH`. This was first introduced in the `MEDIPS` package, via the `MEDIPS.CpGenrich()` function, and is a metric to express the proportion of CG dinucleotides, compared to the reference genome. The relative enrichment `relH` is calculated in the following way: * For each mapped DNA fragment, retrieve the end-to-end DNA spanned on the reference sequence. This is used instead of the reads themselves, as a fragment will often be longer than the read(s) sequenced at its ends, and may therefore contain CGs not covered by the reads. * In these sequences, count how many CGs are present in total, and divide by the total number of base pairs covered. * Divide by the same ratio, calculated for the entire genome. This is therefore a measure of how many more CGs are present in the fragments covered within the bam file, compared to the reference genome. Within `mesa`, this metric is calculated through `makeQset`, as part of the information added to the `libraries` table on the `qseaSet`. As a rough guideline, for the T7-MBD-seq data studied by the package authors, `relH` values of 4-5 were typical, with a minimum of 2.5-3 regarded as acceptable. Typical values may depend on enrichment method and sample type. We can use the `addLibraryInformation` function to add the information in the `libraries` table directly onto the `sampleTable`: ```{r chunk-5} exampleTumourNormal %>% addLibraryInformation() %>% getSampleTable() %>% dplyr::select(sample_name, relH) ``` When a non-enriched ("Input") sample is used to calculate copy number variation, an equivalent calculation is also performed on these non-enriched fragments. This is also stored in the `libraries` table, in columns starting with `input_`: ```{r chunk-6} exampleTumourNormal %>% addLibraryInformation() %>% getSampleTable() %>% dplyr::select(sample_name, input_relH) ``` Values of `input_relH` between 1.3 and 1.6 were observed in non-enriched data studied by the package authors. These values above 1 indicate that more CGs are seen in the sequenced fragments than would be expected from the genomic sequence. This is possibly due to the evolutionary suppression of `CG` dinucleotides outside of specific contexts, meaning they are less common in those areas of the genome that are repetitive and hence less easy to map. This is likely to depend on the particular sequencing methods and the sample type under consideration. We expect to see an increase in the enrichment metric between the enriched and non-enriched data: ```{r chunk-7} exampleTumourNormal %>% addLibraryInformation() %>% getSampleTable() %>% dplyr::select(sample_name, relH_MeCap = relH, relH_Input = input_relH) %>% pivot_longer(matches("relH"), names_prefix = "relH_", names_to = "enrichment", values_to = "relH") %>% ggplot(aes(x = sample_name, y = relH, col = enrichment)) + geom_point() + theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) ``` There is a second relative enrichment metric defined by the `MEDIPS` package, known as `GoGe`, that is also calculated globally in a similar way but considers the ratio between the number of `CG`s present in the fragments compared to the number of `C`s and `G`s independently. This is calculated and stored in the `libraries` information alongside the `relH` metric. ## Fragments without pattern The `libraries` table also records how many fragments do not contain a CG at all. This can be a measure of background noise. ## Hyper-stably methylated regions Another metric used for quality control is based on windows that are known to be commonly methylated across a range of tissues and diseases. For instance, for the human genome, [Edgar et al. (2014)](https://doi.org/10.1186/1756-8935-7-28) identified probes from the Illumina 450k methylation array that were consistently methylated in over 1700 samples. We take windows covered by these probes, calculate beta values for them, and then determine the proportion of these windows that have a beta value above a threshold of 0.8; we call this metric `hyperStableEdgar`. Samples with lower enrichment scores (`relH`) require a greater coverage depth for beta values to reach a given threshold. The metric therefore incorporates both the coverage depth and how well the enrichment has worked. As noted above, windows with low `CpG_density` are more likely to have an NA beta value. To avoid considering windows that typically have a missing beta value across most samples, the calculation of `hyperStableEdgar` is restricted to windows that have a `CpG_density` of at least 5. Samples with insufficient coverage and/or enrichment will still have missing beta values for some windows with `CpG_density` above this minimum value. These missing beta values are treated as zero when calculating `hyperStableEdgar`. For hg38, these hyper-stable methylated probes are provided in a data object, `hg38UltraStableProbes`, and the function `addHyperStableFraction` calculates the `hyperStableEdgar` metric and adds it as a column to the sampleTable. Note that the function will not work with the provided example data, as there is no overlap of the probes with the windows in the example dataset (i.e. `filterByOverlaps(exampleTumourNormal, hg38UltraStableProbes)` has no windows). For the T7-MBD-seq data studied by the package authors, a `hyperStableEdgar` value of at least 0.4 was considered acceptable. Again, suitable values may depend on enrichment method and sample type. ## Summary In the package authors' experience, the three metrics `valid_fragments`, `relH` and `hyperStableEdgar` are the most useful for quality control. For 300 bp windows on the human genome with T7-MBD-seq data, the authors considered a `relH` of at least 2.5–3, a `hyperStableEdgar` value of at least 0.4 (40%), and around 3,000,000 `valid_fragments` acceptable. Samples that met the `relH` and `hyperStableEdgar` thresholds usually also met the `valid_fragments` threshold. These values are provided as guidelines rather than strict thresholds; suitable values will depend on the enrichment protocol (MBD-seq or MeDIP-seq), the sample type and the sequencing depth, so we recommend choosing thresholds from the distribution of these metrics in your own data. ## Session info ```{r session-info, echo=FALSE} sessionInfo() ```