You can start by doing PCA, and seeing if you can determine what variable drives the first 2 or 3 principal components. Your ability to account for more than that is going to depend on how many samples you have. If you have only 8, you probably won't be able to model anything more than disease status. With 100, you could probably include more things, like age, sex, other lab results. But of those other things aren't visible in the PCA, they won't matter much.

Its also difficult to decide what covariates to include and which not to. In general correcting uneccesarily for covariates costs you power becuase the model has more parameters to estimate. One might correct for covariates for two reasons:

1) Avoid false positives due to changes actaully being associated with a covariate, rather than the condition of interest 2) Avoid false negatives due to the uncorrected covariates increasing the heterogeneity.

It is also unlikely that uncorrected covariates will cause false positives as long as the ages/co-morbilities are balanced between the disease conditions. If they are not, then you probably need to talk to an epidemiologist because I forsee that there are all sorts of issues to be addressed in such a design. In general, I don't think that normal frameworks we use for RNAseq DE are well suited for causal inference. That is, they are good as saying "genes different between these two groups of samples that are from different conditions", but less good at saying "gene differ between these two groups of samples BECAUSE they are from different conditions". My understanding is that they were designed with experimental data (where most confounders are balanced between conditions via experimental design), rather than observational data (where they are not) in mind. Related to that is the fact that the models we use are generally "fixed-effects" models, but I think (although I'm sort of straying beyond my statistical comfort zone) that things observed after the fact are "random-effects". Limma can handle random effects to some extent via its duplicateCorrelation function, but using that in this case would be beyond me.

In the case of avoiding false negatives, you'd want to know that the covariates in question really were inducing extra hetrogenity before correcting for them, because as I said before, correcting for unneccessary covariates reduces power, while correcting for the right ones increases power. You might like to look at something like a PCA, or and MDS or heirarchical clustering of the data to see if those covariates seem to be driving sample sepeartion. If they are, thats a good arguement for correcting for them. If they are not, it might reduce rather than increase power.