In perishky/meffil: Efficient algorithms for DNA methylation

library(knitr)
library(Cairo)
opts_chunk$set(warning=FALSE, fig.width=6, fig.height=6, dev="CairoPNG", stop=TRUE)

library(minfi, quietly=TRUE)
library(FlowSorted.CordBlood.450k, quietly=TRUE)
#source("http://bioconductor.org/biocLite.R")
#biocLite("FlowSorted.CordBlood.450k")

library(GEOquery)

EWAS using different cell type references

Download example data set

path <- download.450k.lead.dataset()

Normalize dataset

Create samplesheet

library(meffil)
options(mc.cores=5)
samplesheet <- meffil.create.samplesheet(path)

samples <- read.csv(file.path(path, "samples.csv"), check.names=F, row.names=1)
samplesheet <- data.frame(samplesheet,
                          samples[match(samplesheet$Sample_Name, rownames(samples)),],
                          stringsAsFactors=F, check.names=F)

samplesheet <- samplesheet[which(samplesheet[["sample type"]] == "HM450K"),]

Parameters.

qc.file <- "cord/qc-report.html"
author <- "Sen, et al."
study <- "Cord blood DNA methylation and lead exposure (GSE69633)"
norm.file <- "cord/normalization-report.html"
cell.type.reference <- "gervin and lyle cord blood"

Generate QC objects for each sample and QC report.

qc.objects <- meffil.qc(samplesheet, cell.type.reference=cell.type.reference, verbose=T)

qc.summary <- meffil.qc.summary(qc.objects, verbose=T)
meffil.qc.report(qc.summary,
                 output.file=qc.file,
                 author=author,
                 study=study)

Remove any low quality samples.

if (nrow(qc.summary$bad.samples) > 0)
    qc.objects <- meffil.remove.samples(qc.objects, qc.summary$bad.samples$sample.name)

samplesheet <- samplesheet[match(names(qc.objects), rownames(samplesheet)),]

Check how many principal components to include.

print(meffil.plot.pc.fit(qc.objects, n.cross=3)$plot)

Ten seems about right.

number.pcs <- 10

Normalize dataset and generate normalization report.

norm.objects <- meffil.normalize.quantiles(qc.objects, number.pcs=number.pcs, verbose=T)
norm.dataset <- meffil.normalize.samples(norm.objects,
                                        just.beta=F, 
                                        cpglist.remove=qc.summary$bad.cpgs$name,
                                        verbose=T)

beta <- meffil.get.beta(norm.dataset$M, norm.dataset$U)
pcs <- meffil.methylation.pcs(beta, sites=meffil.get.autosomal.sites("450k"), verbose=T)

parameters <- meffil.normalization.parameters(norm.objects)
norm.summary <- meffil.normalization.summary(norm.objects=norm.objects,
                                             pcs=pcs,
                                             parameters=parameters, verbose=T)
meffil.normalization.report(norm.summary,
                            output.file=norm.file,
                            author=author,
                            study=study)

Estimate cell counts using all available references

counts <- sapply(meffil.list.cell.type.references(), function(reference) {
    cat(date(), reference, "\n")
    do.call(rbind, mclapply(qc.objects, function(qc.object) {
        meffil.estimate.cell.counts(qc.object, cell.type.reference=reference)$counts
    }))
}, simplify=F)

Compare counts between references

r <- sapply(names(counts), function(n1)
       sapply(names(counts), function(n2) {
           sapply(intersect(colnames(counts[[n1]]), colnames(counts[[n2]])),
                  function(type) cor(counts[[n1]][,type], counts[[n2]][,type]))
       }))

types <- unique(unlist(lapply(counts, function(counts) colnames(counts))))
for (type in types) {
    cat("Cell type", type, "-----------------\n")
    for (n1 in colnames(r))
        for (n2 in rownames(r)) 
            if (n1 < n2 && type %in% names(r[n1,n2][[1]]))
                cat(r[n1,n2][[1]][[type]], n1, "--", n2, "\n")    
}

Variables of interest

We will test associations with gender, gestational age, birth weight and blood lead levels. The following list provides the names of potential confounders to be included in EWAS regression models.

covariate.names <- list("gender"=c("socioeconomic score","gestational age","smoke ever","birth weight"),
                        "gestational age"=c("socioeconomic score","gender","smoke ever","birth weight"),
                        "birth weight"=c("socioeconomic score","gestational age","smoke ever","gender"),
                        "pbconc (ng/dl)"=c("socioeconomic score","gestational age","smoke ever","birth weight","gender"))

For the gender EWAS, we will test only CpG sites on autosomes.

features <- meffil.featureset("450k")
autosomal <- features[which(!(features$chromosome %in% c("chrX","chrY"))),]

Test associations

Test associations for each variable of interest (4) and each cell type reference.

results <- sapply(names(covariate.names), function(variable.name) {
    variable <- samplesheet[[variable.name]]
    sapply(names(counts), function(reference) {
        covariates <- data.frame(counts[[reference]],
                                 samplesheet[,covariate.names[[variable.name]]],
                                 stringsAsFactors=F)
        cat(date(), "ewas", variable.name, reference, "\n")
        if (variable.name == "gender")
            beta <- beta[which(rownames(beta) %in% autosomal$name),]
        meffil.ewas(beta, variable=variable, covariates=covariates,
                    isva=F, sva=F, smartsva=T,
                    most.variable=50000,verbose=T)
    }, simplify=F)
}, simplify=F)

