Hi there, I have been using edgeR for quite some time now, and now I am facing with a question that applies to experimental designs that are not the "traditional" control vs treated design, but to cases where multiple treatments are being tested.
In my experimental design, I have 3 different micro plastic treatments (whose effects were assessed: MIX1, MIX2 and MIX3) and one control (i.e. CTL). For each treatment I have 4 independent biological replicates (n=16 total RNASeq libraries). The goal is to assess each of the micro plastic against the control by doing 3 pairwise comparisons.
My question applies to the initial filtering steps. I use filterByExpr() on the raw count table and I provide the group variable (which I input as a factor) to define all the 4 groups, then I set min.count = 5, min.prop = 0.8.

All works like a charm, I then run normalisation and DEA by running all the pairwise comparisons that I mentioned above. However my doubts arise when looking at the results of some pariwise analysis. In example, in the comparison CTL vs MIX1 I get, among the other differentially expressed genes, this one (I am reporting the statistics):

    logFC   logCPM  LR  PValue  FDR
transcript_172821   22.80989877 1.655460041 37.96739989 7.19366E-10 2.91137E-06

This gene has a very low FDR and high logFC. However, when looking back at the original count table I see that the counts were as this:

    MIX1_r1 MIX1_r2 MIX1_r3 MIX1_r4 CTL_r1  CTL_r2  CTL_r3  CTL_r4 MIX2_r1  MIX2_r2 MIX2_r3 MIX2_r4 MIX3_r1 MIX3_r2 MIX3_r3 MIX3_r4
transcript_172821   0   0   0   0   112.941 0   0   0   0   45.9024 6   56.1972 0   0   0   0

I understand that the gene has passed the filtering criteria because it was expressed in one CTL sample and in 3 other samples belonging to MIX2, however I am wondering on the biological meaning of this result: Should I not consider this gene as differentially expressed in the comparison CTL vs MIX1 since it only had reads in 1 sample of the 8 samples (i.e. 4 CTL and 4 MIX1) being compared? So I am wondering if I should apply the filterByExpr separately each time on a count table that only includes the two treatments that I am going to include in the pairwise comparison, rather than applying it on the whole count table (i.e. with all 4 treatments together), to avoid situations like the one above?

I'd like to know your feelings about
Looking forward to hearing from you people

Thanks. Luca

We can't offer biological interpretations because it depends on your experiment and the scientific context.

You could subset the data for each pairwise comparison, but I don't recommend that and I don't do it myself. I would instead revisit how you are doing the differential expression analysis and why your data is so over-dispersed:

1. Leave the filterByExpr settings at the default values. This particular transcript would already have been filtered had you left min.count at the default value of 10. And changing min.prop can't have any effect for your data.
2. You appear to be doing analysis of transcripts (isoforms) rather than of genes. edgeR is not designed to analyse transcript level data unless you account for read assignment ambiguity as described by Baldoni et al 2023.
3. We recommend a quasi-likelihood analysis with glmQLFit etc but you are using the older glmLRT pipeline.
4. If you use glmLRT, then it is essential that you use tagwise dispersion values rather than trended or common. You might already be doing that but, given that this gene is so significant, perhaps not. The negative binomial dispersion for this transcript will be very large.