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 sva paper and here for 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)))

set.seed(123456) #to make siva and sva reproducible you need to set a seed.    
ewas.ret <- meffil.ewas(norm.beta, variable=phen[,1], covariates=covs, isva=F) 
save(ewas.ret, file="myewas.Robj")

# This function sets your default model and your significance threshold. It will show statistics for all four models for all CpGs that pass significance threshold with the default model.

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="median",  ## measure inflation using median
                                          model="sva") ## 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=""))