Question: Can I add more than one effect (example: treatment and sex) to cuffdiff ?
gravatar for TCC
25 days ago by
TCC0 wrote:

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


rna-seq • 110 views
ADD COMMENTlink modified 25 days ago by Charles Warden7.9k • written 25 days ago by TCC0

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.

ADD REPLYlink written 25 days ago by ATpoint40k
gravatar for Charles Warden
25 days ago by
Charles Warden7.9k
Duarte, CA
Charles Warden7.9k wrote:

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

In terms of the relevant commands:


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)


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)


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",, other.groups))) or Likelihood Ratio Test (dds = DESeq(dds, test="LRT", reduced = ~ var2) followed by res = results(dds)).

ADD COMMENTlink written 25 days ago by Charles Warden7.9k
Please log in to add an answer.


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