Dear all, 

Let me ask you all a question. 

If I want to create a volcano plot after using the treat function, What kind of parameters should I use?  Of note, in this case, all genes do not rank using Log2fold change. Should I use  xlab="treat.lfc" and ylab=" ? " . Could anyone please help me?

Thanks in advance.

`> names(tfit)

"coefficients"    "rank"    "assign"    "qr"     "df.residual"

"sigma"   "stdev.unscaled"    "design" "correlation" "genes"

"Amean"    "cov.coefficients"  "contrasts"   "df.prior"   "s2.prior"

"s2.post"    "df.total"   "t"     "p.value"     "treat.lfc" `

>volcanoplot(tfit,coef=1)

>plot(tfit[,6],-log(tfit[,10],10),coef=1, xlab="log2(Fold-Change)", ylab="-log10(P.Value)",main="Volcano plot",cex=0.1,pch=19)

>abline(v=c(-1,1),col="blue")

>abline(h=-log(0.05,10),col="red")

Thanks in advance.

Your question has made me realize that the volcanoplot() function doesn't currently work with treat(). There's no good reason why it shouldn't, so I'll commit a fix.

Note, the fact that p-values and logFCs give different rankings isn't the problem. If p-values and logFC ranked genes in the same order, then people wouldn't do volcano plots.

The code you give for making your own plot isn't correct. What information were you trying to access by tfit[,6] or tfit[,10]? What did you think this code would do?

If you wanted to make your own volcano plot, for contrast 3 for example, you would use:

logFC <- tfit$coefficients[,3]
pvalue <- tfit$p.value[,3]
plot(logFC, log10(pvalue))

If you wanted to make a volcano plot for contrast 1 you would change both the 3s to 1s.