Hi

I think this is correct, but it would be good if some one could just verify this. I am comparing a group of patients before and after vaccination at a number of time points using limma.

PATIENTS <- factor(des$PATIENT_ID) # this is just a list of patient IDs

design <- model.matrix(~0+des$VISIT+PATIENTS) # visit refers to the patient visit number or time point
colnames(design)[1:5] <- c('Visit1', 'Visit2', 'Visit4', 'Visit6', 'Visit8')
matrix <- data.matrix(data3) # this is array data
fit <- lmFit(matrix, design)
contrast.matrix <- makeContrasts(Visit1-Visit2, Visit2-Visit4, Visit4-Visit6, Visit6-Visit8, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
top2 <- topTable(fit2,coef=1,number=Inf,sort.by="P")

Then I am just varying the coef parameter to select different contrasts.

Thanks,

Chris

Looks okay to me. Note that your contrasts are comparing samples from sequential visits, rather than comparing the sample from each visit to the "before vaccination" sample (which I would presume to be that from the first visit). If this is intentional, then it's fine. If you want to get any changes between visits, then you can just set coef=NULL in your final topTable call to do an ANOVA-like contrast.

I am having a problem with the contrast code I posted previously. For coeff=4 i.e. a comparison between visits 6 and visits 8 ALL of my genes or probes are differentially expressed. However, when I just select these visits manually then do a normal paired t-test without contrasts I get no differentially expressed genes. So what is this about? Any clues? I know they should be different due to the variance being shared with the contrast test, but this seems strange. The other contrasts look normal, however.

# This is code for a standard paired t test w/o using contrasts

des2 <- subset(des2, des2$VISIT %in% c('Vaccine Visit 6', 'Vaccine Visit 8'))
tokeep <- as.character(des2$BARCODE)
data2 <- data[,tokeep]
 
PATIENTS <- factor(des2$PATIENT_ID)
 
#set up design
 
design <- model.matrix(~des2$VISIT+PATIENTS)
data3 <- data2
matrix <- data.matrix(data3)
fit <- lmFit(matrix,design)
fit <- eBayes(fit)

top2 <- topTable(fit,coef=2,number=Inf,sort.by="P")