Dear all,

I am analyzing 34 samples of a shRNA screen using edgeR. There are 7450 shRNAs targeting 668 genes.
The experiment is composed of 7 cell lines including 3 controls, 4 treatments and 2 controls & 2 time points. The 4 cell lines which are not controls can be mutated (Mut1 or Mut2).
There are at least 2 replicates for each condition.
Below is the description of the 34 samples:

SampleID    CellLines    Time    Drug    Mut
Sample1    Ctrl1_FH    T1    None    None
Sample2    Ctrl1_FH    T2    None    None
Sample3    Ctrl2_MF    T1    None    None
Sample4    Ctrl2_MF    T2    None    None
Sample5    Ctrl3_SC    T1    None    None
Sample6    Ctrl3_SC    T2    None    None
Sample7    NE4    T1    None    Mut2
Sample8    NE4    T2    None    Mut2
Sample9    NE4    T2    Drug1    Mut2
Sample10    NE4    T2    Drug2    Mut2
Sample11    NE4    T2    Drug3    Mut2
Sample12    NE4    T2    Non-irradiated    Mut2
Sample13    NE4    T2    Irrad    Mut2
Sample14    NE5    T1    None    Mut2
Sample15    NE5    T2    None    Mut2
Sample16    NE5    T2    Drug1    Mut2
Sample17    NE5    T2    Drug2    Mut2
Sample18    NE5    T2    Drug3    Mut2
Sample19    NE5    T2    Non-irradiated    Mut2
Sample20    NE5    T2    Irrad    Mut2
Sample21    NE6    T1    None    Mut1
Sample22    NE6    T2    None    Mut1
Sample23    NE6    T2    Drug1    Mut1
Sample24    NE6    T2    Drug2    Mut1
Sample25    NE6    T2    Drug3    Mut1
Sample26    NE6    T2    Non-irradiated    Mut1
Sample27    NE6    T2    Irrad    Mut1
Sample28    NE7    T1    None    Mut1
Sample29    NE7    T2    None    Mut1
Sample30    NE7    T2    Drug1    Mut1
Sample31    NE7    T2    Drug2    Mut1
Sample32    NE7    T2    Drug3    Mut1
Sample33    NE7    T2    Non-irradiated    Mut1
Sample34    NE7    T2    Irrad    Mut1

For simplicity, I do not use the Mut variable for the beginning.
I try to answer 5 questions:
1) time effect
2) Drug1 effect
3) Drug2 effect
4) Drug3 effect
5) Irrad effect
6) "congelation" effect/ non irradiated vs untreated at T2: for irradiation study, cells were frozen. Controls cells were frozen at the same time points as well to served as controls for irradiated cells. I would like to check the effect of congelation, it is not a biological question, just for information.

I am asking help for defining the contrasts that will allow answering these questions.

Here is the R code I use:

library("edgeR")

rawdata <- read.delim("data_34Samples",header=TRUE, stringsAsFactors=FALSE) 

y <- DGEList(counts=rawdata[,5:38], ,genes=rawdata[,1:4])
selr <- rowSums(cpm(y$counts)>10)>=4

y <- y[selr,]
y <- calcNormFactors(y)

# Set up design matrix for GLM
Time= factor(c("T1","T2","T1","T2","T1","T2" ,
               "T1","T2","T2", "T2", "T2", "T2", "T2",
               "T1","T2","T2", "T2", "T2", "T2", "T2",
               "T1","T2","T2", "T2", "T2", "T2", "T2",
               "T1","T2","T2", "T2", "T2", "T2", "T2"))
CellLines=factor(c("Ctrl1_FH","Ctrl1_FH", "Ctrl2_MF","Ctrl2_MF", "Ctrl3_SC","Ctrl3_SC",
                   "NE4", "NE4", "NE4","NE4", "NE4", "NE4", "NE4",
                   "NE5", "NE5", "NE5", "NE5", "NE5", "NE5", "NE5",
                   "NE6", "NE6", "NE6", "NE6", "NE6", "NE6", "NE6",
                   "NE7", "NE7", "NE7", "NE7", "NE7", "NE7", "NE7"))
Drug= factor(c( "None",  "None", "None", "None", "None", "None",
                "None", "None", "Drug1", "Drug2", "Drug3", "NI", "Irradiated",
                "None", "None", "Drug1", "Drug2", "Drug3", "NI", "Irradiated",
                "None", "None", "Drug1", "Drug2", "Drug3", "NI", "Irradiated",
                "None", "None", "Drug1", "Drug2", "Drug3", "NI", "Irradiated"))
