Is it possible to run DESeq2 with gene-dependent binary variable such as membership in a pathway or transcription factor association? If not, in principle is it possible to modify functions fitNbinomGLMs, fitBeta, etc. so that it becomes possible in a meaningful way, or I miss some issue that doesn't allow this?
Using gene-dependent covariates can make sense - not for DESeq2 per se, but for the subsequent multiple testing analysis. Have a look at the IHW package and its associated paper. In this way, you can take advantage of the fact that the prior probability of being differentially expressed is quite different between the groups, or if you have reason to believe that the tests' power is different.
In fact, this applies not only to multiple testing analyses following DESeq2, but following any method that conducts multiple tests.