Hello all,

I am trying to perform differential methylation analysis with epic array data. I have loaded my idat files into R environment. I filtered and normalized the data with preprocessFunnorm. The next steps are given below. There is something off as the correlation is low and the summary data seems wrong. At the end I also put the beta value matrix. R01C01 and R05C01 are replicates of SD-1. R02C01 and R06C01 are replicates of SD-1R and it goes like that. I don't know if duplicateCorrelation or avereps are better in this case. Also for avereps I tried to use it but it gives me some sort of error.

I would appreciate if you could help me.

> factor <- factor(targets$Condition)
> factor
[1] SD_1    SD_1R  SUP_B15   SUP_B15R SD_1   
[6]  SD_1R  SUP_B15   SUP_B15R
Levels: SUP_B15 SUP_B15R SD_1    SD_1R
> design <- model.matrix(~0+factor)
> colnames(design) <- c("SD_1", "SD_1R", "SUP_B15", "SUP_B15R")
> design
  SD_1 SD_1R SUP_B15 SUP_B15R
1      1        0       0         0
2      0        1       0         0
3      0        0       1         0
4      0        0       0         1
5      1        0       0         0
6      0        1       0         0
7      0        0       1         0
8      0        0       0         1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$factor
[1] "contr.treatment"
> biolrep <- c(1, 1, 2, 2, 3, 3, 4, 4)
> corfit <- duplicateCorrelation(beta, ndups = 1, block = biolrep)
> corfit$consensus.correlation
[1] 0.2776849
> contMatrix <- makeContrasts(
     SD_1R_vs_SD_1 = SD_1R - SD_1,
     SUP_B15R_vs_SUP_B15 = SUP_B15R - SUP_B15,
     SUP_B15_vs_SD_1 = SUP_B15 - SD_1,
     SUP_B15R_vs_SD_1R = SUP_B15R - SD_1R,
     levels = design
 )
> contMatrix
          Contrasts
Levels     SD_1R_vs_SD_1 SUP_B15R_vs_SUP_B15 SUP_B15_vs_SD_1 SUP_B15R_vs_SD_1R
  SD_1                -1                   0              -1                 0
  SD_1R                1                   0               0                -1
  SUP_B15              0                  -1               1                 0
  SUP_B15R             0                   1               0                 1
> fit <- lmFit(beta_val, design, block = biolrep, cor=corfit$consensus.correlation)
> fit2 <- contrasts.fit(fit, contMatrix)
> fit2 <- eBayes(fit2)
> summary(decideTests(fit2))
       SD_1R_vs_SD_1 SUP_B15R_vs_SUP_B15
Down                  412                  414
NotSig             119961               119936
Up                    377                  400
       SUP_B15_vs_SD_1 SUP_B15R_vs_SD_1R
Down                 368                   302
NotSig            120032                120162
Up                   350                   286

> head(beta)
                        R01C01              R02C01              R03C01              R04C01
cg14817997           0.1705210           0.1669450           0.1667081           0.1643278
cg01803908           0.1714407           0.1720723           0.1680338           0.1659238
cg15174812           0.1715732           0.1722038           0.1681296           0.1660180
cg05001044           0.1714453           0.1709194           0.1672138           0.1652023
cg17501828           0.1700388           0.1672769           0.1664959           0.1643007
cg23917638           0.1713424           0.1719742           0.1679110           0.1657868
                        R05C01              R06C01              R07C01              R08C01
cg14817997           0.1676968           0.1614067           0.1732130           0.1682127
cg01803908           0.1716572           0.1651171           0.1748817           0.1704308
cg15174812           0.1716243           0.1651899           0.1750537           0.1705295
cg05001044           0.1712679           0.1641264           0.1738601           0.1692805
cg17501828           0.1663494           0.1618090           0.1727334           0.1677366
cg23917638           0.1715577           0.1650289           0.1747092           0.1703147

You shouldn't be using beta values for a linear model that assumes a unimodal, untruncated distribution (beta values are strictly on (0-1) and tend to be grouped at the extremes of that range). Instead you should be using M-values, which are unimodal and on (-Inf, Inf).