Hi,

We are currently analysing our RNA sequencing data and have decided on EdgeR for calling differentially expressed genes. We have two questions about the analysis. One concerns correcting for the effect of an experimental procedure, the other about the test used to call differentially expressed genes.

We have done RNA sequencing on samples that were either directed through a flowcytometer ('sorted') or not ('unsorted'). The 'sorted' samples were sorted for six different cell types. We know that in EdgeR, we can make a design matrix and thus determine differential gene expression per cell type (and per the unsorted sample). In addition, we could add the replicate number as a variable to take into account (since we always collected a mock and a treatment sample for each cell type at the same time).

We are now wondering whether we could also correct for the 'sorting' effect by adding the extra variable 'sorting' to the design matrix. The other solution to correct for the sorting effect that we were thinking of, is looking for genes that relatively contribute a lot to the clustering by sorted vs unsorted and removing those genes from the analysis. We hope we explained this clear enough. Which of these do you think is the better option?

Second, we were wondering whether we should test with GLM or QL. For us, it seems like the user's guide says the QL method is better in general, (page 21 'the QL F-test is preferred as it reflects the uncertainty in estimating the dispersion for each gene'), but that made us wonder why the GLM method is out there in the first place.

Any help would be greatly appreciated.

Eline Verbon and Ronnie de Jonge

Let's set up an example below. I'm assuming that you have three batches of mock/treatment replicates, where all samples in a single batch (across all tissue types) were generated at the same time, from the same plant, etc.

tissue <- rep(c("root", "root", "type1", "type2", "type3", "type4", "type5"), each=6)
sorting <- rep(c("unsorted", rep("sorted", 6)), each=6)
exposure <- rep(rep(c("exposed", "unexposed"), each=3), 7)
batch <- rep(1:3, 14)

I would suggest a design matrix like this:

group <- factor(paste0(tissue, ".", sorting, ".", exposure))
batch <- factor(batch)
design <- model.matrix(~0 + group + batch)

The inclusion of the sorting factor in group will allow for differences between the whole root sorted and unsorted samples. You can determine the sorting effect by testing for DE between these two groups. However, you shouldn't compare between the sorted cell types and the whole root unsorted samples, because there's no way to disentangle cell type-specific changes from the sorting effect.

Obviously, you can also use this design to test for the treatment effect in each tissue, or for DE between cell types. I don't see a need to correct for the sorting effect here, because you're comparing sorted groups so any sorting effect should cancel out. If it doesn't for some reason, then you can just filter out those DE genes that were also found in the whole root sorted vs. unsorted comparison.

As for the QL vs GLM question, there are several reasons that I can think of/remember:

1. The QL framework is an extension of the GLM machinery, so it makes sense to keep the underlying functions as well as the higher-level extensions. For example, glmQLFTest depends on glmLRT.
2. The LRT still works for routine DE analyses, albeit not as well as the QL methods.
3. They are retained for compatibility with older scripts and packages.
4. In some cases, the LRT actually provides better performance, e.g., for single-cell data where there's plenty of samples (so loss of type I error control from uncertainty in the dispersions is less of an issue) but low, highly overdispersed counts (which compromises the accuracy of the saddlepoint approximation required by the QL methods).

#make design matrix
tissue <- s2c$cell_type #s2c is a table with all the information about the samples, which the order as the order of my sequencing files
sorting <- s2c$sortedness
exposure <- s2c$treatment
replicate <- s2c$replicate

group <- factor(paste0(tissue, ".", sorting, ".", exposure))
replicate <- factor(replicate)
design <- model.matrix(~0 + group + replicate)

#for design matrix see next post: not enough space in 1 comment to add it.

y <- estimateDisp(y,design)
fit <- glmQLFit(y,design)

#perform tests

