Limma to Compare Bulk RNA Seq using makeContrasts and eBayes
1
0
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(), contrasts.fit(), eBayes(), all of which are in the limma package.



Here is my minimal working example.

library(limma)

# 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( contrasts.fit( fit, contrast.matrix ) )

# Want to run but cannot

summary(decideTests(fit2))


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.

Thanks!

RNA-Seq R Limma Bioconductor • 579 views
ADD COMMENT
1
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)

ADD COMMENT

Login before adding your answer.

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