Previously identified associations

For generating the EWAS reports, previously identified associations are provided for comparison with previous analyses. Only the lead ("pbconc") list of associations is from the dataset used here.

published.hits <- list("gender"=c("cg03691818", "cg26921482", "cg17743279", "cg07852945", "cg26355737", "cg25568337", "cg05100634", "cg03608000", "cg02325951", "cg17612569", "cg04874129", "cg08906898", "cg04946709", "cg02989351", "cg12204423", "cg25304146", "cg22345911", "cg01906879", "cg06152526", "cg04190002", "cg06644124", "cg07628841", "cg23001456", "cg26213873", "cg25438440", "cg07816873", "cg24016844", "cg11841231", "cg13323902", "cg12900929"), ## PMID: 26553366 (top 30)
                       "gestational age"=c("cg08943494", "cg11932158", "cg16725984", "cg20334115", "cg18623216", "cg07835443", "cg00220721", "cg04685228", "cg16103712", "cg27518892", "cg22117805", "cg21926626", "cg04347477", "cg08817867", "cg27448161", "cg01154283", "cg02001279", "cg02430430", "cg06870470", "cg13924996", "cg07136133", "cg00481600", "cg13675859", "cg11934771", "cg19744173", "cg12713583", "cg11360522", "cg16834726", "cg18608055", "cg05283597"), ## PMID: 25869828 (top 30)
                       "birth weight"=c("cg20076442", "cg25953130", "cg04521626", "cg14097568", "cg17133774", "cg00654448", "cg00442282", "cg13696490", "cg12044213", "cg08817867", "cg00382138", "cg06870470", "cg25557739", "cg24324628", "cg15783941", "cg14597739", "cg22962123", "cg05851442", "cg23387597", "cg24973755", "cg16219283", "cg25799241", "cg06658067"), ## PMID: 25869828 (bonferroni p < 0.05)
                       "pbconc (ng/dl)"=c("cg01105385", "cg07208333", "cg08945395", "cg20439288", "cg20474370"))## PMID: 26046694, all sites in chr5:67584213-67584451

EWAS reports

Generate a report for each EWAS. Significance thresholds are selected so that the top 100 associations are plotted in the report. The reports can be found in the "cord" directory under the name of the variable of interest (e.g. "cord/gender").

for (variable.name in names(results)) {
    for (reference in names(results[[variable.name]])) {
        cat(date(), "EWAS report", variable.name, reference, "\n")

        ewas.ret <- results[[variable.name]][[reference]]

        ewas.parameters <- meffil.ewas.parameters(max.plots=10)

        ewas.summary <- meffil.ewas.summary(ewas.ret,
                                            beta,
                                            selected.cpg.sites=published.hits[[variable.name]],
                                            parameters=ewas.parameters, verbose=T)
        meffil.ewas.report(ewas.summary,
                           output.file=file.path("cord", sub(" \\(.*", "", variable.name),
                               reference, "report.html"),
                           author="Me",
                           study=paste(variable.name,
                               " in cord blood DNA methylation (GEO:GSE69633)",
                               " adjusted with the '", reference, "' cell type reference.", sep=""))
    }
}

Association overlaps between different analyses

Show size of overlap between analyses using different cell type references. For each, we use the results obtained using confounders obtained using "smartsva" (all potential confounders listed above plus surrogate variables).

compare.top <- function(results, top=100, method) {
    sites <- lapply(results, function(ewas) {
        idx <- order(ewas$p.value[,method], decreasing=F)[1:top]
        rownames(ewas$p.value)[idx]
    })
    sapply(sites, function(x) sapply(sites, function(y) length(intersect(x,y))))
}
comparison.table <- function(results, method) {
    comparisons <- lapply(results, compare.top, method=method)
    for (name in names(results)) {
        comparisons[[name]] <- as.data.frame(comparisons[[name]])
        comparisons[[name]]$reference <- rownames(comparisons[[name]])
        comparisons[[name]]$variable <- name
    }
    comparisons <- do.call(rbind, comparisons)
}
knitr:::kable(comparison.table(results, "smartsva"), row.names=F)

Significance ranks of previously identified associations

Below we show the p-value ranks in our analyses of each previously identified association.

compare.p.value.ranks <- function(results, method) {
    ranks <- sapply(names(results), function(variable.name) {
        do.call(rbind, lapply(results[[variable.name]], function(ewas) {
            idx <- match(published.hits[[variable.name]], rownames(ewas$p.value))
            r <- rank(ewas$p.value[,method])[idx]
            p <- wilcox.test(ewas$p.value[idx,method],
                             ewas$p.value[-idx,method],
                             alternative="less")$p.value
            c(quantile(r), p=p)
        }))
    }, simplify=F)
    for (name in names(results)) {
        ranks[[name]] <- as.data.frame(ranks[[name]])
        ranks[[name]]$reference <- rownames(ranks[[name]])
        ranks[[name]]$variable <- name
    }
    ranks <- do.call(rbind, ranks)
}
knitr:::kable(compare.p.value.ranks(results, "smartsva"), row.names=F)