Hi,
My experiment has the following design:
- 2 groups of patients
- pre-/post-treatment (patient matched)
I think the contrasts.fit/topTable works, but I am not sure how to test gene set using roast under limma with such design. I couldn't work out how to test specific comparison or take into account of patient-matching data. It seems to give me highly significant FDR even though the genes have values are more or less the same across all samples.
Here are some dummy data showing the structure of my data/experimental design and related contrast-fit and roast r script:
expression_data
| sample_1 | sample_2 | sample_3 | sample_4 | sample_5 | sample_6 | sample_7 | sample_8 | sample_9 | sample_10 | sample_11 | sample_12 | |
| Gene1 | 5.7 | 5.8 | 6.1 | 6 | 5.9 | 6 | 6.2 | 5.9 | 6 | 6 | 5.7 | 5.7 |
| Gene2 | 3.4 | 3.9 | 3.6 | 4.2 | 3.4 | 3.9 | 7.3 | 7.6 | 7.2 | 7.6 | 7 | 7.5 |
| Gene3 | 5.6 | 5.5 | 6.2 | 5.9 | 6.3 | 6.3 | 6.5 | 5.2 | 6.4 | 6.4 | 6 | 5.7 |
| Gene4 | 3.7 | 3.6 | 3.6 | 3.9 | 3.4 | 3.5 | 3.5 | 3.5 | 3.6 | 3.6 | 3.6 | 3.6 |
| Gene5 | 6.8 | 6.6 | 6.9 | 6.9 | 5.9 | 6.7 | 6.3 | 5.7 | 6.8 | 6.9 | 6.6 | 6.7 |
| Gene6 | 9.5 | 9.6 | 10.4 | 10.1 | 10.4 | 9.9 | 1.5 | 1.5 | 1.3 | 1.6 | 2.1 | 2.2 |
| Gene7 | 6.3 | 6.6 | 6.5 | 6.5 | 6.5 | 6.5 | 6.6 | 6.4 | 6.5 | 6.4 | 6.6 | 6.3 |
| Gene8 | 1.2 | 1.5 | 1.3 | 1.8 | 1.2 | 1.6 | 5 | 5.4 | 5 | 5.2 | 4.8 | 5.3 |
| Gene9 | 8.3 | 8.4 | 9.4 | 8.8 | 9.3 | 8.6 | 9.3 | 8.5 | 9.3 | 9.1 | 8.5 | 8.5 |
| Gene10 | 1.4 | 1.9 | 1.6 | 2.2 | 1.4 | 1.9 | 5.3 | 5.6 | 5.2 | 5.6 | 5 | 5.5 |
library(limma)
expression_data <- read.delim("expression_data.txt", row.names=1)
study_design <- data.frame(subjectID=c("S1", "S1", "S2", "S2", "S3", "S3", "S4", "S4", "S5", "S5", "S6", "S6"), group=c(rep("Group1", 6), rep("Group2", 6)), treatment=rep(c("ctrl", "treated"), 6))
exprs_mx <- as.matrix(expression_data)
rowMatch <- c(2, 8, 10)
rowMatch_similarValuesAcrossSamples <- c(1, 3, 4, 5)
patientID <- factor(study_design$subjectID)
grp <- factor(study_design$group)
treatment <- factor(study_design$treatment)
grpTreat <- factor(paste(grp, treatment, sep="."))
design <- model.matrix(~0+ grpTreat)
colnames(design) <- levels(grpTreat)
corfit <- duplicateCorrelation (exprs_mx, design, block=patientID)
fit <- lmFit(exprs_mx, design, block=patientID, correlation=corfit$consensus)
groupContrast <- makeContrasts(
group1Treated = Group1.treated - Group1.ctrl,
group2Treated = Group2.treated - Group2.ctrl,
groupDiffTreated = Group2.treated - Group1.treated,
groupDiffBaseline = Group2.ctrl - Group1.ctrl,
levels = design
)
fit2 <- contrasts.fit (fit,groupContrast)
fit2 <- eBayes(fit2)
topTable(fit2, coef = 4) # coef refers to the comparison in groupContrast. i.e. coef = 4 refers to groupDiffBaseline
roast_result <- roast(exprs_mx, rowMatch, design, block=patientID, correlation=corfit$consensus)
roast_resultRandom <- roast(exprs_mx, rowMatch_similarValuesAcrossSamples, design, block=patientID, correlation=corfit$consensus)
I am wondering if anyone know how to use roast on study designs like this?
Thank you very much for your help.
Barbara
Thanks for your reply. I think I am making a mistake somewhere, because when I used contrast & block, everything came up highly significant - I think it is comparing 1 individual to everyone else, but I don't know for sure since there are just 4 comparisons in the design.
Any idea what else I might be able to try?
Thank you