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

Description Usage Arguments Details Value Examples

View source: R/lmmTestAllRegions.R

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
   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.