Help with Limma-model

Question

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),
                         rep("Patient_02",2),
                         rep("Patient_03",2),
                         rep("Patient_04",2)),
                   Time= c(rep(c("t0","t1"),4)),
                   Group=c(rep("Control",4),
                           rep("Treatment",4))
)
                         
#Create levels for model-matrix
time <- factor(meta$Time, labels=c('t0','t1') )
group<- factor(meta$Group, labels=c('Control','Treatment') )
lev=paste(group,".",time,sep="")
lev


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

#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), 
                              levels=design)

#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"))

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