Drug <- relevel(Drug, "None")

design <- model.matrix(~0+CellLines+Time+Drug)
rownames(design) <- colnames(rawdata[,5:38])

# Estimate dispersions
yglm <- estimateDisp(y, design)
sqrt(yglm$common.disp)

# Fit negative bionomial GLM
fit <- glmFit(yglm, design)
colnames(design)

[1] "CellLinesCtrl1_FH" "CellLinesCtrl2_MF"   "CellLinesNE4"         "CellLinesNE5"         "CellLinesNE6"         "CellLinesNE7"        
[7] "CellLinesCtrl3_SC"     "TimeT2"                  "DrugDrug1"           "DrugDrug2"          "DrugIrradiated"          "DrugNI"                 
[13] "DrugDrug3

# Carry out Likelihood ratio tests

#1: time effect
lrt1 <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,1,0,0,0,0,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt1,p.value=0.01, lfc=0))
[,1]
-1 1559
0  3973
1  1239

#2: Drug1 effect
lrt2 <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,-1,1,0,0,0,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt2,p.value=0.01, lfc=0))

#3: Drug2 effect
lrt <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,1,0,0,0,0,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt,p.value=0.01, lfc=0))

#4: Drug3 effect
lrt <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,1,0,0,0,0,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt,p.value=0.01, lfc=0))

#5: Irrad effect
lrt5 <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,0,0,0,1,-1,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt5,p.value=0.01, lfc=0))
   [,1]
-1    0
0  6771
1     0

#6: congelation effect
lrt <- glmLRT(fit, contrast=c(0,0,0,0,0,0,0,1,0,0,0,0,0))
# Show top ranked sgRNAs
summary(de <- decideTestsDGE(lrt,p.value=0.01, lfc=0))

1)  For the first contrast, I think I can compare T2 vs T1 using contrast=c(0,0,0,0,0,0,0,1,0,0,0,0,0). Am I right?Are the information of paired samples between each T1/T2 used?

2-4) For the study of drug effects, I can test the effect of each drug.
I want to compare each cell line treated by one drug with the same cell line at the same time point with no treatment. For example, for drug 1: Samples 9, 16, 23, 30  vs Samples 8, 15, 22, 29.
I don't know how to precise the "controls" 8, 15, 22, 29 with my model.

5) For the study of irrad effect, I guess my contrast is wrong, but I don't know why.

6) For the study of congelation effect, I don't know how to precise that I want to compare "DrugNI" to (Time=T2 & Drug=none & for the pairs only (without usinfg control cellLines) )

Since it seems difficult to write the contrats, my design is probably incorrect.
Any help would be greatly appreciated.
Thank you in advance.

For your first question, yes, that contrast will test for a non-zero average effect of time across all cell lines. You've blocked on the cell line in the model so the pairing will be automatically considered.

For your second to fourth questions, you should be setting up contrasts like below. This will compare the treated time-2 sample of each cell line to the corresponding untreated time-2 sample, and will average the effects across all cell lines. (The coefficient for each drug represents the log-fold change over untreated time-2.)

con.drug1 <- c(0,0,0,0,0,0,0,0,1,0,0,0,0) # drop DrugDrug1
con.drug2 <- c(0,0,0,0,0,0,0,0,0,1,0,0,0) # drop DrugDrug2
con.drug3 <- c(0,0,0,0,0,0,0,0,0,0,1,0,0) # drop DrugDrug3

I don't know how the Drug3 coefficient managed to come after the irradiated/NI coefficients in your code example, I guess it's something to do with different character orders. Just adjust it appropriately in your setting, or use makeContrasts.

For your fifth question, your irradiation contrast looks fine. Because the coefficient ordering is different for me, it would look like:

# compare DrugIrradiated to DrugNI 
con.irrad <- c(0,0,0,0,0,0,0,0,0,0,0,1,-1)

Finally, if you want to compare non-irradiated time-2 samples to the time-2 untreated sample, you could just do:

con.nonir <- c(0,0,0,0,0,0,0,0,0,0,0,0,1) # Drop DrugNI

I should point out that there's no point blocking on mutation status, because the cell line factor is nested within it. This means that if you added Mut as an additive term, you'd end up with linear dependencies in the design matrix. Instead, what you can do is describe mutation-specific responses to time or drug or irradiation, using something like this:

design <- model.matrix(~0+CellLines+Time:Mut+Drug:Mut)

This will account for the possibility that the control cell lines respond differently from the mutated lines when treated with various things. However, the design matrix will be a lot more complicated, so be warned.