Electron-phonon interaction plays an important role in metals and can lead to superconductivity and other instabilities. Previous theoretical studies on superconductivity are largely based on the Migdal-Eliashberg theory, which neglects all the vertex corrections to electron-phonon coupling and breaks down in many unconventional superconductors. Here, we go beyond the Migdal-Eliashberg approximation and develop a nonperturbative Dyson-Schwinger equation approach to deal with the superconducting transition. Remarkably, we take into account all the vertex corrections by solving two coupled Ward-Takahashi identities derived from two global U(1) symmetries and rigorously prove that the fully renormalized electron propagator satisfies a self-closed integral equation that is directly amenable to numerical computations. Our approach works equally well in the weak and strong coupling regimes and provides an efficient method to determine superconducting $T_c$ and other quantities. As an application, our approach is used to investigate the high-$T_c$ superconductivity in one-unit-cell FeSe/SrTiO$_3$.