Statistical Thresholds for RNA-Seq Analysis: Considering the Intersection of Treatment Effects

I'm conducting an RNA-seq experiment involving worms exposed to chemicals A, B, and A+B to observe mRNA changes due to exposure. I'm particularly interested in identifying mRNA changes in the A+B treatment compared to treatments A or B alone, which act as controls.

To analyze this, I'm using edgeR, a statistical method, to compare A vs A+B and B vs A+B separately, and then examining the intersection of the results (intersecting A vs A+B and B vs A+B).

In my approach, I'm using an intercept method where I aim to identify transcripts differentially expressed in both A vs A+B and B vs A+B. My question pertains to the justification of using a p-value threshold higher than 0.05 for each separate analysis (A vs A+B and B vs A+B).

To explain further, I'm considering the Multiplication Rule of Probability. If we have event A with a probability of 1/2 and event B with the same probability of 1/2, the probability of both events occurring simultaneously is calculated as (1/2)*(1/2) = 1/4.

I'm curious if a similar principle can be applied here. Specifically, if I start with a p-value of 0.223 for each separate analysis, could I square that value to obtain the p-value for the intercept, resulting in (0.223)^2 = 0.05? However, I'm uncertain about how using the same A+B reads in both comparisons might affect this calculation.

Thank you!

• Such multiplication only holds if the events are independent. Is there any reason why you cannot have the three groups in one GLM and extract the relevant contrasts from it? Apr 3 at 8:30
• Is there a more suitable multiplication method for this scenario? I've attempted to include all three groups in one GLM and extract relevant contrasts from it. However, it seems that this approach also necessitates intersection. Could you provide guidance on how to structure such a GLM and define appropriate contrasts? Apr 3 at 8:47

I agree with Ian Sudbery that this does not sound like a good idea, multiplying probabilities is only valid when they are independent. I see no reason why you could not do this with a regular GLM that contains all groups.

Here's a worked example how you can approach this in edgeR. First, let's generate a hypothetical RNA-seq experiment, no differential expression yet:

ngenes <- 1E3
groups <- c("a", "b", "a+b")
nsamples <- 1E2

ngroups <- length(groups)
group <- factor(rep(groups, each=nsamples), levels = groups)

set.seed(1)
counts <- matrix(
rnbinom(ngroups*nsamples*ngenes, mu=500, size=10),
ncol = ngroups*nsamples, nrow = ngenes
)


On top of this I'll include a number of DE genes, unique to each group. This may not be the most realistic assumption, you might think that the effect of A also appears in A+B for example, but could easily be extended:

## Add random DE to (masked) rows
nde <- length(rows)
dir <- sample(c(-1,1), nde, replace = TRUE) * runif(nde, 1, 3)
counts[rows,] <- matrix(
nrow=nde, byrow=TRUE
)
counts
}

## DE genes in each group
nde <- 10
de_a <- sample(seq_len(ngenes), nde)
de_b <- sample(seq_len(ngenes), nde)
de_ab <- sample(seq_len(ngenes), nde)

de_counts <- add_de(counts, de_a, group=="a") |>


Just to draw your attention to the genes that are DE in A+B but not in the others:

sort(de_ab)
#> 14 125 150 268 325 476 503 872 939 979


With the example data in place, let's analyze using edgeR.

fit <- edgeR::DGEList(de_counts, groups = group) |>
edgeR::calcNormFactors() |>
edgeR::estimateDisp() |>
edgeR::glmFit(design = model.matrix(~ 0 + group))

## Contrast A+B to the mean of A and B
tt <- edgeR::topTags(edgeR::glmLRT(fit, contrast = c(-0.5, -0.5, 1)))


To highlight which top 10 DE genes were identified in this contrast:

sort(as.integer(rownames(tt)))
#> 14 125 150 268 325 476 503 872 939 979

all(as.integer(rownames(tt)) %in% de_ab)
#> TRUE

• The difference between this and the approach I suggest is which genes are being looked for. In this model, you look for genes where the effect of A+B is significantly different to the average of the A and B effects. In my answer you are looking for genes where the the effect of A+B is significantly different to the sum of the A effects and the B effects. Apr 3 at 14:29
• @IanSudbery indeed, this can be obtained easily in either package by reparametrizing to an overall intercept -- otherwise the internal null model will be nonsense -- and using coef=3 (or whichever index is the A+B effect) in edgeR::glmLRT, edgeR::topTags, limma::contrasts.fit or limma::topTable. Apr 3 at 17:07
• Opps, yes I'm so used to considering limma and edgeR interchagably, I missed that I was talking about limma when the OP was talking about edgeR. Apr 3 at 17:30

It is exactly these sorts of considerations that mean that I don't like intersections when analysing this sort of data. The other problem with intersections is they rely on an idea that things that are not significant can be classed as not different. Which clearly isn't the case.

My instinct here would be to use an interaction model. Something like ~0 + A + B + A:B, and comapring it, via edgeR's pseudoLRT to ~ 0 + A + B.

This will identify genes where the effect of combining A and B is different to sum of the effects of A and B independently. However, I'm not quite sure if this is what you are asking for?

• I will try it, thanks. Just to clarify, in ~0 + A + B + A:B, is it the sample A+B , or is it samples A and B ? Also, here ~ 0 + A + B. also, can i do it in edgeR? Apr 3 at 12:48
• No, A and B are indicator variables (seperate columns in your sample table). You specify this formula when you build your model matrix for limma. Apr 3 at 13:56
• My appologies. Yes, this also works in edgeR - I misread edgeR in your OP as limma (same author). Apr 3 at 17:30