Hello,
I would like to test for differential methylation variability using missMethyl, but with a block design to take sample pairs/quartets into account.
If I wanted to test for differential (mean) methylation with my experimental design, I would use a code like this for example:
design <- model.matrix(~0+A1+A2+B1+B2) blocks <- id dupcor <- duplicateCorrelation(object=mval, design=design, block=blocks) fit <- lmFit(mval, design, blocks, correlation=dupcor$consensus) contrast.matrix <- makeContrasts( contrastA=A2-A1, contrastB=B2-B1, contrast1=B1-A1, contrast2=B2-A2, levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) topTable(fit2, coef="contrastA") topTable(fit2, coef="contrastB") topTable(fit2, coef="contrast1") topTable(fit2, coef="contrast2")
This works fine.
However, when I try to "translate" this into code for missMethyl, the following happens:
> fitvar <- varFit(mval, design, blocks, correlation=dupcor$consensus) Error in varFit(mval, design, blocks, correlation = dupcor$consensus) : unused argument (correlation = dupcor$consensus)
I used blocks because I have samples from the same (or a subset of the same) individuals across all conditions and would like to take this into account. As the block design didn't work, I also tried to simply add the ids (which are factors) as covariates to the design, but then I got the following warning and no results (all NA):
> design <- model.matrix(~0+A1+A2+B1+B2+id) > [...] > fitvar <- varFit(mval, design=design) Warning messages: 1: Partial NA coefficients for 466566 probe(s)
Does anybody know how I could properly test for differential variability taking the information about the individuals into account?
Best,
Simone
> sessionInfo() R version 3.2.1 (2015-06-18) Platform: x86_64-unknown-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux Server 7.3 (Maipo) locale: [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8 [6] LC_MESSAGES=en_GB.UTF-8 LC_PAPER=en_GB.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] limma_3.26.9 missMethyl_1.4.0 RSQLite_1.0.0 DBI_0.5-1