Introduction

MeDeCom is an R-package for reference-free decomposition of heterogeneous DNA methylation profiles. It uses matrix factorization enhanced by constraints and a specially tailored regularization. MeDeCom represents an input $m\times n$ data matrix ($m$ CpGs measured in $n$ samples) as a product of two other matrices. The first matrix has $m$ rows, just as the input data, but the number of columns is equal to $k$. The columns of this matrix can be interpreted as methylomes of the $k$ unknown cell populations underlying the samples and will be referred to as latent methylation components or LMCs. The second matrix has $k$ rows and $n$ columns, and can be interpreted as a matrix of relative contributions (mixing proportions) of each LMC to each sample.

MeDeCom starts with a set of related DNA methylation profiles, e.g. a series of Infinium microarray measurements from a population cohort, or several bisulfite sequencing-based methylomes. The key requirement is that the input data represents absolute DNA methylation measurements in a population of cells and contains values between 0 and 1. MeDeCom implements an alternating scheme which iteratively updates randomly initialized factor matrices until convergence or until the maximum number of iterations has been reached. This is repeated for multiple random initializations and the best solution is returned.

MeDeCom features two tunable parameters. The first one is the number of LMCs $k$, an approximate choice for which should be known from prior information. To enforce the distribution properties of a methylation profile upon LMCs MeDeCom uses a special for of regularization controlled by the parameter $\lambda$. A typical MeDeCom experiment includes testing a grid of values for $k$ and $\lambda$. For each combination of parameter values MeDeCom estimates a cross-validation error. The latter helps select the optimal number of LMCs and the strength of regularization.

Installation

MeDeCom can be installed directly from github using the package devtools:

devtools::install_github("lutsik/MeDeCom")

The master branch of the GitHub repository only compiles on nix-like platforms with a C++11-compatible compiler are supported. Users on Mac OS might run into compilation errors as OpenMP, which is required for compilation, is disabled by default. To enable OpenMP, the following tutorial can be followed https://mac.r-project.org/openmp/ We created a separate branch for installation of MeDeCom on Windows machines. However, parallel processing options are currently not supported for Windows due to some incompatabilites of the R/Windows connection. Thus, executing MeDeCom on a Windows machine takes sustantially longer. Additionally, we creared a Docker image with MeDeCom, which can be used in case of further installation issues https://hub.docker.com/r/mscherer/medecom.

devtools::install_github("lutsik/MeDeCom",ref="windows")

MeDeCom uses stack model for memory to accelerate factorization for smaller ranks. This requires certain preparation during compilation, therefore, please, note the extended compilation time (15 to 20 minutes).

Data preparation

MeDeCom accepts DNA methylation data in several forms. Preferably the user may load and preprocess the data using a general-purpose DNA methylation analysis package RnBeads. A resulting RnBSet object can be directly supplied to MeDeCom. Alternatively, MeDeCom runs on any matrix of type numeric with valid methylation values.

MeDeCom comes with a small example data set obtained by mixing reference profiles of blood cell methylomes in silico. The example data set can be loaded in a usual way:

## load the package
suppressPackageStartupMessages(library(MeDeCom))
##  load the example data sets
data(example.dataset, package="MeDeCom")
## you should get objects D, Tref and Aref
## in your global R environment
ls()

Loaded numeric matrix D contains 100 in silico mixtures and serves as an example input. Columns of matrix Tref contains the methylomes of 5 blood cell types used to generate the mixtures, while matrix Aref provides the mixing proportions.

## matrix D has dimension 10000x100
str(D)
## matrix Tref has dimension 10000x5
str(Tref)
## matrix Aref has dimension 5x100
str(Aref)

Performing a methylome decomposition experiment

MeDeCom can be run directly on matrix D.

