Hi, I am analyzing RNA-seq data of a tree infected with 2 pathogens. The experiment lay-out is like this:

• condition 1(C): tree + water (control)
• condition 2(P1): tree + pathogen1
• condition 3(P1_P2): tree + pathogen1 + pathogen2

I inherited this code and I succeeded using it to call DEGs for condition 2 (contrast = 2-1; 2401 DEGs) and 3 (contrast = 3-1; 5209 DEGs), but when I tried calling DEGs for differences between condition 2 and 3 (contrast = either 2-3 or 3-2), the result is always 0. I expect there are some differences between the two conditions given that, when compared to control, the number of DEGs found in the two conditions are quite different. May you please help enlighten if I have use wrong codes or this is what I can expect if there is no difference between the two conditions?

I also wonder if I really need the line "colnames(design) <- levels(group)" in my code. I tried reading up about a function of colnames, but I don't understand what it actually does.

All suggestions and advice are appreciated. Thanks a lot in advance for your help! :)

Code should be placed in three backticks as shown below

> rawCountTable <- read.table("counts.txt", header=TRUE, sep="\t", row.names="Geneid")
> sampleInfo <- read.table("sampleInfo.csv", header=TRUE, sep=",", row.names="file")
> group <- factor(sampleInfo$condition)
> y <- DGEList(counts=rawCountTable,group=sampleInfo$condition)
> keep <- filterByExpr(y)
> y <- y[keep,,keep.lib.sizes=FALSE]
> y <- calcNormFactors(y, method="TMM")
> design <- model.matrix(~0+group)
$design
         C  P1  P1_P2
1       1   0       0
2       1   0       0
3       1   0       0
4       0   1       0
5       0   1       0
6       0   0       1
7       0   0       1
8       0   0       1
> colnames(design) <- levels(group)
> y <- estimateDisp(y,design,robust=TRUE)
> y$common.dispersion
> fit <- glmQLFit(y,design,robust=TRUE)
> my.contrasts <- makeContrasts(vsP1=P1-C, vsP1_P2=P1_P1-C, vsdiff=(P1_P2-C)-(P1-C), vsdiff2=(P1-C)-(P1_P2-C)), levels=design)
> my.contrasts
  Contrasts
Levels        P1_P2    P1 vsdiff vsdiff2
  control        -1    -1      0       0
  P1              0     1     -1       1
  P1_P2           1     0      1      -1
> qvsdiff <- glmQLFTest(fit, contrast=my.contrasts[,"vsdiff"])
> toptags.qvsdiff <- topTags(qvsdiff, n=nrow(qvsdiff$table),adjust.method="fdr")
> sig005.qvsdiff <- toptags.qvsdiff$table[toptags.qvsdiff$table$FDR<0.05,] 
> sum(toptags.qvsdiff$table$FDR < 0.05)
0

I expect there are some differences between the two conditions given that, when compared to control, the number of DEGs found in the two conditions are quite different.

No, that reasoning doesn't follow. It is perfectly possible to get more DE genes for P1_P2 vs C than P1 vs C but to still get no DE genes for P1_P2 vs P1. It may be that there are lots of genes changing in the same direction in both cases just making the FDR cutoff for P1_P2 vs C and just missing out for P1 vs C,.

May you please help enlighten if I have use wrong codes

Your contrast code is unnecessarily complicated but it all looks correct.

or this is what I can expect if there is no difference between the two conditions?

Yes, this what would be expected if P1_P2 and P1 are both very different to C but similar to one another. It would be logical to explore this visually with plotMDS(). See the edgeR workflow for a case study of a complete edgeR analysis with diagnostic plots.

I also wonder if I really need the line colnames(design) <- levels(group) in my code.

Yes you do, if you want to use the column headings C, P1 and P1_P2 when forming contrasts.

I tried reading up about a function of colnames, but I don't understand what it actually does

It simply sets the columns names. Just looking at the results would immediatel show you what it does:

> design <- model.matrix(~0+group)
> design
> colnames(design) <- levels(group)
> design