Test for differential variability with missMethyl (DiffVar) using block design
1
0
Entering edit mode
Simone ▴ 190
@simone-5854
Last seen 6.4 years ago

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   
missmethyl diffvar design matrix limma block • 1.5k views
ADD COMMENT
0
Entering edit mode
lhuang7 ▴ 50
@lhuang7-7824
Last seen 4.5 years ago
United States

I got the similar warning message when calling varFit under a contrast design matrix setting (multiple comparison groups).

> fitvar <- varFit(mval, design = design)
Warning messages:
1: Partial NA coefficients for 466566 probe(s)

The issue turns out to be a typo within the function varFit.default() in the source code DiffVar.R:

if(is.null(coef)) 
   coef <- c(1, ncol(design))
z <- getLeveneResiduals(data, design = design[,coef], type = type)

As you can see, the coefficients only take the first and the last comparison group. It should be coef <- 1:ncol(design).

I hope the authors can fix it in the next release but meanwhile you could just modify this line of code by yourself and re-install the package from the fixed source file. HTH.

 

 

 

ADD COMMENT

Login before adding your answer.

Traffic: 879 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6