# Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design

I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in in vitro cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps.

> library(edgeR)
> y <- DGEList(counts = KCraw.data[,2:25])
> keep <- rowSums(cpm(y) > 10) >= 3
> targets <- data.frame(rownames=colnames(KCraw.data)[2:25] ,
+                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
+                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
+                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
> Group <- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep=".")) > targets <- cbind(targets,Group=Group) > targets rownames Time FieldN CultureN Group 1 KC1_H2L 2DIC FH CL FH.2DIC.CL 2 KC2_H2L 2DIC FH CL FH.2DIC.CL 3 KC3_H2L 2DIC FH CL FH.2DIC.CL 4 KC4_H2H 2DIC FH CH FH.2DIC.CH 5 KC5_H2H 2DIC FH CH FH.2DIC.CH 6 KC6_H2H 2DIC FH CH FH.2DIC.CH 7 KC7_L2L 2DIC FL CL FL.2DIC.CL 8 KC8_L2L 2DIC FL CL FL.2DIC.CL 9 KC9_L2L 2DIC FL CL FL.2DIC.CL 10 KC10_L2H 2DIC FL CH FL.2DIC.CH 11 KC11_L2H 2DIC FL CH FL.2DIC.CH 12 KC12_L2H 2DIC FL CH FL.2DIC.CH 13 KC13_H5L 5DIC FH CL FH.5DIC.CL 14 KC14_H5L 5DIC FH CL FH.5DIC.CL 15 KC15_H5L 5DIC FH CL FH.5DIC.CL 16 KC16_H5H 5DIC FH CH FH.5DIC.CH 17 KC17_H5H 5DIC FH CH FH.5DIC.CH 18 KC18_H5H 5DIC FH CH FH.5DIC.CH 19 KC19_L5L 5DIC FL CL FL.5DIC.CL 20 KC20_L5L 5DIC FL CL FL.5DIC.CL 21 KC21_L5L 5DIC FL CL FL.5DIC.CL 22 KC22_L5H 5DIC FL CH FL.5DIC.CH 23 KC23_L5H 5DIC FL CH FL.5DIC.CH 24 KC24_L5H 5DIC FL CH FL.5DIC.CH  I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example > y <- DGEList(counts = KCraw.data[keep,2:25], group = Group) > y <- calcNormFactors(y) > > TMM <- KCraw.data[keep,2:25] > for (i in 1:24) { + TMM[,i] <- TMM[,i] / (y$$samples$$lib.size[i] * y$$samples$$norm.factors[i]) * 1e6 + } > > y <- DGEList(counts = TMM,group = Group) > > design <- model.matrix(~0+Group) > colnames(design) <- levels(Group) > y <- calcNormFactors(y,method = "TMM") > y <- estimateDisp(y,design) > fitQL <- glmQLFit(y,design) > fit <- glmFit(y,design) > myKC.contrasts <- makeContrasts( + H2H.H2L = FH.2DIC.CH - FH.2DIC.CL, + L2H.L2L = FL.2DIC.CH - FL.2DIC.CL, + H2H.L2H = FH.2DIC.CH - FL.2DIC.CH, + H2L.L2L = FH.2DIC.CL - FL.2DIC.CL, + H5H.H5L = FH.5DIC.CH - FH.5DIC.CL, + L5H.L5L = FL.5DIC.CH - FL.5DIC.CL, + H5H.L5H = FH.5DIC.CH - FL.5DIC.CH, + H5L.L5L = FH.5DIC.CL - FL.5DIC.CL, + H2H.L2L = FH.2DIC.CH - FL.2DIC.CL, + H5H.L5L = FH.5DIC.CH - FL.5DIC.CL, + H5L.H2L = FH.5DIC.CL - FH.2DIC.CL, + H5H.H2H = FH.5DIC.CH - FH.2DIC.CH, + L5L.L2L = FL.5DIC.CL - FL.2DIC.CL, + L5H.L2H = FL.5DIC.CH - FL.2DIC.CH, + levels=design) > design FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL 1 0 1 0 0 0 0 0 0 2 0 1 0 0 0 0 0 0 3 0 1 0 0 0 0 0 0 4 1 0 0 0 0 0 0 0 5 1 0 0 0 0 0 0 0 6 1 0 0 0 0 0 0 0 7 0 0 0 0 0 1 0 0 8 0 0 0 0 0 1 0 0 9 0 0 0 0 0 1 0 0 10 0 0 0 0 1 0 0 0 11 0 0 0 0 1 0 0 0 12 0 0 0 0 1 0 0 0 13 0 0 0 1 0 0 0 0 14 0 0 0 1 0 0 0 0 15 0 0 0 1 0 0 0 0 16 0 0 1 0 0 0 0 0 17 0 0 1 0 0 0 0 0 18 0 0 1 0 0 0 0 0 19 0 0 0 0 0 0 0 1 20 0 0 0 0 0 0 0 1 21 0 0 0 0 0 0 0 1 22 0 0 0 0 0 0 1 0 23 0 0 0 0 0 0 1 0 24 0 0 0 0 0 0 1 0 attr(,"assign") [1] 1 1 1 1 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$Group
[1] "contr.treatment"

> myKC.contrasts
Contrasts
Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
Contrasts
Levels       L5H.L2H
FH.2DIC.CH       0
FH.2DIC.CL       0
FH.5DIC.CH       0
FH.5DIC.CL       0
FL.2DIC.CH      -1
FL.2DIC.CL       0
FL.5DIC.CH       1
FL.5DIC.CL       0


After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.

> myKC.contrasts <- cbind(myKC.contrasts,
+                         Development = c(1,1,-1,-1,1,1,-1,-1),
+                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
+                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
+ )
> myKC.contrasts
H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
L5H.L2H Development FieldN CultureN
FH.2DIC.CH       0           1      1        1
FH.2DIC.CL       0           1      1       -1
FH.5DIC.CH       0          -1      1        1
FH.5DIC.CL       0          -1      1       -1
FL.2DIC.CH      -1           1     -1        1
FL.2DIC.CL       0           1     -1       -1
FL.5DIC.CH       1          -1     -1        1
FL.5DIC.CL       0          -1     -1       -1


Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast.

> lrt <- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
> topTags(lrt,n=nrow(y\$counts))["GRMZM2G445575",]
Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL
logFC   logCPM        F       PValue          FDR
GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
# FC is a data frame of the logFC of each constrast in columns for each gene in rows
> sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
[1] -6.636197


My first question is if this analysis is a valid way of summarizing the simple effects of each treatment. I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors.

• I don't quite get why you need to manually change the counts with a for loop. It is used as a size factor in the linear model. Now, for the pairwise comparison, what you have in the first part is ok – StupidWolf Jun 16 at 21:26
• For the general culture N etc, run a linear model and run the estimateGLM... etc stuff, it should give you the effects you want. You cannot get everything from one model – StupidWolf Jun 16 at 21:27

• Ideally you could use a package as reprex. But if you use code blocks and remove the > at the beginning of the code, comment the output (so that if I paste it on the console it doesn't error) you'll do good. Another thing to consider is that I cannot load("KC_Raw.RData") because I don't have it, so try to make the example self contained. But the code provided is enough to understand your problem so don't worry if you use the same code used here there. – llrs Jun 18 at 9:29