Is this the correct way of doing a paired t-test in limma using contrasts
2
1
Entering edit mode
chris86 ▴ 420
@chris86-8408
Last seen 4.3 years ago
UCL, United Kingdom

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

 

 

 

 

limma • 2.5k views
ADD COMMENT
2
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 17 hours ago
The city by the bay

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.

ADD COMMENT
1
Entering edit mode

Hi Chris,

just as a side note: instead of

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

you can also specify the name of the contrasts directly like so:

top2 <- topTable(fit2,coef="Visit1-Visit2",number=Inf,sort.by="P")

I usually prefer to do it this way, as it avoids mistakes and makes the code more readable.

 

 

 

ADD REPLY
0
Entering edit mode
chris86 ▴ 420
@chris86-8408
Last seen 4.3 years ago
UCL, United Kingdom

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

 

ADD COMMENT
1
Entering edit mode

The original setup will have more power because you've got more residual d.f. to estimate the variance (and thus more certainty for downstream inferences). This would explain why you get more DE genes compared to the new analysis above. However, it's surprising that all of the genes are DE in your original analysis. This shouldn't occur if you've normalised correctly, as systematic differences in expression between samples should be normalised out. Alternatively, if you have normalised and you still get this problem, perhaps there's just a lot of genes with very small log-fold changes; try using treat to set a log-fold change threshold to get only those genes with large changes in expression.

ADD REPLY
0
Entering edit mode

Hi, sorted it. Actually the data had some quirks I did not realise and the design matrix was messed up as a result. Thanks.

ADD REPLY

Login before adding your answer.

Traffic: 964 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6