Can I add more than one effect (example: treatment and sex) to cuffdiff ?

Question

[TCC]

I am having a problem to analyze the differentially expressed genes by using cuffdiff , because besides the two treatments I also have samples of different sex, and I would like to distinguish beteween them ?

Is it possible to run a 2 step approach using cuffdiff (first with sex effect and repeat with treatemnt) ? If so , which files should I load in the program?

Which alternative programs should I use after aligment to analyze this model with two effects (treatment and sex).

Thanks

[ATpoint]

Do you realize that cuffdiff is a tool for differential transcript rather than gene level analysis? Better use the tools Charles Warden mentions below. Only recommendation I have when using edgeR is to use the QLF pipeline (check edgeR manual) which is the current recommendation of the authors, claiming to be slightly superior to the LRT pipeline.

Answer 1

[Charles Warden]

I don't believe that you can include a second variable for cuffdiff.

However, you can do so for edgeR, limma-voom, DESeq2, and even a regular ANOVA or linear regression.

In fact, you might find those strategies advantageous over cuffdiff (even though cufflinks can be useful for reference guided de novo assemblies):

http://cdwscience.blogspot.com/2019/11/requiring-at-least-some-methods-testing.html

In terms of the relevant commands:

edgeR:

design = model.matrix(~var1 + var2)
fit = glmFit(y, design)
lrt = glmLRT(fit, coef=2)
test.pvalue = lrt$table$PValue

(you can also add y = estimateGLMRobustDisp(y, design) for edgeR-robust)

limma-voom:

design = model.matrix(~var1 + var2)
v = voom(y,design,plot=TRUE)
fit = lmFit(v,design)
fit = eBayes(fit)
pvalue.mat = data.frame(fit$p.value)

DESeq2:

dds = DESeqDataSetFromMatrix(countData = deg.counts, colData = colData, design = ~ var1 + var2)

You can then calculate a p-value via Wald Test (dds = DESeq(dds) followed by res = results(dds, contrast = c("var1", trt.group, other.groups))) or Likelihood Ratio Test (dds = DESeq(dds, test="LRT", reduced = ~ var2) followed by res = results(dds)).