Examining possible extraneous variables in RNA-sequencing data

[SuzuBell]

I have an RNA-sequencing dataset with 48 samples. I have a list of several possible extraneous variable values for each of the 48 samples. Some of these extraneous variables are categorical (lane1 versus lane2, performed on day 1 versus performed on day 2, performed in room 1 versus performed in room 2) and some are quantitative (ng/UL of RNA submitted for RNA-seq, RIN, % values, etc).

The dataset contains two factors (one with 4 levels and one with 2 levels). Hence, there are 8 treatment groups, and there are 6 samples per treatment group. I did all 28 pairwise DEG analyses between the 8 treatment groups using edgeR, DESeq2, and limmaVoom, and all three methods revealed similar and disappointing results. Some treatment pairs had thousands of DEGs, and others had zero DEGs (even if that particular pair should have many in theory).

I am trying to determine if these possible extraneous variables (both categorical and quantitative) may be causing issues. The MDS plots are not too good. Is there a recommended/reliable method for testing and/or accounting for extraneous variables (both categorical and quantitative) in RNA-sequencing DEG analysis? I am most familiar with edgeR, DESeq2, and limmaVoom, but am open to other software.

Thank you for any advice!