library(methodical)
library(DESeq2)
library(TumourMethData)methodical facilitates the rapid analysis of the association between DNA
methylation and expression of associated transcripts. It can calculate
correlation values between individual methylation sites (e.g. CpG sites) or the
methylation of wider genomic regions (e.g. expected promoter regions, gene bodies)
and the expression of associated transcripts. We will demonstrate both cases in
this vignette.
# Installing Methodical
if (!require("BiocManager"))
install.packages("BiocManager")
BiocManager::install("methodical")We now first download a RangedSummarizedExperiment with WGBS data for prostate tumour metastases and associated transcript counts from TumourMethData.
# Download mcrpc_wgbs_hg38 from TumourMethDatasets.
mcrpc_wgbs_hg38_chr11 <- TumourMethData::download_meth_dataset(dataset = "mcrpc_wgbs_hg38_chr11")
# Download transcript counts
mcrpc_transcript_counts_file <- TumourMethData::download_rnaseq_dataset(dataset = "mcrpc_wgbs_hg38_chr11")
mcrpc_transcript_counts = data.frame(data.table::fread(mcrpc_transcript_counts_file), row.names = 1)We now load annotation for our transcripts of interest, in this case all protein-coding transcripts annotated by Gencode.
# Download Gencode annotation
library(AnnotationHub)
transcript_annotation <- AnnotationHub()[["AH75119"]]
# Import the Gencode annotation and subset for transcript annotation for non-mitochondrial protein-coding genes
transcript_annotation <- transcript_annotation[transcript_annotation$type == "transcript"]
transcript_annotation <- transcript_annotation[transcript_annotation$transcript_type == "protein_coding"]
# Remove transcript version from ID
transcript_annotation$ID <- gsub("\\.[0-9]*", "", transcript_annotation$ID)
# Filter for transcripts on chromosome 11
transcript_annotation_chr11 <- transcript_annotation[seqnames(transcript_annotation) == "chr11"]
# Filter for transcripts in the counts table
transcript_annotation_chr11 <- transcript_annotation_chr11[transcript_annotation_chr11$ID %in% row.names(mcrpc_transcript_counts)]
# Set transcript ID as names of ranges
names(transcript_annotation_chr11) <- transcript_annotation_chr11$ID
# Get the TSS for each transcript set ID as names
transcript_tss_chr11 <- resize(transcript_annotation_chr11, 1, fix = "start")
rm(transcript_annotation)
# Get the region surrounding each TSS
tss_proximal_regions_chr11 = methodical::expand_granges(transcript_tss_chr11, upstream = 5000, downstream = 5000)Next we will normalize counts data for these transcripts from prostate cancer metastases samples using DESeq2.
# Subset mcrpc_transcript_counts for protein-coding transcripts
mcrpc_transcript_counts <- mcrpc_transcript_counts[transcript_tss_chr11$ID, ]
# Create a DESeqDataSet from Kallisto counts.
mcrpc_transcript_counts_dds <- DESeqDataSetFromMatrix(countData = mcrpc_transcript_counts,
colData = data.frame(sample = names(mcrpc_transcript_counts)), design = ~ 1)
mcrpc_transcript_counts_dds <- estimateSizeFactors(mcrpc_transcript_counts_dds)
mcrpc_transcript_counts_normalized <- data.frame(counts(mcrpc_transcript_counts_dds, normalized = TRUE))We’ll first demonstrate calculating the correlation values between all transcripts
on chromosome 11 and the methylation of all CpG sites within 5 KB of their TSS
and also the significance of these correlation values. We do this with the
calculateMethSiteTranscriptCors function, which takes methylation
RangedSummarizedExperiment and a table with transcript counts as input.
It also takes a GRanges object with the TSS of interest. The arguments
expand_upstream and expand_downstream define the regions around the TSS
where we want to examine CpG sites. We set these to 5KB upstream and downstream.
We can also provide a subset of samples that we want to use to calculate the
correlation values with the samples_subset parameter.
The default behaviour is to use try to use all samples, but we use
common_mcprc_samples since not all samples in mcrpc_wgbs_hg38_chr11 are
also present in mcrpc_transcript_counts_normalized.
Finally, we can control the memory usage and parallelization with the
max_sites_per_chunk and BPPARAM parameters.
calculateMethSiteTranscriptCors aims to keep the memory footprint low by
splitting the provided tss_grinto chunks of closely-situated TSS and
processing these chunks one at a time. The chunks are limited to that they
contain an approximate maximum number of CpGs within expand_upstream and
expand_downstream of the constituent TSS.
max_sites_per_chunk controls this upper limit of CpG that are read into memory
at once and has a default value of floor(62500000/ncol(meth_rse)) using
approximately 500 MB of RAM.
BPPARAM controls the number of TSS within a chunk processed in parallel. As
the TSS belong to the same chunk, several cores can be used simultaneously
without drastically increasing RAM consumption. The following example
should just take a few minutes.
# Find samples with both WGBS and RNA-seq count data
common_mcprc_samples <- intersect(names(mcrpc_transcript_counts), colnames(mcrpc_wgbs_hg38_chr11))
# Create a BPPARAM class
BPPARAM = BiocParallel::bpparam()
# Calculate CpG methylation-transcription correlations for all TSS on chromosome 11
system.time({transcript_cpg_meth_cors_mcrpc_samples_5kb <- calculateMethSiteTranscriptCors(
meth_rse = mcrpc_wgbs_hg38_chr11,
transcript_expression_table = mcrpc_transcript_counts_normalized,
samples_subset = common_mcprc_samples,
tss_gr = transcript_tss_chr11, tss_associated_gr = tss_proximal_regions_chr11,
cor_method = "pearson", max_sites_per_chunk = NULL, BPPARAM = BPPARAM)})As an example, we’ll examine the CpG methylation transcription correlations for a TSS of GSTP1, a gene known to be important to prostate cancer. We’ll use the TSS for the ENST00000398606 transcript, which is the ENSEMBL canonical and MANE select transcript for GSTP1.
We see that we have five columns in the results table: the name of the CpG site, the name of the associated transcript, the correlation value, the p-value of the correlation and the distance of the CpG site to the TSS of the transcript.
The location of the TSS associated with the results is stored as a GRanges as an
attribute called tss_range, making it easy to retrieve this information.
# Extract correlation results for ENST00000398606 transcript
gstp1_cpg_meth_transcript_cors <- transcript_cpg_meth_cors_mcrpc_samples_5kb[["ENST00000398606"]]
# Examine first few rows of the results
head(gstp1_cpg_meth_transcript_cors)
# Extract TSS for the transcript
attributes(gstp1_cpg_meth_transcript_cors)$tss_rangePlotting the methylation-transcription correlation values around a TSS enable
us to visualize how the relationship between DNA methylation and transcription
changes around the TSS. We can use the plotMethSiteCorCoefs function to
plot values for methylation sites, such as their correlation values or
methylation levels. It takes a data.frame as input which should have a column
called “meth_site” which gives the location of CpG sites as a character string
which can be coerced to a GRanges (e.g. “seqname:start”). We also provide the
name a numeric column which contains the values we want to plot, which is “cor”
here.
Values are displayed on the y-axis and the genomic location
on the x-axis. This location can be either the actual genomic location on the
relevant sequence or else the distance relative to the TSS, depending on
whether reference_tss is set to TRUE or FALSE.
# Plot methylation-transcription correlation values for GSTP1 showing
# chromosome 11 position on the x-axis
meth_cor_plot_absolute_distance <- plotMethSiteCorCoefs(
meth_site_cor_values = gstp1_cpg_meth_transcript_cors,
reference_tss = FALSE,
ylabel = "Correlation Value",
title = "Correlation Between GSTP1 Transcription and CpG Methylation")
print(meth_cor_plot_absolute_distance)
# Plot methylation-transcription correlation values for GSTP1 showing
# distance to the GSTP1 on the x-axis
meth_cor_plot_relative_distance <- plotMethSiteCorCoefs(
meth_site_cor_values = gstp1_cpg_meth_transcript_cors,
reference_tss = TRUE,
ylabel = "Correlation Value",
title = "Correlation Between GSTP1 Transcription and CpG Methylation")
print(meth_cor_plot_relative_distance)We see that there is a group of negatively correlated CpGs immediately surrounding the TSS and positively correlated CpGs further away, especially between 2,500 and 1,250 base pairs upstream.
The returned plots are ggplot objects and can be manipulated in any way that a regular ggplot object can.
# Add a dashed line to meth_cor_plot_relative_distance where the correlation is 0
meth_cor_plot_relative_distance + ggplot2::geom_hline(yintercept = 0, linetype = "dashed")We’ll now demonstrate calculating the correlation values between transcript
expression and methylation of wider regions associated them instead of just
single CpG sites, e.g. promoters, gene bodies or enhancers. We use the
calculateRegionMethylationTranscriptCors function for this, which takes a
RangedSummarizedExperiment with methylation values and a table with transcript
expression values, like calculateMethSiteTranscriptCors. However,
calculateRegionMethylationTranscriptCors takes a GRanges object with wider
regions which are associated with the transcripts. More than one region may be
associated with a given transcript, for example we could analyse the correlation
between transcript expression and methylation of its exons separately for each
exon.
Here we’ll investigate the relationship between gene body methylation and transcription for all transcripts on chromosome 11.
# Calculate gene body methylation-transcription correlations for all TSS on chromosome 11
system.time({transcript_body_meth_cors_mcrpc_samples <- calculateRegionMethylationTranscriptCors(
meth_rse = mcrpc_wgbs_hg38_chr11,
transcript_expression_table = mcrpc_transcript_counts_normalized,
samples_subset = common_mcprc_samples,
genomic_regions = transcript_annotation_chr11,
genomic_region_transcripts = transcript_annotation_chr11$ID,
region_methylation_summary_function = colMeans,
cor_method = "pearson", BPPARAM = BPPARAM)})
# Show the top of the results table
head(transcript_body_meth_cors_mcrpc_samples)
# We'll filter for significant correlation values and plot their distribution
transcript_body_meth_cors_significant <-
dplyr::filter(transcript_body_meth_cors_mcrpc_samples, q_val < 0.05)
ggplot(transcript_body_meth_cors_significant, aes(x = cor)) +
geom_histogram() + theme_bw() +
scale_y_continuous(expand = expansion(mult = c(0, 0.2), add = 0)) +
labs(x = "Correlation Values", y = "Number of Significant Correlation Values")
Within minutes, we calculated all the correlation values for the methylation of CpG sites within 5KB of TSS on chromosome 11 and expression of the associated transcripts as well as the p-values for the correlations. We then analysed the relationship between mean gene body methylation and transcript expression, first calculating the mean methylation of CpGs overlapping gene bodies on chromosome 11 and then computing the correlation values, again in short time using just a few cores.
Thus, in under an hour we could easily analyse the correlations for all TSS in the genome on a desktop computer or even a laptop using just a few cores.
sessionInfo()