ANOVA Limma
1
1
Entering edit mode
georgina.fqw ▴ 10
@georginafqw-23788
Last seen 3.1 years ago

Hello,

I would like to run ANOVA on my seq data using limma. I have paired samples (normal vs disease) and 3 different types of disease A,B,C. I would like to use block=patient_block for paired design. Here is my code so far. Thank you for help.

Create DGEList object

y <- DGEList(counts=counts, samples=samples)
keep <- filterByExpr(y, group=y$samples$Disease)
disease <- y[keep,, keep.lib.sizes=FALSE]

Normalization

disease <- calcNormFactors(disease, method = "TMM")

Block factor: patient

patient_block<-as.factor(as.character(samples$Patient))
model <- model.matrix(~0 + Disease+Age+Sex, data=samples)
z <- voom(disease, model, plot = T)
dupcor<-duplicateCorrelation(z,model,block=patient_block)
dupcor$consensus.correlation

Fit a mixed effect model using limma

fit <- lmFit(z,model,block=patient_block,correlation=dupcor$consensus.correlation)
fit2 <- eBayes(fit)
topTable(fit2)

Where samples contains Disease information (either A, B, C or normal), Patient ID, Age, Sex And counts contains matrix count (RNAseq data)

and ANOVA??? for comparison AvBvCvNormal

topTable(fit2, coef=c("DiseaseA", "DiseaseB", "DiseaseC", "DiseaseN"))

OR should I use

contrast.matrix<-makeContrasts(Normal_A=Normal-DiseaseA, Normal_B=Normal-DiseaseB, Normal_C=Normal-DiseaseC, levels=colnames(coef(fit)))

fit3 <- contrasts.fit(fit, contrast.matrix)

fit3 <- eBayes(fit3)

summa.fit3 <- decideTests(fit3)

topTable(fit3, coef=c("Normal_A", "Normal_B", "Normal_C")

which one, if any is correct?

ANOVA limma • 1.2k views
ADD COMMENT
1
Entering edit mode
@james-w-macdonald-5106
Last seen 24 minutes ago
United States

You should probably do this instead:

fit <- voomLmFit(y, model, patient_block)
fit2 <- eBayes(fit)

Which will deal with estimating the correlation structure and weights better than you are (you have to iterate).

And this is wrong

topTable(fit2, coef=c("DiseaseA", "DiseaseB", "DiseaseC", "DiseaseN"))

because you have fit a cell means model, and those coefficients simply estimate the mean expression for each group. Which is an uninteresting quantity. Instead you want to know the difference between groups, for which you always have to use contrasts.fit if you have a cell means model (because it makes the between groups comparisons).

ADD COMMENT
0
Entering edit mode

Hi James,

Thank you for your response. So if I understand correctly:

fit <- voomLmFit(disease, model, patient_block)
contrast.matrix<-makeContrasts(Normal_A=Normal-DiseaseA, Normal_B=Normal-DiseaseB, Normal_C=Normal-DiseaseC, levels=colnames(coef(fit)))
fit <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit)
topTable(fit2, coef=c("Normal_A", "Normal_B", "Normal_C")

Is that what you mean?

ADD REPLY
0
Entering edit mode

Oh yeah, sorry. You do want the DGEList with the normalization factors. What you are doing is what I meant to say you should do.

ADD REPLY

Login before adding your answer.

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