Black hole solutions in general relativity are simple. The frequency spectrum of linear perturbations around these solutions (i.e., the quasinormal modes) is also simple, and therefore it is a prime target for fundamental tests of black hole spacetimes and of the underlying theory of gravity. The following technical calculations must be performed to understand the imprints of any modified gravity theory on the spectrum: 1. Identify a healthy theory; 2. Find black hole solutions within the theory; 3. Compute the equations governing linearized perturbations around the black hole spacetime; 4. Solve these equations to compute the characteristic quasi normal modes. In this work (the first of a series) we assume that the background spacetime has spherical symmetry, that the relevant physics is always close to general relativity, and that there is no coupling between the perturbation equations. Under these assumptions, we provide the general numerical solution to step 4. We provide publicly available data files such that the quasi normal modes of any spherically symmetric spacetime can be computed (in principle) to arbitrary precision once the linearized perturbation equations are known. We show that the isospectrality between the even- and odd-parity quasinormal mode spectra is fragile, and we identify the necessary conditions to preserve it. Finally, we point out that new modes can appear in the spectrum even in setups that are perturbatively close to general relativity.