Hello,

I have a question regarding differential expression (DE) with limma/voom in an RNA-Seq experiment of ~250 samples. I already posted this question on Biostars and it was suggested that I also post here.

I have 4 different point mutations that are mutually exclusive and I wish to identify DE isoforms specific for each mutation. Additionally, I have reason to believe that the mutations cause similar differential expression for some isoforms and I also wish to identify those.

My questions are:

• What is the best approach to identify the differences and the similarities of those mutations?
• How do I adress the issue of confounding?

Specifically, I have 3, 4, 6 and 17 samples for each mutation respectively, as well as ~200 samples as a control group.

The approach I have come up with so far is to:

1) perform a DE analysis for each mutation comparing it to the control group (which does not include the other mutations), like this:

design.matrix <- model.matrix(~ factor(mut1), data)

design.matrix <- model.matrix(~ factor(mut2), data)... and so on

2) after doing this for each of the 4 mutations, combine them into a single binary variable and perform a DE analysis comparing all of them against the control group:

design.matrix <- model.matrix(~ factor(all.mutations), data

My way of thinking is that since the mutations are mutually exclusive, I can compare each mutation against the control group (which does not include the other mutations) in order to identify DE isoforms specific to each mutation. Afterwards by combing them, I hope to highlight isoforms that are similarly DE for all mutations. If I check the overlap of the last analysis with the preceding 4 I should be able to identify at least some isoforms that are affected in a common way. Is it maybe sufficient to just compare the overlap of the first 4 analyses without the 2nd step?

As to the 2nd question. Is there a rule of thumb as to how many confounders I can add in a limma analysis while avoiding overfitting? Since in a differential expression analysis we don't really have "events" I am unsure how to determine the number of confounders I can adjust for. Especially for the mutation with the smallest subset (only 3 mutations) I am unsure if the relatively large control group of ~200 samples permitts me to adjust for multiple confounders.

Any input is welcome, thanks in advance!

Stefan

For your first question: your suggested approach is inefficient from both a statistical and computational perspective. Statistically inefficient as it does not use information across the entire dataset to model variability; computationally inefficient because you're repeating the entire analysis for each comparison. Rather, you should do something like:

all.mutations <- rep(c("A", "B", "C", "D", "con"), c(3, 4, 6, 17, 200))
all.mutations <- relevel(factor(all.mutations), "con") # set control as baseline
design <- model.matrix(~all.mutations)

... and fit one model for the entire dataset. Each non-intercept coefficient here will represent the log-fold change of a mutation against control, so it's easy to determine the effect of each mutation relative to the control. You can also quantify the differences between mutations by directly testing for DE between the corresponding groups.

Identification of genes that respond "similarly" across mutations is harder, depending on how you define "similar". If you only care about the fact that they are DE (in the same direction), then you can just intersect the DE lists from the relevant mutation-control comparisons. If you want also same magnitude of log-fold change across mutations... this is difficult as you're asking for genes where the null hypothesis (of no difference between mutations) is true. Standard hypothesis tests aren't built to answer that question. The best you can do is to check that the genes are not significant in the mutation-mutation comparisons.

For your second question: I assume you're referring to blocking factors rather than confounding factors, as your model becomes overfitted by definition once you have one factor that is confounded with the mutation status. As long as you have enough residual d.f. for variance estimation, limma will still work in terms of controlling type I error. Of course, unnecessary increases to the complexity of the model will also reduce power as the estimates of the model parameters become more uncertain.

I would not worry about the number of additional factors, but consider their impact on the downstream analysis. This involves some thinking about the nature of the biological system, and how certain factors might affect the differential expression results. For example, I always block on processing batch in my DE analyses, as I know that batch effects exist in my datasets. Yes, this means you'll have to think about the experiment (the horror!), which is always difficult - but blindly using automated model selection is arguably worse, see some comments here. Unfortunately, the only way to make analysis easy is to have a good experimental design in the first place.