I am following a course on QFT that is based upon canonical quantization and not path integrals.
When calculating scattering amplitudes we compute, at the relevant order, the corresponding matrix element, which might look like
$$\int \textrm{d}^4x\textrm{d}^4y\,\langle 0|\hat{a}_p\hat{a}_q T\{\bar{\psi}_x\phi_x\psi_x\bar{\psi}_y\phi_y\psi_y\}\hat{a}^\dagger_r\hat{a}^\dagger_s|0\rangle .$$
In order to compute this matrix element, we simply replace the time ordered product by the normal order of contractions, and then we see we need to fully contract the fields with the external states, as is explained in Peskin and Schroeder and in Computing S-Matrix Elements from Feynman Diagrams.
We then come up with Feynman Rules, define propagators as vacuum expectation values and retrieve the expression for the amplitude.
However, I have just read that both propagators and Feynman Rules can be recovered from the Lagrangian of the theory. How is this related to the usual definition of propagators $$\langle 0|T\{\phi_x\phi_y\}|0\rangle=D(x-y)$$ and the Feynman rules derived from the Wick's theorem?
Also, please, could you provide me with a good reference to derive Feynman Rules from Wick's Theorem? What I have found so far is just a list of rules to apply, but I would like to understand where they come from.
Edit: Just to clarify the question. I am familiar with the definition of propagator $\langle 0|T\{\phi_x\phi_y\}|0\rangle=D(x-y)$. How is this related to the fact that propagator can also be retrieved from second order terms in the lagrangian? View for instance https://cds.cern.ch/record/319569/files/AT00000309.pdf.
Second part of the question is: since Feynman Rules are just a way to calculate Wick contrations in the reasoning I exposed in the first part of my question, I do not see why they can also be derived from the Lagrangian, as explained in the same link. In P&S, Feynman Rules are obtained as a way to compute the matrix element of the fully contracted fields and initial/final states.