for(x in 1:length(unique(s2c$cell_type))) {
  #get and save fold change and FDR for each gene in each sample
  coeffi <- c(rep(0,17)) #to indicate which which groups (=columns) to compare in the design matrix)
  #not sure what to do with the replicate effects 
  
  #find the samples in the info table in which this cell_type is present
  samplesWithThisCellTypeWCS <- which(grepl(unique(s2c$cell_type)[x], s2c$cell_type) & grepl("WCS", s2c$treatment))
  samplesWithThisCellTypeNB <- which(grepl(unique(s2c$cell_type)[x], s2c$cell_type) & grepl("NB", s2c$treatment))
  
  #find which of the first 14 columns of the design matrix contains a '1' in the rows considering to the correct samples (the first occurrence is enough) --> group containing that cell type + treatment
  columnCellTypeWCS <- which(design[samplesWithThisCellTypeWCS[1],]==1)
  columnCellTypeNB <- which(design[samplesWithThisCellTypeNB[1],]==1)

  #make the coeffi 1 for the treatment column, -1 for the mock column 
  coeffi <- replace(coeffi, columnCellTypeWCS, 1)
  coeffi <- replace(coeffi, columnCellTypeNB, -1)
  
  ql <- glmQLFTest(fit, contrast=coeffi)
  results <- topTags(ql, n=Inf)

   #there is more code here to save all results and the up- and down-regulated genes with a certain fold change and FDR

}

Dear Aaron Lun (and others reading along),

We have continued working on this since our discussion here on the web. Here I first explain how we used your method and run into problems when trying to validate it. Next, I explain how I tried to correct for sorting myself and the problems I run into then. We are probably not understanding some details of your (simpler) method and don't need our own method. However, I included it so show our reasoning, which might help you understand where we go wrong trying to use your method. 

It became quite a bit of text, I hope you will understand what I say.

To correct for our sorting effect we followed your suggestion to make a design matrix.

Afterwards, we make the fit object and calculate the genes changed due to our treatment. (As you said, I could make the contrasts by using makeContrasts, but my, indeed more complex, method only works, so for now I have kept it like this. I will adjust later.)

The code works (no bugs and we end up getting fold changes per gene for each of our cell types and the controls), but we're still struggling to find out how to determine whether the correction works and adjusts the values in the sorted samples as needed. In other words, how we can convince ourselves and others that the sorting effect is adjusted for and our results show the effect of our treatment (almost) only.

We figured we could plot our fold changes of all the genes in the sorted control on the x-axis and the unsorted control on the y-axis (because as you mentioned, we cannot compare with the cell types directly, too many confounding variables). We did this both with correcting for sorting and without (in the latter case, removing 'sorting' from 'group'). We expected the plot without the correction to have points that fell less close to the identity line then the points in the plot with correction. However, we got two exactly the same graphs (with the same R squared). To us, this seems to indicate that the sorting correction did not do anything. 

What do you think? Are we missing something here?

Ps. We collected 1 replicate of 2 cell types or controls per day, with and without bacteria. So the experiment was performed over many different days and we do not have batches we can adjust for (or very many of them). Hope that makes the batch-thing clear. I am not sure what you mean with a 'target file', but if you explain I would be happy to include it.

To correct for the sorting effect myself, I

1) took the counts per gene of the DGELists of the untreated sorted and unsorted control samples. I calculated the fold change due to sorting with EdgeR and saved this in a table.
2) divided the counts of all the sorted treated samples by the fold change calculated from the treated control samples
3) repeat step 1 and 2 for the treated control samples

The next step is to make a new DGEList to perform EdgeR on again (and normalize). So, something like this (this is the code I used to generate the DGEList all the way at the beginning of my script, which I would have to adapt):

However, here I run into problems because I started with a DGEList to calculate the fold changes due to sorting with EdgeR. But that DGEList is already filtered, which leads to problems when trying to normalize the new DGEList. I could remove the filtering before correcting for the sorting effect, but then the changes due to sorting are not calculated as they should be.

(sorry for the long messages, tried to make it as clear as possible)