Citation

If you use borealis in your research, please cite Oliver et al. (2022).

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] 20

list_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  24793719      * |        12        12  0.578665 2.49863e-16
##             pVal    isHypo   effSize
##        <numeric> <logical> <numeric>
##   [1] 0.00281942     FALSE  0.421335
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

# 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  24793719      * |        12        12  0.578665 2.49863e-16
##             pVal    isHypo   effSize               sampleID      pAdj
##        <numeric> <logical> <numeric>            <character> <numeric>
##   [1] 0.00281942     FALSE  0.421335 vignette_borealis_pa..  0.552606
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

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.

# 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 rows of data in combined GRange object?
length(combined_files)

## [1] 3709

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

#Number of unique samples
length(unique(combined_files_df$sampleID))

## [1] 20

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

## [1] 241

#generate summaries for mu and theta
summary(mu_theta$mu)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.028   0.068   0.122   0.332   0.780   0.921

summary(mu_theta$theta)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.0144  0.0243  0.1417

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

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

#Summarize pvals
summary(depth_pvals_eff$pVal)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.001   0.366   0.608   0.602   0.874   1.000

#Summarize corrected pvals
summary(depth_pvals_eff$pAdj)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.019   1.000   1.000   0.975   1.000   1.000

#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.043   0.264   0.560   1.000

#Summarize effect size
summary(depth_pvals_eff$effSize)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -0.483  -0.069  -0.040  -0.017   0.021   0.595

# 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
## 4 vignette_borealis_patient_72 51

# 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
## 365    chr14 24780569 24780569     1      * 17 28 0.117 0.0567 0.000592  FALSE
##     effSize                     sampleID  pAdj
## 365    0.49 vignette_borealis_patient_72 0.019

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 genome

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.

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

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.

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

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

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,

# 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$NoSites

And 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)] <- 0

Now 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
## 85        LTB4R      45          37          0.822

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!