Borealis outlier methylation detection
30 December 2022
Package
borealis 1.2.2
Contents
- 1 Introduction
- 2 Installation
- 3 Running Borealis
- 4 Basic post-processing and analysis of Borealis results
- 4.1 Read in entire cohort’s results
- 4.2 Generating summary metrics across all samples
- 4.2.1 How many CpG sites worth of data do we have across all samples combined?
- 4.2.2 How many unique samples and unique CpG sites are we analyzing data from?
- 4.2.3 Distribution of the mean methylation values and variability per CpG Site
- 4.2.4 Summary of read depth distributions
- 4.2.5 Summary of uncorrected p-values
- 4.2.6 Summary of corrected p-values
- 4.2.7 Summary of methylation fraction across sites
- 4.2.8 Summary of effect sizes
- 4.2.9 Outliers and most significant CpG site
- 4.3 Annotating outputs with epigenetic features
- 4.4 Summarizing single-site data across epigenetic features
- Session info
1 Introduction
Borealis is a package for the detection of outlier methylation at single CpG-site resolution, where a cohort of 3 to 100+ samples can be processed and each sample is analyzed versus the rest of the cohort to identify outlier hypermethylated or hypomethylated CpG sites. This form of one vs many analysis differs from traditional case vs control group analyses and has been successful in domains such as rare genetic disease driver identification. Furthermore, the ability of Borealis to identify single CpG-site differences offers a higher resolution view of methylation, since increasing numbers of studies demonstrate single-site methylation aberrations in disease.
This vignette provides an introduction to some basic and advanced operations with Borealis using a region of chromosome 14 in a cohort of 20 individuals being investigated for causes of rare genetic disease. In real-world use, Borealis expects outputs from the aligner Bismark which is available on Github.
After completing this vignette users should be able to conduct in-depth methylation analysis and discover biological signals in their data. After learning basic functionality, creating summary metrics and looking for outlier samples, users will optionally learn how to annotate CpGs with epigenetic feature context using annotatr. Finally users will learn how to summarize CpG data across epigenetic features.
1.1 Citation
If you use borealis in your research, please cite Oliver et al. (2022).
2 Installation
The release version of borealis is available via Bioconductor and can be installed as follows:
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("borealis")The development version of borealis can be obtained via the Github repository.
It is easiest to install development versions with the CRANpkg(“devtools”) package as follows:
devtools::install_github('GarrettJenkinson/borealis')Changelogs for development releases will be detailed on GitHub releases.
3 Running Borealis
Now let’s load the test data included with Borealis. This represents a specific region from chromosome 14 (hg19) for 20 individuals with undiagnosed rare disease.
The entire Borealis pipeline can be run with the single runBorealis command.
library("borealis")
# Set data locations
outdir <- tempdir()
extdata <- system.file("extdata","bismark", package="borealis")
# Run borealis
results <- runBorealis(extdata,nThreads=2, chrs="chr14", suffix=".gz",
outprefix=file.path(outdir,"vignette_borealis_"),
modelOutPrefix=file.path(outdir,"vignette_CpG_model"))Now let’s quickly make sure we generated one output per sample:
# Read in the name of all files for analysis
borealis_in <- dir(outdir,pattern="*DMLs.tsv")
length(borealis_in)## [1] 20Looks good? Let’s continue!
4 Basic post-processing and analysis of Borealis results
4.1 Read in entire cohort’s results
Now you should have successfully loaded the provided methylation data, run Borealis and created a list of its output files.
Let’s now have a look at the data and generate some summary metrics.
First we’ll read in the data for each of the 20 patients and create GRanges for each:
# Read in list of Borealis output files and create a dataframe for each
for (file in borealis_in) {
name <- sub("_DMLs.tsv", "", file)
assign(name,GenomicRanges::makeGRangesFromDataFrame(
read.table(file.path(outdir,file), header=TRUE,
stringsAsFactors=FALSE), start.field="pos", end.field="pos.1",
seqnames.field="chr", keep.extra.columns=TRUE))
}
# Create a list of the newly created dataframes
list_object_names <- ls(pattern="borealis_patient")
listGRanges <- lapply(list_object_names, get)Let’s check we have 20 new GRange objects and confirm the naming convention:
length(list_object_names)## [1] 20list_object_names[1]## [1] "vignette_borealis_patient_70_chr14"Looks good - now let’s have a look at one row to get familiar with the fields of data:
listGRanges[[1]][1,]## GRanges object with 1 range and 7 metadata columns:
## seqnames ranges strand | x n mu theta
## <Rle> <IRanges> <Rle> | <integer> <integer> <numeric> <numeric>
## [1] chr14 24802152 * | 11 18 0.21535 0.0502731
## pVal isHypo effSize
## <numeric> <logical> <numeric>
## [1] 0.0156873 FALSE 0.395761
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsTo explain the data you see in brief:
- x = number of methylated CpGs at a given position
- n = the toal number of reads at that position
- mu = the mean fraction of reads methylated at this position across all samples
- theta = the disperson of the methylation values at this position
- pVal= probability of deviation from normal methylation at this position
- isHypo = if the site is hypomethylated at this position
- Yes = the site is defined as hypomethylated
- No = the site is defined as hypermethylated
- NA = the change is not significant enough to call hyper/hypo
- effSize = the effect size i.e. how large a deviation from normal methylation
Okay - now let’s start summarizing some of this information across samples.
4.2 Generating summary metrics across all samples
First we will add two new columns to each GRange. One column will be the sample ID for each patient, to enable us to distinguish between results later. The second will contain p-values adjusted for multiple comparisons at each CpG site per-patient.
We decouple the generation of adjusted p-values from the main package since users may wish to use uncorrected p-values or generate adjusted p-values in a specific fashion depending on their analysis.
# Add sample ID and a corrected p-value to each and output as new files (.padj)
for (i in seq_along(listGRanges)) {
sample=sub("_chr.*", "", list_object_names[i])
listGRanges[[i]]$sampleID <- sample
listGRanges[[i]]$pAdj <- p.adjust( listGRanges[[i]]$pVal, method="BH")
}Let’s have a look at the first entry to see those newly added columns:
listGRanges[[1]][1,]## GRanges object with 1 range and 9 metadata columns:
## seqnames ranges strand | x n mu theta
## <Rle> <IRanges> <Rle> | <integer> <integer> <numeric> <numeric>
## [1] chr14 24802152 * | 11 18 0.21535 0.0502731
## pVal isHypo effSize sampleID pAdj
## <numeric> <logical> <numeric> <character> <numeric>
## [1] 0.0156873 FALSE 0.395761 vignette_borealis_pa.. 1
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengthsLooking good. We now have adjusted p-values and sample identifiers.
Now we will create a single dataframe containing data for all samples and generate our summary metrics.
# Create a single GRanges obect with data for all samples and a dataframe for summaries
combined_files <- Reduce(c,listGRanges)
combined_files_df<-data.frame(combined_files)4.2.1 How many CpG sites worth of data do we have across all samples combined?
# How many rows of data in combined GRange object?
length(combined_files)## [1] 3756Lets create a table to summarize useful metrics about the CpG sites
# Create table of unique positions and mu/theta values
mu_theta <- unique(subset(combined_files_df,
select=-c(x,n,pVal, isHypo, pAdj, effSize, sampleID)))4.2.2 How many unique samples and unique CpG sites are we analyzing data from?
#Number of unique samples
length(unique(combined_files_df$sampleID))## [1] 20#Number of unique CpG sites
nrow(unique(mu_theta))## [1] 2474.2.3 Distribution of the mean methylation values and variability per CpG Site
#generate summaries for mu and theta
summary(mu_theta$mu)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.028 0.068 0.129 0.341 0.789 0.921summary(mu_theta$theta)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0141 0.0238 0.1417Now let’s look at the distribitions of other important values.
# Create table of unique positions and depth/p-val/padj for each position in
# each case
depth_pvals_eff <- unique(combined_files_df)4.2.4 Summary of read depth distributions
#Summarize read depths
summary(depth_pvals_eff$n)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.0 15.0 22.0 24.6 33.0 71.04.2.5 Summary of uncorrected p-values
#Summarize pvals
summary(depth_pvals_eff$pVal)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001 0.360 0.623 0.605 0.881 1.0004.2.6 Summary of corrected p-values
#Summarize corrected pvals
summary(depth_pvals_eff$pAdj)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.019 1.000 1.000 0.976 1.000 1.0004.2.7 Summary of methylation fraction across sites
#Summarize fraction of methylation
summary(depth_pvals_eff$x/depth_pvals_eff$n)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.045 0.266 0.571 1.0004.2.8 Summary of effect sizes
#Summarize effect size
summary(depth_pvals_eff$effSize)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.483 -0.068 -0.039 -0.015 0.022 0.5954.2.9 Outliers and most significant CpG site
Now lets see if any samples are extreme outliers in terms of number of CpG sites called significant (p <= 0.05):
# Detection of outliers based on number of significant sites
# Count significant CpG sites per patient
signif_only <- subset(combined_files_df, pVal <= 0.05)
signif_counts <- dplyr::count(as.data.frame(signif_only),sampleID)
# Calculate the percentiles of the number of significant sites
sig_quantiles <- quantile(signif_counts$n,
probs=c(0.025, 0.05, 0.95, 0.975, 0.999))
# Check if nay patients are above the 99.9th percentile
subset(signif_counts, n >= sig_quantiles["99.9%"])## sampleID n
## 3 vignette_borealis_patient_72 54We can see that one patient appears to be an extreme outlier.
Let’s also have a look at which site is most significant.
# What is the most significant CpG site between all samples?
subset(combined_files_df, pVal == min(combined_files$pVal))## seqnames start end width strand x n mu theta pVal isHypo
## 371 chr14 24780569 24780569 1 * 17 28 0.117 0.0567 0.000592 FALSE
## effSize sampleID pAdj
## 371 0.49 vignette_borealis_patient_72 0.0192Depending on your use-case, this might be enough information for you and you could stop here, or you could extract the top 100 outlier sites.
Another useful step is contextualizing CpG sites on the basis of which epigenetic features they reside in. We provide an example of how to do so in the following section. We will do this for the single sample containing the most significant CpG site for the sake of brevity in the vignette.
4.3 Annotating outputs with epigenetic features
In this section we will add further biological context to our results by annotating epigenetically relevant features. We will use the annotatr package but you can use your favorite tool instead. The release version of annotatr is available via Bioconductor and can be installed as follows:
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("annotatr")Now let’s load the package and read in the data we created in the previous steps for the patient with the most significant CpG site.
We will annotate with annotations from annotatr.
Let’s will start by defining the annotations we want to utilize.
#Assign approproate genome version based on alignments
genome.version <- "hg19"
# Select annnotation classes we want from annotatr (can be user-customized)
myAnnots <- c('_genes_promoters', '_genes_5UTRs', '_genes_exons',
'_genes_3UTRs','_cpg_islands')Now let’s annotate for our patient of interest:
#Read in patient 72 Grange data for annotation
dmrs.gr<-subset(combined_files,
sampleID == "vignette_borealis_patient_72")
# Annotate using annotatr
myAnnots <- paste(genome.version,myAnnots,sep="")
annots.all.gr <- annotatr::build_annotations(genome = genome.version,
annotations = myAnnots)
allAnnot <- annotatr::annotate_regions(regions=dmrs.gr,
annotations=annots.all.gr, ignore.strand=TRUE, minoverlap=0)Now let’s look at that most significant site again - with annotations.
# Extract the annotated site with the smallest p-value
subset(allAnnot, pVal == min(allAnnot$pVal))$annot## GRanges object with 9 ranges and 5 metadata columns:
## seqnames ranges strand | id tx_id
## <Rle> <IRanges> <Rle> | <character> <character>
## [1] chr14 24779705-24780704 + | promoter:51107 uc001wos.3
## [2] chr14 24779855-24780576 - | 5UTR:72873 uc001woo.3
## [3] chr14 24779855-24780576 - | 5UTR:72876 uc001wop.3
## [4] chr14 24779861-24781259 + | exon:468887 uc001wor.3
## [5] chr14 24779805-24781259 + | exon:468889 uc010alo.3
## [6] chr14 24779871-24780947 + | exon:468890 uc021rrp.1
## [7] chr14 24779855-24780576 - | exon:481685 uc001woo.3
## [8] chr14 24779855-24780576 - | exon:481692 uc001wop.3
## [9] chr14 24779875-24780932 * | island:17101 <NA>
## gene_id symbol type
## <character> <character> <character>
## [1] 1241 LTB4R hg19_genes_promoters
## [2] 27141 CIDEB hg19_genes_5UTRs
## [3] 27141 CIDEB hg19_genes_5UTRs
## [4] 56413 LTB4R2 hg19_genes_exons
## [5] 56413 LTB4R2 hg19_genes_exons
## [6] <NA> <NA> hg19_genes_exons
## [7] 27141 CIDEB hg19_genes_exons
## [8] 27141 CIDEB hg19_genes_exons
## [9] <NA> <NA> hg19_cpg_islands
## -------
## seqinfo: 93 sequences (1 circular) from hg19 genomeThis CpG site overlaps multiple genes and features in Patient 72 but perhaps most interesting is the LTB4R promoter, since promoters are so strongly linked to control of gene expression.
Let’s use a handy Borealis plotting function to investigate this site further.
# Use Borealis plotting function to investigate this site further
plotCpGsite("chr14:24780288", sampleOfInterest="patient_72",
modelFile=file.path(outdir,"vignette_CpG_model_chr14.csv"),
methCountFile=file.path(outdir,
"vignette_CpG_model_rawMethCount_chr14.tsv"),
totalCountFile=file.path(outdir,
"vignette_CpG_model_rawTotalCount_chr14.tsv"))## $`chr14:24780288`
Figure 1: plotCpGsite function demo
The graph above shows us the expected methylation profile under our model, and how far our site deviates from the expected methylation profile. As you can see this site stands out. It is hypermethylated.
Now what if we wanted to know if this site is surrounded by other hypermethylated CpG sites? We could do this using the data we have already generated, but an advanced approach that will be useful for further analysis is to summarize all our data across epigenetic features.
That’s what we will do next.
4.4 Summarizing single-site data across epigenetic features
First lets set a p-value we will use to determine significance in our feature-summarized data. In our summarization we will make use of corrected p-values and ignore uncorrected p-values. Your analysis can optionally use an alternative approach.
padjThresh <- 0.05Okay - now let’s do the summarization. You can use this code as-is or easily customize it to create summary metrics specific to your application.
# Calculate how may CpGs per annotatr feature and store in dataframe
allAnnot <- as.data.frame(allAnnot)
featureids <- allAnnot$annot.id
featurecnts <- as.data.frame(table(featureids))
colnames(featurecnts) <- c("annot.id", "NoSites")Since the inputs are individual CpGs we are just creating a count of how many individual sites exist in every unique feature. Let’s look:
head(featurecnts)## annot.id NoSites
## 1 3UTR:40890 4
## 2 3UTR:40894 2
## 3 5UTR:70956 7
## 4 5UTR:70963 7
## 5 5UTR:72870 10
## 6 5UTR:72873 38Okay, now we want to figure out how many of the CpGs in each feature have a significant corrected p-value, meaning they appear abnormal. We will subset the sample data on the basic of the p-value and create counts,
# Calculate how many sites per feature pass p-value threshold
# Add data to new summary dataframe
signifonly <- subset(allAnnot, pAdj<=padjThresh)
signifonly <- signifonly$annot.id
signifonlycnt <- as.data.frame(table(signifonly))
colnames(signifonlycnt) <- c("annot.id", "signifCount")
featurecnts <- merge(featurecnts, signifonlycnt, by.x="annot.id",
by.y="annot.id", all.x=TRUE)Now for convenience let’s figure out what fraction of sites in each feature are significant.
# What fraction of sites per feature pass p-value threshold?
featurecnts$fractionSignif <- featurecnts$signifCount/featurecnts$NoSitesAnd finally lets merge the new and orignal data for output:
# Let's combine the data for final output
locations <- subset(allAnnot, select=c("annot.id", "annot.seqnames",
"annot.start", "annot.end"))
featurecnts <- merge(unique(locations), featurecnts, by="annot.id")
genemap <- unique(cbind(allAnnot$annot.symbol, allAnnot$annot.id,
allAnnot$annot.tx_id,allAnnot$annot.width,
allAnnot$sampleID))
colnames(genemap) <- c("annot.symbol", "annot.id", "annot.tx_id", "annot.width",
"SampleID")
summarized <- merge(featurecnts, genemap, by="annot.id")
summarized$signifCount[is.na(summarized$signifCount)] <- 0
summarized$fractionSignif[is.na(summarized$fractionSignif)] <- 0Now let’s have a look at that promoter region!
# Select the LTB4R promoter region
subset(summarized, select=c(annot.symbol, NoSites, signifCount, fractionSignif),
(annot.symbol=="LTB4R" & grepl("promoter", annot.id)))## annot.symbol NoSites signifCount fractionSignif
## 80 LTB4R 47 37 0.787You can see that the outlier site from earlier isn’t alone - multiple other sites are affected in the same region. In fact, in that LTB4R promoter 78.7% of the CpG sites we have data for are called significant by Borealis! That’s 37 CpG sites and all are called hypermethylated. This could be a very relevant finding. Maybe you can now used the code provided to look at the data in different ways. What other genes overlap this region? What features are affected for each?
Thank you for trying Borealis and for working through this vignette. We hope it was helpful. This is as far as we will go but hopefully you have all you need to expand the analysis and apply it to your own data!
Good luck!
Session info
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] org.Hs.eg.db_3.16.0
## [2] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
## [3] GenomicFeatures_1.50.3
## [4] AnnotationDbi_1.60.0
## [5] GenomicRanges_1.50.2
## [6] GenomeInfoDb_1.34.5
## [7] IRanges_2.32.0
## [8] S4Vectors_0.36.1
## [9] borealis_1.2.2
## [10] Biobase_2.58.0
## [11] BiocGenerics_0.44.0
## [12] BiocStyle_2.26.0
##
## loaded via a namespace (and not attached):
## [1] snow_0.4-4 AnnotationHub_3.6.0
## [3] BiocFileCache_2.6.0 plyr_1.8.8
## [5] splines_4.2.2 gamlss_5.4-10
## [7] BiocParallel_1.32.5 ggplot2_3.4.0
## [9] digest_0.6.31 foreach_1.5.2
## [11] htmltools_0.5.4 magick_2.7.3
## [13] fansi_1.0.3 magrittr_2.0.3
## [15] memoise_2.0.1 BSgenome_1.66.1
## [17] doParallel_1.0.17 tzdb_0.3.0
## [19] limma_3.54.0 Biostrings_2.66.0
## [21] readr_2.1.3 matrixStats_0.63.0
## [23] R.utils_2.12.2 prettyunits_1.1.1
## [25] colorspace_2.0-3 blob_1.2.3
## [27] rappdirs_0.3.3 gamlss.data_6.0-2
## [29] xfun_0.36 dplyr_1.0.10
## [31] crayon_1.5.2 RCurl_1.98-1.9
## [33] jsonlite_1.8.4 annotatr_1.24.0
## [35] bsseq_1.34.0 survival_3.4-0
## [37] iterators_1.0.14 glue_1.6.2
## [39] gtable_0.3.1 zlibbioc_1.44.0
## [41] XVector_0.38.0 DelayedArray_0.24.0
## [43] Rhdf5lib_1.20.0 HDF5Array_1.26.0
## [45] scales_1.2.1 DBI_1.1.3
## [47] Rcpp_1.0.9 xtable_1.8-4
## [49] progress_1.2.2 DSS_2.46.0
## [51] bit_4.0.5 httr_1.4.4
## [53] ellipsis_0.3.2 farver_2.1.1
## [55] pkgconfig_2.0.3 XML_3.99-0.13
## [57] R.methodsS3_1.8.2 sass_0.4.4
## [59] dbplyr_2.2.1 locfit_1.5-9.6
## [61] utf8_1.2.2 labeling_0.4.2
## [63] reshape2_1.4.4 tidyselect_1.2.0
## [65] rlang_1.0.6 later_1.3.0
## [67] munsell_0.5.0 BiocVersion_3.16.0
## [69] tools_4.2.2 cachem_1.0.6
## [71] cli_3.5.0 generics_0.1.3
## [73] RSQLite_2.2.20 evaluate_0.19
## [75] stringr_1.5.0 fastmap_1.1.0
## [77] yaml_2.3.6 knitr_1.41
## [79] bit64_4.0.5 purrr_1.0.0
## [81] KEGGREST_1.38.0 nlme_3.1-161
## [83] sparseMatrixStats_1.10.0 mime_0.12
## [85] R.oo_1.25.0 xml2_1.3.3
## [87] biomaRt_2.54.0 compiler_4.2.2
## [89] filelock_1.0.2 curl_4.3.3
## [91] png_0.1-8 interactiveDisplayBase_1.36.0
## [93] tibble_3.1.8 bslib_0.4.2
## [95] stringi_1.7.8 highr_0.10
## [97] lattice_0.20-45 Matrix_1.5-3
## [99] permute_0.9-7 vctrs_0.5.1
## [101] pillar_1.8.1 lifecycle_1.0.3
## [103] rhdf5filters_1.10.0 BiocManager_1.30.19
## [105] jquerylib_0.1.4 data.table_1.14.6
## [107] cowplot_1.1.1 bitops_1.0-7
## [109] httpuv_1.6.7 rtracklayer_1.58.0
## [111] R6_2.5.1 BiocIO_1.8.0
## [113] bookdown_0.31 promises_1.2.0.1
## [115] codetools_0.2-18 MASS_7.3-58.1
## [117] gtools_3.9.4 assertthat_0.2.1
## [119] rhdf5_2.42.0 SummarizedExperiment_1.28.0
## [121] rjson_0.2.21 withr_2.5.0
## [123] regioneR_1.30.0 GenomicAlignments_1.34.0
## [125] Rsamtools_2.14.0 GenomeInfoDbData_1.2.9
## [127] parallel_4.2.2 hms_1.1.2
## [129] grid_4.2.2 rmarkdown_2.19
## [131] DelayedMatrixStats_1.20.0 MatrixGenerics_1.10.0
## [133] shiny_1.7.4 gamlss.dist_6.0-5
## [135] restfulr_0.0.15Oliver, Gavin R., Garrett Jenkinson, Rory J. Olson, Laura E. Shultz-Rogers, and Eric W. Klee. 2022. “Detection of Outlier Methylation from Bisulfite Sequencing Data with Novel Bioconductor Package Borealis.” bioRxiv. https://doi.org/10.1101/2022.05.19.492700.