To my understanding, unlike free fields, interacting fields cannot be expanded in terms of Fourier modes, with the Fourier coefficients representing creation and annihilation operators. Then is it possible to quantize an interacting field theory like the self-interacting $\phi^4$ theory or the QED? If yes, what is the prescription and how does the quanta of the interacting $\phi^4$-theory be related to the quanta of the free scalar field?
There is a well defined formalism for this, which is called perturbative quantum field theory. It follows the spirit of doing perturbation theory in quantum mechanics, only here the perturbative expansion is done about some adequately fixed field configuration.
The problem with doing a mode-expansion of fields for any general Lagrangian is a difficult one as we don't get decoupled equations in the momenta space for the fields, so we can't really do mode-expansions for the fields easily. Instead, the thing of paramount importance here is the time ordered correlation function of fields in the interaction picture (the interaction picture is in a sense a combination of the Schrodinger and Heisenberg pictures). It can then be shown that any such time ordered correlation function of the fields in interaction picture can be represented as an series expansion in orders of the coupling constant of the theory (for $\lambda \phi^4 $ theory the coupling constant is $\lambda$), which can be represented by disconnected Feynman diagrams. We can then construct the S-Matrix or calculate decay rates easily enough once we know how to compute these terms, for which the Feynman diagram approach is the most convenient one.