We numerically study the evolution of the vibrational density of states $D(omega)$ of zero-temperature glasses when their kinetic stability is varied over an extremely broad range, ranging from poorly annealed glasses obtained by instantaneous quenches from above the onset temperature, to ultrastable glasses obtained by quenching systems thermalised below the experimental glass temperature. The low-frequency part of the density of states splits between extended and quasi-localized modes. Extended modes exhibit a boson peak crossing over to Debye behaviour ($D(omega) sim omega^2$) at low-frequency, with a strong correlation between the two regimes. Quasi-localized modes instead obey $D(omega) sim omega^4$, irrespective of the glass stability. However, the prefactor of this quartic law becomes smaller in more stable glasses, and the corresponding modes become more localized and sparser. Our work is the first numerical observation of quasi-localized modes in a regime relevant to experiments, and it establishes a direct connection between glass stability and soft vibrational motion in amorphous solids.