Hello Bioconductor community,
I am currently analyzing methylation array data using the EPIC array. The study is case-control, and my goal is to identify differentially methylated positions (DMPs) and differentially methylated regions (DMRs).
However, I'm unsure whether I need to set up a contrast for the case-control comparison while keeping the intercept at 0,
Or as the second group of codes if I should proceed using the limma package without specifying a contrast which is going to be the reference group.
My goal is to get DMPs with adjusting for age and gender.
The codes are below
S_Group <- factor(sample_sheet$HTN) # Case and Control
age <- sample_sheet$AGE # Continuous variable
gender <- factor(sample_sheet$GENDER) # Male and Female
# 1A. Adjusting for Age and gender
# Design Matrix
design <- model.matrix(~0 + S_Group + age + gender , data=sample_sheet)
colnames(design)
# Naming the columns
colnames(design) <- gsub("genderM", "sexMale", colnames(design))
# Fitting the Linear Model
fit <- lmFit(mVals, design)
str(design)
# Contrast Matrix (Case vs Control)
contMatrix <- makeContrasts(CaseVsControl = S_Group1 - S_Group2, levels = design)
contMatrix
# Applying the Contrast
fit2 <- contrasts.fit(fit, contMatrix)
fit2 <- eBayes(fit2)
# Match the rownames of m_vals to the IlmnID column of annEPICv2 to extract relevant annotation
annEPICv2Sub <- annEPICv2[match(rownames(mVals), annEPICv2$IlmnID),]
# Extract top differentially methylated positions (DMPs)
DMPs1 <- topTable(fit2, num=Inf, coef=1, genelist=annEPICv2Sub)
head(DMPs1)
Yesterday, I realized that the approach I was using with the limma package was performing a pairwise comparison. To address this, I used the following code:
design <- model.matrix(~ S_Group + age + gender , data=sample_sheet)
fit <- lmFit(mVals, design)
fit <- eBayes(fit)
topTable(fit)
Thank you