We apply a novel approach based on physics-informed neural networks to the computation of quasinormal modes of black hole solutions in modified gravity. In particular, we focus on the case of Einstein-scalar-Gauss-Bonnet theory, with several choices of the coupling function between the scalar field and the Gauss-Bonnet invariant. This type of calculation introduces a number of challenges with respect to the case of General Relativity, mainly due to the extra complexity of the perturbation equations and to the fact that the background solution is known only numerically. The solution of these perturbation equations typically requires sophisticated numerical techniques that are not easy to develop in computational codes. We show that physics-informed neural networks have an accuracy which is comparable to traditional numerical methods in the case of numerical backgrounds, while being very simple to implement. Additionally, the use of GPU parallelization is straightforward thanks to the use of standard machine learning environments.