Help with Limma-model

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

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)



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

• Thank you for your answer. Checking the distro might be very useful, indeed! Sep 23 at 21:22