Running EWAS

Four regression models

Meffil has four different EWAS models implemented:

1. none=no covariates
2. all=user covariates
3. sva=surrogate variables in addition to user covariates
4. isva=independent surrogate variables in addition to user covariates

See here for a link to the isva paper.

Running EWAS

load(beta_file) # methylation betas: CpGs in rows, samples in cols
                # assume that it loads a matrix called 'norm.beta'
covs <- read.table(covs_file, he=T, stringsAsFactors=FALSE) #covariate file: samples in rows, covs in cols
                # assume that the rownames of 'covs' correspond to the column names of 'norm.beta'
phen <- read.table(phen_file, he=T, stringsAsFactors=FALSE) #phenotype file: samples in rows, phenotypes in cols 
                # assume that the rownames of 'phen' correspond to the column names of 'norm.beta'
                # assume that the first column of 'phen' is the phenotype of interest for the EWAS

phen<-na.omit(phen)

norm.beta <- norm.beta[, colnames(norm.beta) %in% rownames(phen)]
phen <- phen[match(colnames(norm.beta), rownames(phen)), , drop=FALSE]
stopifnot(identical(rownames(phen), colnames(norm.beta)))
    
covs <- covs[match(colnames(norm.beta), rownames(covs)), , drop=FALSE]
stopifnot(identical(rownames(covs), colnames(norm.beta)))

most.variable <- min(nrow(norm.beta), 20000) ## number most variable probes for ISVA
    
ewas.ret <- meffil.ewas(norm.beta, variable=phen[,1], covariates=covs, 
                        winsorize.pct = NA, ## do not handle outliers by winsorization
                        most.variable = most.variable) 

ewas.parameters <- meffil.ewas.parameters(sig.threshold=1e-7,  ## EWAS p-value threshold
                                          max.plots=100, ## plot at most 100 CpG sites
                                          qq.inflation.method="regression",  ## measure inflation using regression
                                          model="sea") ## select default EWAS model

ewas.summary<-meffil.ewas.summary(ewas.ret,norm.beta,parameters=ewas.parameters)                              

meffil.ewas.report(ewas.summary, output.file=paste(report_file,".ewas.report.html",sep=""))