It is crucial to select the values of parameters $k$ and $\lambda$ to test. A choice of $k$ is often dictated by prior knowledge about the methylomes. Precise value of lambda has to be selected for each data set independently. A good start is a logarithmic grid of lambda values. It is important to include $\lambda=0$ into the grid, as this particular case the regularization is effectively absent making MeDeCom similar to other NMF-based deconvolution algorithms.

medecom.result<-runMeDeCom(D, Ks=2:10, lambdas=c(0,10^(-5:-1)))

MeDeCom is based upon an alternating optimization heuristic and requires a lot of computation. The processing of the data matrix can take several hours. One can speed up the run by decreasing the number of cross-validation folds and random initializations, and increasing the number of computational cores.

medecom.result<-runMeDeCom(D, 2:10, c(0,10^(-5:-1)), NINIT=10, NFOLDS=10, ITERMAX=300, NCORES=9)
cat("
[Main:] checking inputs
[Main:] preparing data
[Main:] preparing jobs
[Main:] 3114 factorization runs in total
[Main:] runs 2755 to 2788 complete
[Main:] runs 2789 to 2822 complete
[Main:] runs 2823 to 2856 complete
[Main:] runs 2857 to 2890 complete
......
[Main:] finished all jobs. Creating the object
")
data(example.MeDeComSet)

This can, however, lead to decomposition slightly different from the one presented given below.

The results of a decomposition experiment are saved to an object of class MeDeComSet. The contents of an object can be conveniently displayed using the print functionality.

medecom.result

Exploring the decomposition results

Parameter selection

The first key step is parameter selection. It is important to carefully explore the obtained results and make a decision about the most feasible parameter values, or about extending the parameter value grids to be tested in refinement experiments.

MeDeCom provides a cross-validation error (CVE) for each tested parameter combination.

plotParameters(medecom.result)

A lineplot helping to select parameter $\lambda$ can be produced by specifying a fixed value for $k$:

plotParameters(medecom.result, K=5, lambdaScale="log")

Cross-validation error has a minimum at $\lambda=10^{-2}$ so this value is preferred.

Latent methylation components (LMCs)

A matrix of LMCs can be extracted using getLMCs:

lmcs<-getLMCs(medecom.result, K=5, lambda=0.01)
str(lmcs)

LMCs can be seen as measured methylation profiles of purified cell populations. MeDeCom provides for several visualization methods for LMCs using the function plotLMCs which operates directly on MeDeComSet objects.

Clustering

For instance, standard hierarchical clustering can be visualized using:

plotLMCs(medecom.result, K=5, lambda=0.01, type="dendrogram")

A two-dimensional embedding with MDS is also obtainable:

plotLMCs(medecom.result, K=5, lambda=0.01, type="MDS")

Input data can be included into the MDS plot to enhance the interpretation.

plotLMCs(medecom.result, K=5, lambda=0.01, type="MDS", D=D)

Matching LMCs to reference profiles

In many cases reference methylomes exists, which are relevant for the data set in question. For our example analysis matrix Tref contains the reference type profiles which were in silico mixed. MeDeCom offers several ways to visualize the resulting LMCs together with the reference methylation profiles. The reference methylomes can be included into a joint clustering analysis:

plotLMCs(medecom.result, K=5, lambda=0.01, type="dendrogram", Tref=Tref, center=TRUE)

Furthermore, a similarity matrix of LMCs vs reference profiles can be visualized as a heatmap.

plotLMCs(medecom.result, K=5, lambda=0.01, type="heatmap", Tref=Tref)

Correlation coefficient values and asterisks aid the interpretation. The values are displayed in the cells which contain maximal values column-wise. The asterisks mark cells which have the highest correlation value in the respective rows. Thus, a value with asterisk corresponds to a mutual match, i.e. LMC unambiguously matching a reference profile.

In this example analysis each LMC uniquely matches one of the reference profiles. The matching of

Function matchLMCs offers several methods for matching LMCs to reference profiles.

perm<-matchLMCs(lmcs, Tref)

LMC enrichment analysis

MeDeCom provides functions to perform enrichment analysis on the sites that are particularly hypo-/hypermethylated in an LMC. These sites can then be used for GO and LOLA enrichment analysis. Importantly, genomic annotations of the LMC sites is required to be specified. We thus recommend to use the DecompPipeline package for processing, but the annotation can also be specified manually using a data.frame that looks as follows:

      Chromosome   Start     End Strand CpG GC CGI Relation SNPs
30365       chr1 1036375 1036376      +   2 59        Shelf <NA>
42681       chr1 1184537 1184538      +   2 58     Open Sea <NA>
45091       chr1 1218625 1218626      +   5 66       Island <NA>
51615       chr1 1292773 1292774      +   3 64       Island <NA>
52001       chr1 1295504 1295505      +   9 65        Shore <NA>
52003       chr1 1295507 1295508      +   9 65        Shore <NA>

The required columns are Chromosome, Start, End, and Strand. Using this data.frame (called df in the following), enrichment analysis can be performed using:

lmc.lola.enrichment(medecom.result,anno.data=df,K=5,lambda=0.001,diff.threshold = 0.5, region.type = "tiling")

Please note that df needs to have the same number of rows than the methylation matrix used as input to MeDeCom. CpGs are first aggregated over the region.type specified, then the regions are selected that have a difference larger than diff.threshold. The list of available region types is published here https://rnbeads.org/regions.html.

Mixing proportions

A matrix of mixing proportions is obtained using getProportions:

prop<-getProportions(medecom.result, K=5, lambda=0.001)
str(prop)

Visualization of the complete proportion matrix

A complete matrix of propotions can be visualized as a stacked barplot:

plotProportions(medecom.result, K=5, lambda=0.01, type="barplot")

or a heatmap:

plotProportions(medecom.result, K=5, lambda=0.01, type="heatmap")

The heatmap can be enhanced by clustering the columns:

plotProportions(medecom.result, K=5, lambda=0.01, type="heatmap", heatmap.clusterCols=TRUE)

or adding color code for the samples:

sample.group<-c("Case", "Control")[1+sample.int(ncol(D))%%2]
plotProportions(medecom.result, K=5, lambda=0.01, type="heatmap", sample.characteristic=sample.group)

Visualization of selected LMC proportions

rownames(Aref)<-colnames(Tref)
plotProportions(medecom.result,  K=5, lambda=0.01, type="lineplot", lmc=2, Aref=Aref, ref.profile=2)

Advanced usage

Running MeDeCom on a compute cluster

MeDeCom experiments require a lot of computational time. On the other hand most of the factorization runs are independent and, therefore, can be run in parallel. Thus, a significant speedup can be achieved when running MeDeCom in an HPC environment. MeDeCom can be easily adapted to most of the popular schedulers. There are, however, several prerequisites:

The example below is for the cluster operated by Son of Grid Engine (SoGE). To be able to run on a SoGE cluster MeDeCom needs to know:

These settings should be stored in a list object:

sge.setup<-list(
R_bin_dir="/usr/bin",
host_pattern="*",
mem_limit="5G"
)

This object should be supplied to MeDeCom as the argument cluster.settings. It is also important to specify a valid temporary directory, which is available to all execution nodes.

medecom.result<-runMeDeCom(D, Ks=2:10, lambdas=c(0,10^(-5:-1)), N_COMP_LAMBDA=1, NFOLDS=5, NINIT=10, 
temp.dir="/cluster_fs/medecom_temp",
cluster.settings=sge.setup)

MeDeCom will start the jobs and will periodically monitor the number of remaining ones.

cat("
[Main:] checking inputs
[Main:] preparing data
[Main:] preparing jobs
[Main:] 3114 factorization runs in total
[Main:] 3114 jobs remaining
....
[Main:] finished all jobs. Creating the object
")

R session

Here is the output of sessionInfo() on the system on which this document was compiled:

sessionInfo()


lutsik/MeDeCom documentation built on Feb. 15, 2023, 11:32 a.m.