Question: Limma to Compare Bulk RNA Seq using makeContrasts and eBayes
0
gravatar for crystalize123
9 months ago by
crystalize1230 wrote:

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 limma bioconductor R • 204 views
ADD COMMENTlink modified 9 months ago by h.mon31k • written 9 months ago by crystalize1230
1
gravatar for Papyrus
9 months ago by
Papyrus540
Papyrus540 wrote:

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 COMMENTlink written 9 months ago by Papyrus540
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1211 users visited in the last hour