---
title: "Borealis outlier methylation detection"
author:
- name: Gavin Oliver
  affiliation: Mayo Clinic
- name: Garrett Jenkinson
  affiliation: Mayo Clinic
  email: jenkinson.william@mayo.edu
- name: Eric Klee
  affiliation: Mayo Clinic
package: borealis
bibliography: borealis.bib
output:
    BiocStyle::html_document
vignette: |
    %\VignetteIndexEntry{Borealis outlier methylation detection}
    %\VignetteEncoding{UTF-8}
    %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 80
---
```{r, echo=FALSE}
options(digits=3)
```

# 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 `r Githubpkg("/FelixKrueger/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 `r Biocpkg("annotatr")`. Finally users will learn how to
summarize CpG data across epigenetic features.

## Citation

If you use borealis in your research, please cite @borealis.

# Installation

The release version of `r Biocpkg("borealis")` is available via Bioconductor and
can be installed as follows:

```{r, eval=FALSE}
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("borealis")
```

The development version of `r Githubpkg("GarrettJenkinson/borealis")` can be
obtained via the Github repository.

It is easiest to install development versions with the CRANpkg("devtools")
package as follows:

``` r
devtools::install_github('GarrettJenkinson/borealis')
```

Changelogs for development releases will be detailed on GitHub releases.

# 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.

```{r, message=FALSE, warning=FALSE}
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:
```{r}
# Read in the name of all files for analysis
borealis_in <- dir(outdir,pattern="*DMLs.tsv")
length(borealis_in)
```

Looks good? Let's continue!

# Basic post-processing and analysis of Borealis results

## 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:

```{r, message=FALSE, warning=FALSE}
# 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:

```{r}
length(list_object_names)
list_object_names[1]
```

Looks good - now let's have a look at one row to get familiar with the fields
of data:

```{r}
listGRanges[[1]][1,]
```
To 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.

## 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. 

```{r}
# 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:

```{r}
listGRanges[[1]][1,]
```

Looking 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.

```{r}
# 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)
```

### How many CpG sites worth of data do we have across all samples combined?

```{r}
# How many rows of data in combined GRange object?
length(combined_files)
```

Lets create a table to summarize useful metrics about the CpG sites

```{r}
# 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)))
```

### How many unique samples and unique CpG sites are we analyzing data from?

```{r}
#Number of unique samples
length(unique(combined_files_df$sampleID))

#Number of unique CpG sites
nrow(unique(mu_theta))
```

### Distribution of the mean methylation values and variability per CpG Site

```{r}
#generate summaries for mu and theta
summary(mu_theta$mu)
summary(mu_theta$theta)
```

Now let's look at the distribitions of other important values.

```{r}
# Create table of unique positions and depth/p-val/padj for each position in 
# each case
depth_pvals_eff <- unique(combined_files_df)
```

### Summary of read depth distributions

```{r}
#Summarize read depths
summary(depth_pvals_eff$n)
```

### Summary of uncorrected p-values

```{r}
#Summarize pvals
summary(depth_pvals_eff$pVal)
```

### Summary of corrected p-values

```{r}
#Summarize corrected pvals
summary(depth_pvals_eff$pAdj)
```

### Summary of methylation fraction across sites

```{r}
#Summarize fraction of methylation
summary(depth_pvals_eff$x/depth_pvals_eff$n)
```

### Summary of effect sizes

```{r}
#Summarize effect size
summary(depth_pvals_eff$effSize)
```

### 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):

```{r}
# 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%"])
```

We can see that one patient appears to be an extreme outlier.

Let's also have a look at which site is most significant.

```{r}
# What is the most significant CpG site between all samples?
subset(combined_files_df, pVal == min(combined_files$pVal))
```

--------------------------------------------------------------------------------

Depending 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.

## 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
`r Biocpkg("annotatr")` package but you can use your favorite tool instead. The
release version of `r Biocpkg("annotatr")` is available via Bioconductor and can
be installed as follows:

```{r, eval=FALSE}
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 `r Biocpkg("annotatr")`.

Let's will start by defining the annotations we want to utilize.

```{r, message=FALSE, results='hide', warning=FALSE}
#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:

```{r, message=FALSE, results='hide', warning=FALSE}
#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.

```{r}
# Extract the annotated site with the smallest p-value
subset(allAnnot, pVal == min(allAnnot$pVal))$annot
```

This 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.

```{r sg, fig.wide=TRUE, fig.cap="plotCpGsite function demo"}
# 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"))
```

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.

## 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.

```{r}
padjThresh <- 0.05
```

Okay - 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.



```{r}
# 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:

```{r}
head(featurecnts)
```

Okay, 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,

```{r}
# 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.

```{r}
# What fraction of sites per feature pass p-value threshold?
featurecnts$fractionSignif <- featurecnts$signifCount/featurecnts$NoSites
```

And finally lets merge the new and orignal data for output:
```{r}
# 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)] <- 0
```

Now let's have a look at that promoter region!

```{r}
# Select the LTB4R promoter region
subset(summarized, select=c(annot.symbol, NoSites, signifCount, fractionSignif),
        (annot.symbol=="LTB4R" & grepl("promoter", annot.id)))
```

You 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 {.unnumbered}

```{r sessionInfo, echo=FALSE}
sessionInfo()
```
