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