Statistical model for Joint Analysis of Methylation and Genotype Data

Single nucleotide polymorphisms (SNPs) can create and destroy CpGs. As methylation occurs mostly at CpGs, such CpG-SNPs can directly affect methylation measurements.

Recall that enrichment-based methylation methods measure total methylation in a vicinity of a CpG. By creating or destroying a CpG, CpG-SNPs introduce a variation in the total methylation in a vicinity of the CpG which can greatly reduce our power to detect case-control differences.

RaMWAS can account for a possible effect of CpG-SNPs by testing for joint significance of β₁ and β₂ the following model:

μ_i = β₀ + outcome * β₁ + outcome * SNP_i * β₂ + SNP_i * β₃ + γ * cvrt + ϵ

where

• μ_i – methylation measurement for i-th CpG.
• outcome – phenotype of interest.
• SNP_i – the SNP values (0/1/2 or dosages for imputed genotype) at the i-th CpG.
• cvrt – covariates and the principal components.
• ϵ – noise.

Input data

For CpG-SNPs analysis RaMWAS requires the usual input (see steps 4 and 5) with an additional SNP matrix.

The SNP data must have the same dimensions as the CpG score matrix, i.e. it must be available for the same set of samples and the same set of locations. Data preparation may include finding the closest SNP for every CpG and exclusion of CpGs without any SNPs in vicinity.

library(ramwas)

# work in a temporary directory
dr = paste0(tempdir(), "/simulated_matrix_data")
dir.create(dr, showWarnings = FALSE)
cat(dr,"\n")

## /tmp/RtmpE42VVP/simulated_matrix_data

nsamples = 200
nvariables = 100000

covariates = data.frame(
    sample = paste0("Sample_",seq_len(nsamples)),
    sex = seq_len(nsamples) %% 2,
    age = runif(nsamples, min = 20, max = 80),
    batch = paste0("batch",(seq_len(nsamples) %% 3))
)
pander(head(covariates))

temp = cumsum(sample(20e7 / nvariables, nvariables, replace = TRUE) + 0)
chr      = as.integer(temp %/% 1e7) + 1L
position = as.integer(temp %% 1e7)

locmat = cbind(chr = chr, position = position)
chrnames = paste0("chr", 1:10)
pander(head(locmat))

fmloc = fm.create.from.matrix(
            filenamebase = paste0(dr, "/CpG_locations"),
            mat = locmat)
close(fmloc)
writeLines(con = paste0(dr, "/CpG_chromosome_names.txt"), text = chrnames)

fmm = fm.create(paste0(dr,"/Coverage"), nrow = nsamples, ncol = nvariables)
fms = fm.create(paste0(dr,"/SNPs"), nrow = nsamples, ncol = nvariables,
                size = 1, type = "integer")

# Row names of the matrices are set to sample names
rownames(fmm) = as.character(covariates$sample)
rownames(fms) = as.character(covariates$sample)

# The matrices are filled, 2000 variables at a time
byrows = 2000
for( i in seq_len(nvariables/byrows) ){ # i=1
    ind = (1:byrows) + byrows*(i-1)

    snps = rbinom(n = byrows * nsamples, size = 2, prob = 0.2)
    dim(snps) = c(nsamples, byrows)
    fms[,ind] = snps

    slice = double(nsamples*byrows)
    dim(slice) = c(nsamples, byrows)
    slice[,  1:225] = slice[,  1:225] + covariates$sex / 50 / sd(covariates$sex)
    slice[,101:116] = slice[,101:116] + covariates$age / 16 / sd(covariates$age)
    slice = slice +
                ((as.integer(factor(covariates$batch))+i) %% 3) / 200 +
                snps / 1.5 +
                runif(nsamples*byrows) / 2
    fmm[,ind] = slice
}
close(fms)
close(fmm)

SNP-CpG analysis

Let us test for association between CpG scores and and the sex covariate (modeloutcome parameter) correcting for batch effects (modelcovariates parameter). Save top 20 results (toppvthreshold parameter) in a text file.

param = ramwasParameters(
    dircoveragenorm = dr,
    covariates = covariates,
    modelcovariates = "batch",
    modeloutcome = "sex",
    toppvthreshold = 20,
    fileSNPs = "SNPs"
)

The CpG-SNP analysis:

ramwasSNPs(param)

The QQ-plot shows better enrichment with significant p-values.

For comparison, we also perform the usual MWAS for these CpGs without regard for SNPs.

ramwas5MWAS(param)

The QQ-plot shows much weaker signal for the standard MWAS.

The top finding are saved in the text files Top_tests.txt for both analyses:

# Get the directory with testing results
toptbl = read.table(
                paste0(pfull$dirSNPs, "/Top_tests.txt"),
                header = TRUE,
                sep = "\t")
pander(head(toptbl,10))

Note that CpG-SNP analysis tests for joint significance of β₁ and β₂ and thus uses F-test, while regular MWAS uses t-test.

pfull = parameterPreprocess(param)
toptbl = read.table(
                paste0(pfull$dirmwas, "/Top_tests.txt"),
                header = TRUE,
                sep = "\t")
pander(head(toptbl,10))