4.2.1 Methylation data

The BPRMeth package provides methods for reading files for methylation, annotation and expression data, however they do not cover all possible data formats available for describing biological datasets. The user can implement his own methods for reading files with different formats, provided that he can create an object similar to what is described below. First, we load the preprocessed methylation data as follows:

library(BPRMeth)  # Load package
met_region <- encode_met

The met_region object contains the following important slots:

• met: A list containing the methylation region, where each element of the list is a different genomic region. More specifically, each methylation promoter region is an \(L_{i} \times 3\) dimensional matrix, where \(L_{i}\) denotes the number of CpGs found in region \(i\). The columns contain the following information:
  • 1st column: Contains the locations of CpGs relative to central location. Note that the actual locations are scaled to the (-1, 1) region.
  • 2nd column: Contains the total reads of each CpG.
  • 3rd column: Contains the methylated reads each CpG.
• anno: Corresponding annotation data for each genomic region as GRanges object. The anno slot is only required for downstream analysis, e.g. to predict gene expression levels from methylation profiles, since we need to map the genomic regions with the gene IDs.

Now let’s observe the format of the met_region object. For example, we can access the \(20^{th}\) promoter methylation region as follows

head(met_region$met[[20]])

##         [,1] [,2] [,3]
## [1,] -0.8787   42    1
## [2,] -0.8779   42    0
## [3,] -0.8770   42    2
## [4,] -0.8759   42    0
## [5,] -0.8754   42    0
## [6,] -0.8667    9    1

Below are the annotation data

head(met_region$anno)

## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames          ranges strand |              id    centre
##          <Rle>       <IRanges>  <Rle> |     <character> <numeric>
##   [1]    chr12   562529-576529      + | ENSG00000139044    569529
##   [2]    chr12   745147-759147      + | ENSG00000177406    752147
##   [3]    chr12 1051888-1065888      - | ENSG00000002016   1058888
##   [4]    chr12 1696331-1710331      - | ENSG00000171823   1703331
##   [5]    chr12 2020870-2034870      - | ENSG00000151062   2027870
##   [6]    chr12 2897118-2911118      + | ENSG00000004478   2904118
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

and the total number of genomic regions are

NROW(met_region$anno)

## [1] 339

4.2.2 Read methylation and annotation data (Do not run)

In the previous section for simplicity we provided a pre-processed object which we then can directly use for inferring methylation profiles. The steps required to create this object are really simple and follow the diagram in Figure 2. First we need to read the methylation data using the read_met function. Also we need read annotation data using the read_anno file. Note that the annotation file can contain any genomic context: from promoters and gene bodies to enhancers, Nanog regulatory regions and CTCF regions; hence BPRMeth can be used for a plethora of analyses that want to take spatial genomic correlations into account. Finally, the create_region_object will create the methylation regions object which is the main object for storing methylation data. The following lines of code provide an example of creating this object ( Do not run )

# Read bulk methylation data
met_dt <- read_met(file = "met_file", type = "bulk_seq")
# Read annotation file, which will create annotation regions as well
anno_dt <- read_anno(file = "anno_file")
# Create methylation regions that have CpG coverage
met_region <- create_region_object(met_dt = met_dt, anno_dt = anno_dt)

4.2.3 Expression data

If the goal is to relate methylation profiles to gene expression levels, we need also to store expression data. Again, here we have pre-processed expression data which can be loaded as follows

expr <- encode_expr

The structure of the expression data is rather simple, it is a data.table with two columns as shown below

head(expr)

##                 id      expr
## 1: ENSG00000230032 -3.321928
## 2: ENSG00000226210  3.833477
## 3: ENSG00000225143 -1.967380
## 4: ENSG00000234465 -3.321928
## 5: ENSG00000206114 -3.321928
## 6: ENSG00000235591 -3.321928

To create the expression data (which are \(log_{2}\) transformed) we could call ( Do not run )

# Read expression data and log2 transform them
expr <- read_expr(file = "expr_file", log2_transf = TRUE)

4.5.1 Predicting gene expression

To quantitatively predict expression at each region, we construct a regression model by taking as input the higher-order methylation features learned from the infer_profiles_vb. In addition to these features, we consider two supplementary sources of information: (1) the goodness of fit in RMSE and (2) the CpG density. For our analysis an SVM regression model is considered. We will use \(70\%\) of the data for training, and we will test the model’s performance on the remaining \(30\%\).

All the aforementioned steps are assembled in the predict_expr wrapper function:

# Set seed for reproducible results
set.seed(22)
# Perform predictions using methylation profiles
pred_profile <- predict_expr(prof_obj = fit_profiles, expr = expr,
                             anno = met_region$anno, model_name = "svm",
                             fit_feature = "RMSE", cov_feature = TRUE, 
                             is_summary = FALSE)
# Perform predictions using mean methylation rate
pred_mean <- predict_expr(prof_obj = fit_mean, expr = expr, 
                          anno = met_region$anno, model_name = "svm",
                          is_summary = FALSE)

We can now compare the Pearson’s correlation coefficient \(r\) for both models and observe that the higher-order methylation features achieve test correlations twice as large compared to mean methylation rate. Note that someone should use cross-validaion to get a better estimate of the correlation performance.

print(paste("Profile r: ", round(pred_profile$test_errors$pcc, 2)))

## [1] "Profile r:  0.66"

print(paste("Mean r:", round(pred_mean$test_errors$pcc, 2)))

## [1] "Mean r: 0.31"

Below in Figure 4 we show a scatter plot of the predicted and measured expression values for the and of the K562 cell line.

g1 <- plot_predicted_expr(pred_obj = pred_profile, 
                          title = "Methylation profile")
g2 <- plot_predicted_expr(pred_obj = pred_mean, 
                          title = "Methylation rate")
g <- cowplot::plot_grid(g1, g2, ncol = 2, nrow = 1)
print(g)

Figure 4: Relationship between DNA methylation patterns and gene expression
Scatter plots of predicted versus measured (log2-transformed) gene expression values using profiles (left) and rates (right); each shaded blue dot represents a different gene, The red line indicates the linear fit between the predicted and measured expression values.

4.5.2 Clustering methylation profiles

Another application of the BPRMeth model is to use the higher-order methylation features to cluster DNA methylation patterns across genomic regions and examine whether distinct methylation profiles are associated to different gene expression levels. To cluster methylation profiles, a mixture modelling approach is considered and we perform inference either via EM algorithm (MLE estimates) or using the mean-field variational inference (Variational Bayes VB) model. In this example we will use the VB model to perform variational inference.

The BPRMeth package provides the cluster_profiles_vb function for performing Bayesian clustering, where the user needs to provide the number of clusters \(K\), the methylation regions \(X\), the observation model and a basis object. For efficiency we run on 20 VB iterations.

# Set seed for reproducible results
set.seed(12) 
# Create basis object
basis_obj <- create_rbf_object(M = 3)
# Set number of clusters K
K <- 4
# Perform clustering
cl_obj <- cluster_profiles_vb(X = met_region$met, K = K, model = "binomial", 
                              alpha_0 = .5, beta_0 = .1,
                              basis = basis_obj, vb_max_iter = 20)

Figure 5 shows the fitted methylation profiles for each cluster.

cluster_plot <- plot_cluster_profiles(cluster_obj = cl_obj)
print(cluster_plot)

Figure 5: Clustering methylation profiles across promoter-proximal regions