lmmTestAllRegions: Fit mixed model to test association between a continuous...

In lissettegomez/coMethDMR: Accurate identification of co-methylated and differentially methylated regions in epigenome-wide association studies

Description

Fit mixed model to test association between a continuous phenotype and methylation values in a list of genomic regions

Usage

lmmTestAllRegions(
  betas,
  region_ls,
  pheno_df,
  contPheno_char,
  covariates_char,
  modelType = c("randCoef", "simple"),
  genome = c("hg19", "hg38"),
  arrayType = c("450k", "EPIC"),
  outFile = NULL,
  outLogFile = NULL,
  nCores_int = 1L,
  ...
)

Arguments

betas	data frame or matrix of beta values for all genomic regions, with row names = CpG IDs, column names = sample IDs. This is often the genome-wide array data.
region_ls	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.
pheno_df	a data frame with phenotype and covariates, with variable Sample indicating sample IDs.
contPheno_char	character string of the main effect (a continuous phenotype) to be tested for association with methylation values in each region
covariates_char	character vector for names of the covariate variables
modelType	type of mixed model, can be randCoef for random coefficient mixed model, or simple for simple linear mixed model.
genome	Human genome of reference hg19 or hg38
arrayType	Type of array, can be "450k" or "EPIC"
outFile	output .csv file with the results for the mixed model analysis
outLogFile	log file for mixed models analysis messages
nCores_int	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.

Details

This function implements a mixed model to test association between methylation values in a genomic region with a continuous phenotype.

When randCoef is selected, the model is

methylation M value ~ contPheno_char + covariates_char + (1|Sample) + (contPheno_char|CpG). The last term specifies both random intercept and slope for each CpG.

When simple is selected, the model is

methylation M value ~ contPheno_char + covariates_char + (1|Sample)

In our simulation studies, we found both models are conservative, so p-values are estimated from normal distributions instead of t-distributions.

For the results of mixed models, note that

(1) When mixed model failed to converge, p-value for mixed model is set to 1.

(2) When 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 variabilities in data to estimate the random effects. In this case, mixed model reduces to a fixed effects model. The p-values for these regions are still valid.

Value

(1) output file: a .csv file with 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).

(2) log file: a .txt file that includes messages for mixed model fitting

Examples

   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
   )


   results <- lmmTestAllRegions(
     betas = betasChr22_df,
     region_ls = coMeth_ls,
     pheno_df,
     contPheno_char = "stage",
     covariates_char = "age.brain",
     modelType = "randCoef",
     arrayType = "450k"
     # generates a log file in the current directory
     # outLogFile = paste0("lmmLog_", Sys.Date(), ".txt")
   )

lissettegomez/coMethDMR documentation built on April 25, 2021, 1:10 p.m.