mesa analyses DNA methylation measured by enrichment sequencing, such as
MeDIP-seq (methylated DNA immunoprecipitation) and MBD-seq (methyl-CpG
binding domain capture), including low-input samples such as cell-free DNA.
These assays pull down methylated fragments instead of converting
unmethylated cytosines, so the signal is a read count per genomic window
whose relation to methylation depends on the local CpG density. mesa builds
on the qseaSet class from qsea, which models that relation, and adds
functions to build, quality-control, explore and compare qseaSet objects
in a tidyverse style.
A typical analysis follows these steps:
Build a qseaSet from BAM files with makeQset(), which tiles the
genome into windows, counts fragments and estimates copy number (see
addHMMcopyCNV()). addNormalisation() fits the CpG-density-dependent
enrichment used to estimate beta values.
Check quality with getSampleQCSummary() and
plotCorrelationMatrix(), and subset samples or windows with
filter(), filterWindows() or subsetQset().
Explore sample structure with principal components (getPCA()) or
UMAP (getUMAP()).
Compare groups of samples with calculateDMRs(), which fits a
negative binomial model per window and returns differentially methylated
regions (DMRs), and summarise them with summariseDMRsByContrast().
Extract counts, normalised reads per million (nrpm) or beta values
per window with getDataTable(), for example over the DMRs after
filterByOverlaps().
Annotate windows or DMRs with nearby genes and genomic features with
annotateWindows() and plotGenomicFeatureDistribution().
Two small qseaSet objects are shipped for examples:
exampleTumourNormal (human, paired tumour/normal samples) and
exampleMouse (mouse cell line).
vignette("introduction", package = "mesa"): overview and quick start.
vignette("generation", package = "mesa"): building qseaSet objects.
vignette("data-and-qc", package = "mesa"): data tables and quality
control.
vignette("pca", package = "mesa"): PCA and UMAP.
vignette("differentially-methylated-regions", package = "mesa"): DMRs
and their annotation.
Useful links:
# A qseaSet of 10 samples (5 tumour/normal pairs) over windows on chr7
exampleTumourNormal
#> qseaSet
#> =======================================
#> 10 Samples:
#> Colon1_N, Colon1_T, Colon2_N, Colon2_T, Lung1_N, Lung1_T, Lung2_N, Lung2_T, Lung3_N, Lung3_T
#> 819 Regions in 22 chromosomes of BSgenome.Hsapiens.NCBI.GRCh38
getSampleTable(exampleTumourNormal)
#> sample_name group patient type tumour tissue age gender
#> Colon1_N Colon1_N Colon1_N Colon1 NormalColon Normal Colon 80 F
#> Colon1_T Colon1_T Colon1_T Colon1 CRC Tumour Colon 80 F
#> Colon2_N Colon2_N Colon2_N Colon2 NormalColon Normal Colon 67 F
#> Colon2_T Colon2_T Colon2_T Colon2 CRC Tumour Colon 67 F
#> Lung1_N Lung1_N Lung1_N Lung1 NormalLung Normal Lung 64 F
#> Lung1_T Lung1_T Lung1_T Lung1 LUAD Tumour Lung 64 F
#> Lung2_N Lung2_N Lung2_N Lung2 NormalLung Normal Lung 71 F
#> Lung2_T Lung2_T Lung2_T Lung2 LUSC Tumour Lung 71 F
#> Lung3_N Lung3_N Lung3_N Lung3 NormalLung Normal Lung 72 M
#> Lung3_T Lung3_T Lung3_T Lung3 LUAD Tumour Lung 72 M
#> stage
#> Colon1_N <NA>
#> Colon1_T III
#> Colon2_N <NA>
#> Colon2_T IV
#> Lung1_N <NA>
#> Lung1_T IIA
#> Lung2_N <NA>
#> Lung2_T IA
#> Lung3_N <NA>
#> Lung3_T IB
# Beta values (methylation estimates) per window and sample
exampleTumourNormal %>%
getDataTable(normMethod = "beta") %>%
head()
#> Generating table of beta values for 819 regions across 10 samples.
#> # A tibble: 6 × 14
#> seqnames start end CpG_density Colon1_N Colon1_T Colon2_N Colon2_T Lung1_N
#> <fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 7 2.50e7 2.50e7 5.51 NA 0.347 NA 0.369 0.358
#> 2 7 2.50e7 2.50e7 13.4 0.0565 0.0171 0.0275 0.0165 0.0171
#> 3 7 2.50e7 2.50e7 5.57 NA 0.204 NA 0.363 0.353
#> 4 7 2.50e7 2.50e7 5.10 NA 0.382 NA 0.406 0.393
#> 5 7 2.50e7 2.50e7 8.38 0.122 0.0695 0.918 0.818 0.901
#> 6 7 2.50e7 2.50e7 17.9 0.0240 0.0660 0.265 0.250 0.161
#> # ℹ 5 more variables: Lung1_T <dbl>, Lung2_N <dbl>, Lung2_T <dbl>,
#> # Lung3_N <dbl>, Lung3_T <dbl>
# Differentially methylated regions between tumour and normal samples
exampleTumourNormal %>%
calculateDMRs(variable = "tumour", contrasts = "first") %>%
head()
#> Calculating the first possible contrast on the tumour column.
#> Generating table of nrpm values for 819 regions across 10 samples.
#> Fitting initial GLM on 309 windows, without using parallelisation.
#> Performing contrast Normal_vs_Tumour
#> Contrasts Complete
#> selecting regions from contrast Normal_vs_Tumour
#> selected 23/309 windows
#> adding test results from Normal_vs_Tumour
#> # A tibble: 6 × 11
#> seqnames start end CpG_density Normal_vs_Tumour_log2FC
#> <fct> <int> <int> <dbl> <dbl>
#> 1 7 27102901 27103200 10.8 -1.94
#> 2 7 27103501 27103800 13.1 -1.88
#> 3 7 27103801 27104100 7.83 -1.33
#> 4 7 27108301 27108600 15.7 -1.18
#> 5 7 27108601 27108900 16.1 -0.933
#> 6 7 27144601 27144900 11.1 -1.47
#> # ℹ 6 more variables: Normal_vs_Tumour_adjPval <dbl>,
#> # Normal_vs_Tumour_deltaBeta <dbl>, Normal_beta_means <dbl>,
#> # Tumour_beta_means <dbl>, Normal_nrpm_means <dbl>, Tumour_nrpm_means <dbl>