The gravitational wave signature from core-collapse supernovae (CCSNe) is dominated by quadrupolar oscillation modes of the newly born proto-neutron star (PNS), and could be detectable at galactic distances. We have developed a framework for computing the normal oscillation modes of a PNS in general relativity, including, for the first time, the presence of an accretion flow and a surrounding stalled accretion shock. These new ingredients are key to understand PNS oscillation modes, in particular those related to the standing-accretion-shock instability (SASI). Their incorporation is an important step towards accurate PNS asteroseismology. For this purpose, we perform linear and adiabatic perturbations of a spherically symmetric background, in the relativistic Cowling approximation, and cast the resulting equations as an eigenvalue problem. We discretize the eigenvalue problem using collocation Chebyshev spectral methods, which is then solved by means of standard and efficient linear algebra methods. We impose boundary conditions at the accretion shock compatible with the Rankine-Hugoniot conditions. We present several numerical examples to assess the accuracy and convergence of the numerical code, as well as to understand the effect of an accretion flow on the oscillation modes, as a stepping stone towards a complete analysis of the CCSNe case.