In samarafk/MLML2R: Maximum Likelihood Estimation of DNA Methylation and Hydroxymethylation Proportions
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
BiocStyle::markdown()
Introduction
In a given CpG site from a single cell we will either have a $C$ or a $T$ after DNA processing conversion methods, with a different interpretation for each of the available methods. This is a binary outcome and we assume a Binomial model and use the maximum likelihood estimation method to obtain the estimates for hydroxymethylation and methylation proportions.
$T$ reads are referred to as converted cytosine and $C$ reads are referred to as unconverted cytosine. Conventionally, $T$ counts are also referred to as unmethylated counts, and $C$ counts as methylated counts. In case of Infinium Methylation arrays, we have intensities representing the methylated (M) and unmethylated (U) channels that are proportional to the number of unconverted and converted cytosines ($C$ and $T$, respectively). The most used summary from these experiments is the proportion $\beta=\frac{M}{M+U}$, commonly referred to as \textit{beta-value}, which reflects the methylation level at a CpG site. Naïvely using the difference between betas from BS and oxBS as an estimate of 5-hmC (hydroxymethylated cytosine), and the difference between betas from BS and TAB as an estimate of 5-mC (methylated cytosine) can many times provide negative proportions and instances where the sum of 5-C (unmodified cytosine), 5-mC and 5-hmC proportions is greater than one due to measurement errors.
\Biocpkg{MLML2R} package allows the user to jointly estimate hydroxymethylation and methylation consistently and efficiently.
The function \Rfunction{MLML} takes as input the data from the different methods and returns the estimated proportion of methylation, hydroxymethylation and unmethylation for a given CpG site. Table 1 presents the arguments of the \Rfunction{MLML} and Table 2 lists the results returned by the function.
The function assumes that the order of the rows and columns in the input matrices are consistent. In addition, all the input matrices must have the same dimension. Usually, rows represent CpG loci and columns are the samples.
Arguments | Description
--------------------|------------------------
\Robject{G.matrix} | Unmethylated channel (Converted cytosines/ T counts) from TAB-conversion (reflecting 5-C + 5-mC).
\Robject{H.matrix} | Methylated channel (Unconverted cytosines/ C counts) from TAB-conversion (reflecting True 5-hmC).
\Robject{L.matrix} | Unmethylated channel (Converted cytosines/ T counts) from oxBS-conversion (reflecting 5-C + 5-hmC).
\Robject{M.matrix} | Methylated channel (Unconverted cytosines/ C counts) from oxBS-conversion (reflecting True 5-mC).
\Robject{T.matrix} | Methylated channel (Unconverted cytosines/ C counts) from standard BS-conversion (reflecting 5-mC+5-hmC).
\Robject{U.matrix} | Unmethylated channel (Converted cytosines/ T counts) from standard BS-conversion (reflecting True 5-C).
: MLML function and random variable notation.
Value | Description
------|------------------
\Robject{mC} | maximum likelihood estimate for the 5-mC proportion
\Robject{hmC} | maximum likelihood estimate for the 5-hmC proportion
\Robject{C} | maximum likelihood estimate for the 5-mC proportion
\Robject{methods} | the conversion methods used to produce the MLE
: Results returned from the \Rfunction{MLML} function
Worked examples
Simulated data
\Biocpkg{MLML2R} includes small example datasets for illustration.
We simulated counts from Binomial model with true proportions of methylation, hydroxymethylation and unmethylated being $0.3$, $0.2$, and $0.5$, respectively. For instance, \verb|MethylatedBS_sim| is a matrix of simulated counts from BS corresponding to 100 CpGs and 2 samples. Similarly we simulated matrices with methylated and unmethylated counts for all the three methods: BS, oxBS and TAB. The rows and columns in the input matrices are consistent.
BS and oxBS methods
Load the package:
library(MLML2R)
When only two methods are available, the default option returns the exact constrained maximum likelihood estimates using the the pool-adjacent-violators algorithm (PAVA) [@ayer1955].
results_exactBO1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim)
Maximum likelihood estimate via EM-algorithm approach [@Qu:MLML] is obtained with the option \verb|iterative=TRUE|. In this case, the default (or user specified) \verb|tol| is considered in the iterative method.
results_emBO1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,iterative=TRUE)
When only two methods are available, we highly recommend the default option \Rcode{iterative=FALSE} since the difference in the estimates obtained via EM and exact constrained is very small, but the former requires more computational effort:
all.equal(results_emBO1$hmC,results_exactBO1$hmC,scale=1)
library(microbenchmark)
mbmBO1 = microbenchmark(
EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim),
EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
iterative=TRUE),
times=10)
mbmBO1
BS and TAB methods
Using PAVA:
results_exactBT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Using EM-algorithm:
results_emBT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
all.equal(results_emBT1$hmC,results_exactBT1$hmC,scale=1)
mbmBT1 = microbenchmark(
EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim),
EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,
iterative=TRUE),
times=10)
mbmBT1
oxBS and TAB methods
Using PAVA:
results_exactOT1 <- MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Using EM-algorithm:
results_emOT1 <- MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
all.equal(results_emOT1$hmC,results_exactOT1$hmC,scale=1)
mbmOT1 = microbenchmark(
EXACT = MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim),
EM = MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,
iterative=TRUE),
times=10)
mbmOT1
BS, oxBS and TAB methods
When data from the three methods are available, the default otion in the \Rfunction{MLML} function returns the constrained maximum likelihood estimates using an approximated solution for Lagrange multipliers method.
results_exactBOT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Maximum likelihood estimate via EM-algorithm approach [@Qu:MLML] is obtained with the option \verb|iterative=TRUE|. In this case, the default (or user specified) \verb|tol| is considered in the iterative method.
results_emBOT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
We recommend the default option \Rcode{iterative=FALSE} since the difference in the estimates obtained via EM and the approximate exact constrained is very small, but the former requires more computational effort:
all.equal(results_emBOT1$hmC,results_exactBOT1$hmC,scale=1)
mbmBOT1 = microbenchmark(
EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim),
EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim,
L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,
G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,
iterative=TRUE),
times=10)
mbmBOT1
Publicly available data: GSE63179
We will use the dataset from @10.1371/journal.pone.0118202, which consists of eight DNA samples from the same DNA source treated with oxBS-BS and hybridized to the Infinium 450K array.
When data is obtained through Infinium Methylation arrays, we recommend the use of the \Biocpkg{minfi} package [@minfi], a well-established tool for reading, preprocessing and analysing DNA methylation data from these platforms. Although our example relies on \Biocpkg{minfi} and other \Bioconductor{} tools, \Biocpkg{MLML2R} does not depend on any packages. Thus, the user is free to read and preprocess the data using any software of preference and then import the intensities (or $T$ and $C$ counts) for the methylated and unmethylated channel (or converted and uncoverted cytosines) into \R{} in matrix format.
To start this example we will need the following packages:
library(minfi)
library(GEOquery)
It is usually best practice to start the analysis from the raw data, which in the case of the 450K array is a \verb|.IDAT| file.
The raw files are deposited in GEO and can be downloaded by using the \Rfunction{getGEOSuppFiles}. There are two files for each replicate, since the 450k array is a two-color array. The \verb|.IDAT| files are downloaded in compressed format and need to be uncompressed before they are read by the \Rfunction{read.metharray.exp} function.
getGEOSuppFiles("GSE63179")
untar("GSE63179/GSE63179_RAW.tar", exdir = "GSE63179/idat")
list.files("GSE63179/idat", pattern = "idat")
files <- list.files("GSE63179/idat", pattern = "idat.gz$", full = TRUE)
sapply(files, gunzip, overwrite = TRUE)
The \verb|.IDAT| files can now be read:
rgSet <- read.metharray.exp("GSE63179/idat")
To access phenotype data we use the \Rfunction{pData} function. The phenotype data is not yet available from the \Robject{rgSet}.
pData(rgSet)
In this example the phenotype is not really relevant, since we have only one sample: male, 25 years old. What we do need is the information about the conversion method used in each replicate: BS or oxBS. We will access this information automatically from GEO:
geoMat <- getGEO("GSE63179")
pD.all <- pData(geoMat[[1]])
pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1.1",
"characteristics_ch1.2","characteristics_ch1.3")]
pD
This phenotype data needs to be merged into the methylation data. The following commands guarantee we have the same replicate identifier in both datasets before merging.
sampleNames(rgSet) <- sapply(sampleNames(rgSet),function(x)
strsplit(x,"_")[[1]][1])
rownames(pD) <- pD$geo_accession
pD <- pD[sampleNames(rgSet),]
pData(rgSet) <- as(pD,"DataFrame")
rgSet
The \Robject{rgSet} object is a class called \Rclass{RGChannelSet} used for two color data (green and a red channel). The input in the \Rfunction{MLML} funcion is \Rclass{MethylSet}, which contains the methylated and unmethylated signals. The most basic way to construct a \Rclass{MethylSet} is using the function \Rfunction{preprocessRaw}. Here we chose the function \Rfunction{preprocessNoob} [@noob] for background correction and construction of the \Rclass{MethylSet}.
MSet.noob<- preprocessNoob(rgSet)
After the preprocessed steps we can use \Rfunction{MLML} from the \Biocpkg{MLML2R} package.
The BS replicates are in columns 1, 3, 5, and 6. The remaining columns are from the oxBS treated replicates.
MethylatedBS <- getMeth(MSet.noob)[,c(1,3,5,6)]
UnMethylatedBS <- getUnmeth(MSet.noob)[,c(1,3,5,6)]
MethylatedOxBS <- getMeth(MSet.noob)[,c(7,8,2,4)]
UnMethylatedOxBS <- getUnmeth(MSet.noob)[,c(7,8,2,4)]
In this example we only have two methods, therefore we can choose between the EM-algorithm and the exact constrained maximum likelihood estimates (using PAVA).
Estimates via the EM-algorithm:
results_emPD1 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS,
L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS,
iterative = TRUE)
The exact constrained MLE (using PAVA):
results_exactPD1 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS,
L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS)
Publicly available data: GSE73895
We will use the dataset from @pmid27886174, which consists of 30 DNA samples from glioblastoma tumors treated with oxBS-BS and hybridized to the Infinium 450K array.
The raw files are deposited in GEO and can be downloaded and read into \R{} by doing:
getGEOSuppFiles("GSE73895")
untar("GSE73895/GSE73895_RAW.tar", exdir = "GSE73895/idat")
idatFiles <- list.files("GSE73895/idat", pattern = "idat.gz$", full = TRUE)
sapply(idatFiles, gunzip, overwrite = TRUE)
rgSet <- read.metharray.exp("GSE73895/idat")
We need to identify the samples from different methods: BS-conversion, oxBS-conversion. We obtain this information from GEO:
geoMat <- getGEO("GSE73895")
pD.all <- pData(geoMat[[1]])
pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1",
"characteristics_ch1.2","characteristics_ch1.3")]
head(pD)
Keeping only some of the variables from phenotype data:
names(pD)[c(1,3,4,5)] <- c("method","gender","survival_months","age_years")
pD$gender <- sub("^gender: ", "", pD$gender)
pD$age_years <- as.numeric(sub("^subject age: ", "", pD$age_years))
pD$survival_months <- as.numeric(sapply(pD$survival_months, function(x)
strsplit(as.character(x),":")[[1]][2]))
pD$method <- sapply(pD$method, function(x) strsplit(as.character(x),"_")[[1]][3])
We now need to merge this pheno data into the methylation data. The following are commands to make sure we have the same row identifier in both datasets before merging.
sampleNames(rgSet) <- sapply(sampleNames(rgSet),function(x)
strsplit(x,"_")[[1]][1])
rownames(pD) <- pD$geo_accession
pD <- pD[sampleNames(rgSet),]
pData(rgSet) <- as(pD,"DataFrame")
rgSet
Preprocessing and preparing input matrices for \Rfunction{MLML}:
MSet.noob<- preprocessNoob(rgSet)
BS_index <- which(pData(rgSet)$method=="BS")
oxBS_index <- which(pData(rgSet)$method=="oxBS")
MethylatedBS <- getMeth(MSet.noob)[,BS_index]
UnMethylatedBS <- getUnmeth(MSet.noob)[,BS_index]
MethylatedOxBS <- getMeth(MSet.noob)[,oxBS_index]
UnMethylatedOxBS <- getUnmeth(MSet.noob)[,oxBS_index]
In this example we only have two methods, therefore we can choose between the EM-algorithm and the exact constrained maximum likelihood estimates (using PAVA).
Estimates via the EM-algorithm:
results_emPD2 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS,
L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS,
iterative = TRUE)
The exact constrained MLE (using PAVA):
results_exactPD2 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS,
L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS)
References
samarafk/MLML2R documentation built on Oct. 19, 2019, 6:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
BiocStyle::markdown()
Introduction
In a given CpG site from a single cell we will either have a $C$ or a $T$ after DNA processing conversion methods, with a different interpretation for each of the available methods. This is a binary outcome and we assume a Binomial model and use the maximum likelihood estimation method to obtain the estimates for hydroxymethylation and methylation proportions.
$T$ reads are referred to as converted cytosine and $C$ reads are referred to as unconverted cytosine. Conventionally, $T$ counts are also referred to as unmethylated counts, and $C$ counts as methylated counts. In case of Infinium Methylation arrays, we have intensities representing the methylated (M) and unmethylated (U) channels that are proportional to the number of unconverted and converted cytosines ($C$ and $T$, respectively). The most used summary from these experiments is the proportion $\beta=\frac{M}{M+U}$, commonly referred to as \textit{beta-value}, which reflects the methylation level at a CpG site. Naïvely using the difference between betas from BS and oxBS as an estimate of 5-hmC (hydroxymethylated cytosine), and the difference between betas from BS and TAB as an estimate of 5-mC (methylated cytosine) can many times provide negative proportions and instances where the sum of 5-C (unmodified cytosine), 5-mC and 5-hmC proportions is greater than one due to measurement errors.
\Biocpkg{MLML2R} package allows the user to jointly estimate hydroxymethylation and methylation consistently and efficiently.
The function \Rfunction{MLML} takes as input the data from the different methods and returns the estimated proportion of methylation, hydroxymethylation and unmethylation for a given CpG site. Table 1 presents the arguments of the \Rfunction{MLML} and Table 2 lists the results returned by the function.
The function assumes that the order of the rows and columns in the input matrices are consistent. In addition, all the input matrices must have the same dimension. Usually, rows represent CpG loci and columns are the samples.
Arguments | Description
--------------------|------------------------
\Robject{G.matrix} | Unmethylated channel (Converted cytosines/ T counts) from TAB-conversion (reflecting 5-C + 5-mC).
\Robject{H.matrix} | Methylated channel (Unconverted cytosines/ C counts) from TAB-conversion (reflecting True 5-hmC).
\Robject{L.matrix} | Unmethylated channel (Converted cytosines/ T counts) from oxBS-conversion (reflecting 5-C + 5-hmC).
\Robject{M.matrix} | Methylated channel (Unconverted cytosines/ C counts) from oxBS-conversion (reflecting True 5-mC).
\Robject{T.matrix} | Methylated channel (Unconverted cytosines/ C counts) from standard BS-conversion (reflecting 5-mC+5-hmC).
\Robject{U.matrix} | Unmethylated channel (Converted cytosines/ T counts) from standard BS-conversion (reflecting True 5-C).
: MLML function and random variable notation.
Value | Description ------|------------------ \Robject{mC} | maximum likelihood estimate for the 5-mC proportion \Robject{hmC} | maximum likelihood estimate for the 5-hmC proportion \Robject{C} | maximum likelihood estimate for the 5-mC proportion \Robject{methods} | the conversion methods used to produce the MLE : Results returned from the \Rfunction{MLML} function
Worked examples
Simulated data
\Biocpkg{MLML2R} includes small example datasets for illustration.
We simulated counts from Binomial model with true proportions of methylation, hydroxymethylation and unmethylated being $0.3$, $0.2$, and $0.5$, respectively. For instance, \verb|MethylatedBS_sim| is a matrix of simulated counts from BS corresponding to 100 CpGs and 2 samples. Similarly we simulated matrices with methylated and unmethylated counts for all the three methods: BS, oxBS and TAB. The rows and columns in the input matrices are consistent.
BS and oxBS methods
Load the package:
library(MLML2R)
When only two methods are available, the default option returns the exact constrained maximum likelihood estimates using the the pool-adjacent-violators algorithm (PAVA) [@ayer1955].
results_exactBO1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim)
Maximum likelihood estimate via EM-algorithm approach [@Qu:MLML] is obtained with the option \verb|iterative=TRUE|. In this case, the default (or user specified) \verb|tol| is considered in the iterative method.
results_emBO1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim,iterative=TRUE)
When only two methods are available, we highly recommend the default option \Rcode{iterative=FALSE} since the difference in the estimates obtained via EM and exact constrained is very small, but the former requires more computational effort:
all.equal(results_emBO1$hmC,results_exactBO1$hmC,scale=1)
library(microbenchmark) mbmBO1 = microbenchmark( EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim), EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, iterative=TRUE), times=10) mbmBO1
BS and TAB methods
Using PAVA:
results_exactBT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Using EM-algorithm:
results_emBT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
all.equal(results_emBT1$hmC,results_exactBT1$hmC,scale=1)
mbmBT1 = microbenchmark( EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim), EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim, iterative=TRUE), times=10) mbmBT1
oxBS and TAB methods
Using PAVA:
results_exactOT1 <- MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Using EM-algorithm:
results_emOT1 <- MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
all.equal(results_emOT1$hmC,results_exactOT1$hmC,scale=1)
mbmOT1 = microbenchmark( EXACT = MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim), EM = MLML(L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim, iterative=TRUE), times=10) mbmOT1
BS, oxBS and TAB methods
When data from the three methods are available, the default otion in the \Rfunction{MLML} function returns the constrained maximum likelihood estimates using an approximated solution for Lagrange multipliers method.
results_exactBOT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim)
Maximum likelihood estimate via EM-algorithm approach [@Qu:MLML] is obtained with the option \verb|iterative=TRUE|. In this case, the default (or user specified) \verb|tol| is considered in the iterative method.
results_emBOT1 <- MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim,iterative=TRUE)
We recommend the default option \Rcode{iterative=FALSE} since the difference in the estimates obtained via EM and the approximate exact constrained is very small, but the former requires more computational effort:
all.equal(results_emBOT1$hmC,results_exactBOT1$hmC,scale=1)
mbmBOT1 = microbenchmark( EXACT = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim), EM = MLML(T.matrix = MethylatedBS_sim , U.matrix = UnMethylatedBS_sim, L.matrix = UnMethylatedOxBS_sim, M.matrix = MethylatedOxBS_sim, G.matrix = UnMethylatedTAB_sim, H.matrix = MethylatedTAB_sim, iterative=TRUE), times=10) mbmBOT1
Publicly available data: GSE63179
We will use the dataset from @10.1371/journal.pone.0118202, which consists of eight DNA samples from the same DNA source treated with oxBS-BS and hybridized to the Infinium 450K array.
When data is obtained through Infinium Methylation arrays, we recommend the use of the \Biocpkg{minfi} package [@minfi], a well-established tool for reading, preprocessing and analysing DNA methylation data from these platforms. Although our example relies on \Biocpkg{minfi} and other \Bioconductor{} tools, \Biocpkg{MLML2R} does not depend on any packages. Thus, the user is free to read and preprocess the data using any software of preference and then import the intensities (or $T$ and $C$ counts) for the methylated and unmethylated channel (or converted and uncoverted cytosines) into \R{} in matrix format.
To start this example we will need the following packages:
library(minfi) library(GEOquery)
It is usually best practice to start the analysis from the raw data, which in the case of the 450K array is a \verb|.IDAT| file.
The raw files are deposited in GEO and can be downloaded by using the \Rfunction{getGEOSuppFiles}. There are two files for each replicate, since the 450k array is a two-color array. The \verb|.IDAT| files are downloaded in compressed format and need to be uncompressed before they are read by the \Rfunction{read.metharray.exp} function.
getGEOSuppFiles("GSE63179") untar("GSE63179/GSE63179_RAW.tar", exdir = "GSE63179/idat") list.files("GSE63179/idat", pattern = "idat") files <- list.files("GSE63179/idat", pattern = "idat.gz$", full = TRUE) sapply(files, gunzip, overwrite = TRUE)
The \verb|.IDAT| files can now be read:
rgSet <- read.metharray.exp("GSE63179/idat")
To access phenotype data we use the \Rfunction{pData} function. The phenotype data is not yet available from the \Robject{rgSet}.
pData(rgSet)
In this example the phenotype is not really relevant, since we have only one sample: male, 25 years old. What we do need is the information about the conversion method used in each replicate: BS or oxBS. We will access this information automatically from GEO:
geoMat <- getGEO("GSE63179") pD.all <- pData(geoMat[[1]]) pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1.1", "characteristics_ch1.2","characteristics_ch1.3")] pD
This phenotype data needs to be merged into the methylation data. The following commands guarantee we have the same replicate identifier in both datasets before merging.
sampleNames(rgSet) <- sapply(sampleNames(rgSet),function(x) strsplit(x,"_")[[1]][1]) rownames(pD) <- pD$geo_accession pD <- pD[sampleNames(rgSet),] pData(rgSet) <- as(pD,"DataFrame") rgSet
The \Robject{rgSet} object is a class called \Rclass{RGChannelSet} used for two color data (green and a red channel). The input in the \Rfunction{MLML} funcion is \Rclass{MethylSet}, which contains the methylated and unmethylated signals. The most basic way to construct a \Rclass{MethylSet} is using the function \Rfunction{preprocessRaw}. Here we chose the function \Rfunction{preprocessNoob} [@noob] for background correction and construction of the \Rclass{MethylSet}.
MSet.noob<- preprocessNoob(rgSet)
After the preprocessed steps we can use \Rfunction{MLML} from the \Biocpkg{MLML2R} package.
The BS replicates are in columns 1, 3, 5, and 6. The remaining columns are from the oxBS treated replicates.
MethylatedBS <- getMeth(MSet.noob)[,c(1,3,5,6)] UnMethylatedBS <- getUnmeth(MSet.noob)[,c(1,3,5,6)] MethylatedOxBS <- getMeth(MSet.noob)[,c(7,8,2,4)] UnMethylatedOxBS <- getUnmeth(MSet.noob)[,c(7,8,2,4)]
In this example we only have two methods, therefore we can choose between the EM-algorithm and the exact constrained maximum likelihood estimates (using PAVA).
Estimates via the EM-algorithm:
results_emPD1 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS, L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS, iterative = TRUE)
The exact constrained MLE (using PAVA):
results_exactPD1 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS, L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS)
Publicly available data: GSE73895
We will use the dataset from @pmid27886174, which consists of 30 DNA samples from glioblastoma tumors treated with oxBS-BS and hybridized to the Infinium 450K array.
The raw files are deposited in GEO and can be downloaded and read into \R{} by doing:
getGEOSuppFiles("GSE73895") untar("GSE73895/GSE73895_RAW.tar", exdir = "GSE73895/idat") idatFiles <- list.files("GSE73895/idat", pattern = "idat.gz$", full = TRUE) sapply(idatFiles, gunzip, overwrite = TRUE) rgSet <- read.metharray.exp("GSE73895/idat")
We need to identify the samples from different methods: BS-conversion, oxBS-conversion. We obtain this information from GEO:
geoMat <- getGEO("GSE73895") pD.all <- pData(geoMat[[1]]) pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1", "characteristics_ch1.2","characteristics_ch1.3")] head(pD)
Keeping only some of the variables from phenotype data:
names(pD)[c(1,3,4,5)] <- c("method","gender","survival_months","age_years") pD$gender <- sub("^gender: ", "", pD$gender) pD$age_years <- as.numeric(sub("^subject age: ", "", pD$age_years)) pD$survival_months <- as.numeric(sapply(pD$survival_months, function(x) strsplit(as.character(x),":")[[1]][2])) pD$method <- sapply(pD$method, function(x) strsplit(as.character(x),"_")[[1]][3])
We now need to merge this pheno data into the methylation data. The following are commands to make sure we have the same row identifier in both datasets before merging.
sampleNames(rgSet) <- sapply(sampleNames(rgSet),function(x) strsplit(x,"_")[[1]][1]) rownames(pD) <- pD$geo_accession pD <- pD[sampleNames(rgSet),] pData(rgSet) <- as(pD,"DataFrame") rgSet
Preprocessing and preparing input matrices for \Rfunction{MLML}:
MSet.noob<- preprocessNoob(rgSet) BS_index <- which(pData(rgSet)$method=="BS") oxBS_index <- which(pData(rgSet)$method=="oxBS") MethylatedBS <- getMeth(MSet.noob)[,BS_index] UnMethylatedBS <- getUnmeth(MSet.noob)[,BS_index] MethylatedOxBS <- getMeth(MSet.noob)[,oxBS_index] UnMethylatedOxBS <- getUnmeth(MSet.noob)[,oxBS_index]
In this example we only have two methods, therefore we can choose between the EM-algorithm and the exact constrained maximum likelihood estimates (using PAVA).
Estimates via the EM-algorithm:
results_emPD2 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS, L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS, iterative = TRUE)
The exact constrained MLE (using PAVA):
results_exactPD2 <- MLML(T.matrix = MethylatedBS , U.matrix = UnMethylatedBS, L.matrix = UnMethylatedOxBS, M.matrix = MethylatedOxBS)
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.