Limma to Compare Bulk RNA Seq using makeContrasts and eBayes
Entering edit mode
2.2 years ago

After a day of googling, I've decided that it'd be better to ask the question here.

So the experiment is I have bulk RNA seq data from 3 patients: A, B, C. And their RNA seq data is obtained for pre-treatment, treatment cycle 1, treatment cycle 2, treatment cycle 3.

So in total I have 12 samples of bulk RNA seq:

  • A.PreTreat -> A.Cycle1 -> A.Cycle2 -> A.Cycle3

  • B.PreTreat -> B.Cycle1 -> B.Cycle2 -> B.Cycle3

  • C.PreTreat -> C.Cycle1 -> C.Cycle2 -> C.Cycle3

I want to get a differential gene list between different cycles (i.e. cycle 3 to pretreatment, cycle 3 to cycle 2) using model.matrix(), lmFit(), makeContrasts(),, eBayes(), all of which are in the limma package.

Here is my minimal working example.


# Already normalized expression set: rows are genes, columns are the 12 samples`

normalized_expression <- matrix(data=sample(1:100), nrow=10, ncol=12)
colnames(normalized_expression) <- c("A.PreTreat", "A.Cycle1", "A.Cycle2", "A.Cycle3", 
                                     "B.PreTreat", "B.Cycle1", "B.Cycle2", "B.Cycle3", 
                                     "C.PreTreat", "C.Cycle1", "C.Cycle2", "C.Cycle3")

patient_and_treatment <- factor(colnames(normalized_expression), levels = colnames(normalized_expression))

design.matrix <- model.matrix(~0 + patient_and_treatment)
colnames(design.matrix) <- patient_and_treatment
fit <- lmFit(normalized_expression, design.matrix)

# I want to get a contrast matrix to get differential genes between cycle 3 treatment and pre-treatment in all patients

contrast.matrix <- makeContrasts("A.Cycle3+B.Cycle3+C.Cycle3-A.PreTreat-B.PreTreat-C.PreTreat",
                                levels = levels(patient_and_treatment))

# Outputs Error of no residual degree of freedom

fit2 <- eBayes( fit, contrast.matrix ) )

# Want to run but cannot


So far I am stuck on no residual degree of freedom error.

I am not even sure if this is the statistically right way in limma to address my question of getting differential gene list between cycle 3 treatment to pre-treatment in all patients.

Any help will be greatly appreciated.


RNA-Seq R Limma Bioconductor • 579 views
Entering edit mode
2.2 years ago
Papyrus ★ 2.1k

If you're analyzing RNA-seq with limma take into account that the recommended workflow is to use the voom method (or eBayes with trend = T)


Login before adding your answer.

Traffic: 856 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