Struggle with the replicates of epic array data
1
0
Entering edit mode
@ceaa2928
Last seen 4 days ago
Türkiye

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
limma epicarray IlluminaHumanMethylationEPICanno.ilm10b2.hg19 • 272 views
ADD COMMENT
1
Entering edit mode
@james-w-macdonald-5106
Last seen 4 days ago
United States

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).

0
Entering edit mode

Thanks you so much for the instructions. I also tried M values as you suggested but sadly the result was the same. Apparently, the problem was the normalization method. I was using preprocessFunnorm. I changed it to preprocessQuantile and the problem was solved. I don't know why funnorm did not work in my case

ADD REPLY

Login before adding your answer.

Traffic: 701 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