Hello everyone,

I try to implement LIMMA for paired samples in order to compare gene expression before/after treatment. But...i'm not confident in my results since every single gene expression seems to be significantly associeted with the treatment!

Here is my methodology:

mat<-t(df)

where df is a dataframe with: -1000 columns for different gene expression -+ 1 column 'id' (10 different id, repeated 2 times) -+ 1 column 'ttt' (0=before, 1=after treatment) -20 rows of observations (10 patients observed twice)

So once I get my matrix:

subject_id<-df$id A vector with 10 different id repeated 2 times

ttt<-as.factor(df$ttt) A vector with the treatment status (before/after)

Here is my design matrix:

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

colnames(design) <- levels(ttt)

As asked in limma() description i specify the correlations:

corfit <- duplicateCorrelation(mat,design,block=subject_id)

Fitting of the model:

fit<-lmFit(mat,design,block =subject_id,correlation = corfit$consensus)

fit2<-eBayes(fit)

Getting adjusted pvalues for every gene comparison before/after treatment:

topTable(fit2, adjust.method='BH',number=1000)

But every single gene expression seems to be relevant, so i'm very doubtful about my methodology Any idea? Many thanks!

The model you have fit (without an intercept) simply computes the mean for each group, which isn't an interesting measure. Instead you want to know if there are any differences between groups, so you need to construct a contrasts matrix and use contrasts.fit to make the contrasts (before the eBayes step).

You can either use makeContrasts to create the contrasts matrix, or just generate one yourself, using examples in the limma User's Guide to help figure out how one does such a thing.