When aiming to identify differentially expressed genes, thousands of simultaneous hypothesis tests are performed, which could be biased by the presence of unmeasured confounders. In the context of linear models, surrogate variable models and related approaches have been developed to control for the effect of confounding factors with considerable success. However, in recent years, differential expression testing has been dramatically expanded to include a variety of genomic readouts for which the linear model rarely holds. A Poisson, negative binomial or Bernoulli model is likely more appropriate. Inspired by this advancement we develop a solution for multivariate generalized linear models in the presence of arbitrary confounding effects. We establish consistency and asymptotic normality of our proposed test statistic. Numerical experiments demonstrate that the proposed method controls the false discovery rate and is more powerful than alternative methods. By comparing single-cell RNA- seq counts from Lupus and control samples, we demonstrate the suitability of adjusting confounding effects when significant covariates are absent from the model. If time permits we will explore how proximal causing learning methods could provide an alternative approach for removing the effects of unmeasured confounders.