Initial QC and summary

Without any further processing, we can create an interactive html report containing basic summary statistics of the methrix object with methrix_report function.

The report can be accessed here: Initial reports

methrix::methrix_report(meth = meth, output_dir = getwd(), prefix = "initial")

In the report, we can check genome-wide and chromosome based statistics on coverage and methylation, in order to identify potential quality issues, for example:

samples with too high or too low genome-wide coverage or on a specific chromosome. In this case, there is also a possibility, that the sample is affected by copy number alterations. Both methylation levels and coverage therefore shows higher variation among cancer samples.
samples with altered beta level distribution.
A “bump” on the beta level density plot at 0.5 indicates the presence of single nucleotide polymorphisms (SNPs) -> see later.

Filtering

Coverage based filtering

To ensure the high quality of our dataset, the sites with very low (untrustworthy methylation level) or high coverage (technical problem) should not be used. We can mask these CpG sites. Please note, that the DSS we were using for Differential Methylation Region (DMR) calling, doesn’t need the removal of the lowly covered sites, because it takes it into account by analysis. However, using a not too restrictive threshold won’t interfere with DMR calling.

meth <- methrix::mask_methrix(meth, low_count = 2, high_quantile = 0.99)

Later we can remove those sites that are not covered by any of the samples.

meth <- methrix::remove_uncovered(meth)

We might want to remove the sites that are sparsely covered. The function coverage_filter has multiple options on how to perform this. We can set a threshold on the coverage and use either min_samples (minimum number of samples) or prop_samples (minimum proportion of the samples) to define the minimum size of samples that should excess the threshold for each CpG side. bsseq requires

#meth <- methrix::coverage_filter(m=meth, cov_thr = 5, min_samples = 2)
#or
meth <- methrix::coverage_filter(m=meth, cov_thr = 5, min_samples = 2, group="Day")

SNP filtering

SNPs, mostly the C > T mutations can disrupt methylation calling. Therefore it is essential to remove CpG sites that overlap with common variants, if we have variation in our study population. For example working with human, unpaired data, e.g. treated vs. untreated groups. During filtering, we can select the population of interest and the minimum minor allele frequency (MAF) of the SNPs. SNP filtering is currently implemented for hg19 and hg38 assemblies.

if(!requireNamespace("GenomicScores")) {
BiocManager::install("GenomicScores")
}
library(GenomicScores)

meth <- methrix::remove_snps(m = meth, keep = TRUE)

There is a possibility to keep the sites possibly overlapping with SNPs. We can use this data for detecting possible sample mismatches. For this, we need to plot the most variable sites (probably the ones that are variable in our study population) and plot them on a heatmap.

if (!requireNamespace("pheatmap", quietly = TRUE))
    install.packages("pheatmap")
snps <- meth[[2]]
meth <- meth[[1]]
pheatmap::pheatmap(get_matrix(order_by_sd(snps)[1:min(5000, length(snps)),]))

SNP removal in mouse data: Usually, there is no need to remove SNPs from mouse data, since for many experiments, the same strain of mice is used. However, there are cases, when different crosses of strains are used, for which, the level of heterozygosity might be different. SNP data for mice can be downloaded from here: https://www.sanger.ac.uk/science/data/mouse-genomes-project. The SNPs can be filtered out using the region filter function.

Visualization and QC after filtering

Methrix report

After filtering, it worth running a report again, to see if any of the samples were so lowly covered that it warrants an action (e.g. removal of the sample)

Important: Don’t forget to set the recal_stats to TRUE, since the object changed since reading in. The output directory has to be different from last time, to avoid using the intermediate files calculated during the last run.

methrix::methrix_report(meth = meth, output_dir = getwd(), recal_stats = TRUE, prefix="processed")

The report can be found here: processed_methrix_reports.html

There are additional possibilities to visualize the study. We can look at the density of coverage both sample- and group-wise:

methrix::plot_coverage(m = meth, type = "dens")

methrix::plot_coverage(m = meth, type = "dens", pheno = "Day", perGroup = TRUE)

We can visualize te beta value distribution as a violin plot or density plot. These plots (as well as the coverage plot) use the most variable 25000 CpG sites to ensure fast processing. plot_density and plot_violin has the option to use data only from a restricted region.

methrix::plot_density(m=meth)

methrix::plot_density(m=meth, ranges = GRanges("chr22:48007313-49007313"))

Dimensionality reduction

Methrix offers principal component analysis (PCA) to conduct and visualize dimensionality reduction. As a first step, the PCA model has to be calculated. As for plotting, a given number of sites are selected, either random or based on variance (var option). This ensures that the calculations remain feasible.

mpca <- methrix::methrix_pca(m= meth, top_var = 10000, n_pc = 20)

At visualization, we can provide the methrix object to allow color or shape annotation of groups or samples.

methrix::plot_pca(mpca, m=meth, col_anno = "Day", shape_anno = "Replicate")

methrix offers the possibility of region based filtering. With this option, selected regions, e.g. promoters can be visualized. We will use the annotatr package to assess basic annotation categories.

if(!requireNamespace("annotatr")) {
BiocManager::install("annotatr", update = F)
}

library(annotatr)
promoters = build_annotations(genome = 'hg19', annotations = "hg19_genes_promoters")
promoters

mpca <- methrix::methrix_pca(m= meth, top_var = 10000, ranges = promoters, n_pc = 3, do_plot = F)
methrix::plot_pca(mpca, m=meth, col_anno = "Day", shape_anno = "Replicate")