I'm starting to use Limma to find differential expressed sites within my data and would like to ask if my approach was correct.

The data has a control and a treatment set, each containing data of 48 patients, which were measured before treatment (t0) and after treatment (t1). So I have Control-t0 and Control-t1 as well as Treatment-t0 and Treatment-t1, whereas Control-t0 and Treatment-t0 should not show any effect. Now I would like to know if there is a significant change in Treatment compared to Control. Therefore, I think, I have to adjust for the baseline (t0), right?

The code below is probably the most verbose option Limma provides, but I wanted to use the long method first before I use the shortcuts. So would the Interaction-coefficient I defined within my cont-matrix be the difference between Control and Treatment adjusted for the baseline? Also, I had to turn of adjusting for multiple testing here as my example data is random

#Data.frame with columns representing samples and rows different sites/genes
data_df <- data.frame(replicate(8, sample(100,1000, rep=TRUE)))
colnames(data_df) <- c("Patient_01_t0", "Patient_01_t1",
                       "Patient_02_t0", "Patient_02_t1",
                       "Patient_03_t0", "Patient_03_t1",
                       "Patient_04_t0", "Patient_04_t1")
#Metadata for data_df
meta <- data.frame(Sample= colnames(data_df),
                   ID= c(rep("Patient_01",2),
                   Time= c(rep(c("t0","t1"),4)),
#Create levels for model-matrix
time <- factor(meta$Time, labels=c('t0','t1') )
group<- factor(meta$Group, labels=c('Control','Treatment') )

design <- model.matrix(~0+lev)

#Create contrast matrix with all factors of interest
cont_matrix <- makeContrasts( t1_VS_t0_Ctrl = levControl.t1 - levControl.t0,
                              t1_VS_t0_Treat = levTreatment.t1 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t0= levControl.t0 - levTreatment.t0,
                              Ctrl_vs_Treat_at_t1= levControl.t1 - levTreatment.t1,
                              Interaction = (levControl.t1-levControl.t0)-(levTreatment.t1-levTreatment.t0), 

#Calculate different expressed sites
fit <- lmFit(data_df, design) 
fit <- contrasts.fit(fit, cont_matrix) 
fit.eBayes <- eBayes(fit, robust=TRUE)

summary(decideTests(fit.eBayes, adjust.method="none"))

1 Answer 1


It looks good to me. One think I could suggest is to see the distribution of your data/model (with limma/voom).

After you get the list of your genes with:

top.table <- topTable(fit.eBayes, sort.by = "P", n = Inf)

Maybe you can check your multidimensional scaling (MDS):

plotMDS(your_data, col = as.numeric(group))

Then you can also check with voom

y <- voom(data, model.matrix, plot = T)

I hope this helps

  • $\begingroup$ Thank you for your answer. Checking the distro might be very useful, indeed! $\endgroup$
    – Lukas
    Commented Sep 23, 2021 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.