How to model a pre and post treatment in limma
Entering edit mode
AST ▴ 50
Last seen 7 weeks ago

I have RNA-seq data that has three groups - Control (C), Disease1 (D1), and Disease2 (D2). The patients of D1 and D2 were treated with a drug (T) and data was generated on pre- and post-treatment disease groups i.e. I have 5 grouping factors - C, D1-pre, D1-post, D2-pre, D2-post. I am a bit confused while modeling this nested data to identify DE genes in limma. I have two questions in mind:

  1. What are DE genes in D1 and D2 compared to control, and
  2. Whether the drug has a differential effect on the expression of DE genes in D1 compared to D2.

I tried to make a model matrix as:

patient <- data.frame(
    patient$patient <- factor(paste(rep(1:5, times = 5),"_",patient$diagnosis))
    diagnosis = factor(rep(c("D1", "D2","Control"), each=10, length.out=25)),
    treat = factor(rep(c("Pre", "Post"), each=5, times=3, length.out=25)))
patient$group <- factor(paste(patient$diagnosis,patient$treat))
design <- model.matrix(~0+group,patient)
colnames(design) <- c("C", "D1post", "D1pre", "D2post", "D2pre")

I am lost somewhere here, I am unable to decide if should I go for

fit <- lmFit(eset, design) #eset is the data matrix with expression levels

or make a contrast model, something like:

contrast.matrix <- makeContrasts(D1pre-C, D2pre-C, D1post-D1pre, D2post-D2pre, levels=design)
fit2 <-, contrast.matrix)

I will be thankful if someone could please help me in modeling.

limma • 294 views
Entering edit mode

Do you have pre and post samples from the same patients? In other words, do the D1-pre and D1-post groups contain the same patients or different patients?

Entering edit mode

Yes, 'pre' and 'post' patients in each disease group are the same. I am sorry for my laziness. I have updated the same in my question. Thank you for replying on Sunday. Now that you have pointed out, can I use these patients' samples per group as random errors?

Entering edit mode
Last seen 2 hours ago
WEHI, Melbourne, Australia

This seems to be a multi-level experiment, so you could follow the limma User's Guide on that topic.

Entering edit mode

I read the user guide. However, since I am not looking for a pair-wise comparison, I am a bit confused with the modeling. I am thinking in lines, something like:

des1 <- design[,c(1:3)]
fit <- lmFit(eset, des1)
fit.e <- eBayes(fit, trend=TRUE, robust=TRUE)
topTable(fit, coef=2)

and a similar model for D2. I am not sure how correct these models will be to identify 'true' DEGs in post-treatment groups for both diseases.


Login before adding your answer.

Traffic: 423 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6