Fit mixed model to test association between a continuous phenotype and methylation values in a list of genomic regions
data frame or matrix of beta values for all genomic regions, with row names = CpG IDs and column names = sample IDs. This is often the genome-wide array data.
a list of genomic regions; each item is a vector of CpG IDs
within a genomic region. The co-methylated regions can be obtained by
function CoMethAllRegions.
a data frame with phenotype and covariates, with variable
Sample indicating sample IDs.
character string of the main effect (a continuous phenotype) to be tested for association with methylation values in each region
character vector for names of the covariate variables
type of mixed model; can be randCoef for random
coefficient mixed model or simple for simple linear mixed model.
Human genome of reference: hg19 or hg38
Type of array: "450k" or "EPIC"
Whether strand can be ignored, default is TRUE
output .csv file with the results for the mixed model analysis
log file for mixed models analysis messages
Number of computing cores to be used when executing code in parallel. Defaults to 1 (serial computing).
Dots for additional arguments passed to the cluster constructor.
See CreateParallelWorkers for more information.
If outFile is NULL, this function returns a data frame
as described below. If outFile is specified, this function writes a
data frame in .csv format with the following information to the disk:
location of the genomic region (chrom, start, end), number of CpGs
(nCpGs), Estimate, Standard error (StdErr) of the test
statistic, p-value and False Discovery Rate (FDR) for association between
methylation values in each genomic region with phenotype (pValue).
If outLogFile is specified, this function also writes a .txt file
that includes messages for mixed model fitting to the disk.
This function implements a mixed model to test association between methylation M values in a genomic region with a continuous phenotype. In our simulation studies, we found both models shown below are conservative, so p-values are estimated from normal distributions instead of Student's t distributions.
When modelType = "randCoef", the model is:
M ~ contPheno_char + covariates_char + (1|Sample) + (contPheno_char|CpG).
The last term specifies random intercept and slope for each CpG. When
modelType = "simple", the model is
M ~ contPheno_char + covariates_char + (1|Sample).
For the results of mixed models, note that if the mixed model failed to converge, p-value for mixed model is set to 1. Also, if the mixed model is singular, at least one of the estimated variance components for intercepts or slopes random effects is 0, because there isn't enough variability in the data to estimate the random effects. In this case, the mixed model reduces to a fixed effects model. The p-values for these regions are still valid.
data(betasChr22_df)
data(pheno_df)
CpGisland_ls <- readRDS(
system.file(
"extdata",
"CpGislandsChr22_ex.rds",
package = 'coMethDMR',
mustWork = TRUE
)
)
coMeth_ls <- CoMethAllRegions(
dnam = betasChr22_df,
betaToM = TRUE,
CpGs_ls = CpGisland_ls,
arrayType = "450k",
rDropThresh_num = 0.4,
returnAllCpGs = FALSE
)
#> snapshotDate(): 2021-05-18
#> see ?sesameData and browseVignettes('sesameData') for documentation
results_df <- lmmTestAllRegions(
betas = betasChr22_df,
region_ls = coMeth_ls,
pheno_df = pheno_df,
contPheno_char = "stage",
covariates_char = "age.brain",
modelType = "randCoef",
arrayType = "450k",
ignoreStrand = TRUE,
# generates a log file in the current directory
# outLogFile = paste0("lmmLog_", Sys.Date(), ".txt")
)
#> snapshotDate(): 2021-05-18
#> see ?sesameData and browseVignettes('sesameData') for documentation
#> Analyzing region chr22:18268062-18268249.
#> Analyzing region chr22:18324579-18324769.
#> Analyzing region chr22:18531243-18531447.
#> For future calls to this function, perhaps specify a log file.
#> Set the file name of the log file with the outLogFile argument.