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?

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).