How can the Feynman rules be read off the Lagrangian?

I am reading Peskin. In his functional methods chapter he says that

(i) "Once the quadratic terms in the Lagrangian are properly understood" and (ii) "The propagators of the theory are computed" then "the vertices can be read directly from the Lagrangian as the coefficients of the cubic and higher order terms."

What does this mean? In particular: (1) What does it mean that the quadratic terms are properly understood? How can one improperly understand a quadratic? What does this mean?

(2) What does it mean that the vertices can be read directly from the Lagrangian as a coefficient? For example, (2a) how can one determine what the vertex itself looks like? And (2b) in $\phi^4$ theory, the coefficient is $- \lambda/4!$, whilst the Feynman rule for the vertex is $-i\lambda \neq - \lambda/4!$.

Keep in mind that it is nearly impossible to explain how perturbative QFT calculations follow from Lagrangians such that the answer is both relatively short and detailed. So I am going to write an introductory answer. If you want more details on any of its part, you can look up textbooks, or you can let me know in the comments, in which case I will consider updating this answer.

Suppose your model has $n$ quantum fields (they can be organized as Poincare multiplets or all be scalars, for what follows it doesn't matter). The generic expression for the quadratic term in the Lagrangian is thus

$$\mathcal{L}_2 = \frac{1}{2} \left( K_{ab} \partial_{\mu} \phi^{a} \partial^{\mu} \phi^{b} - M_{ab} \phi^{a} \phi^{b} \right).$$

(Actually, if some of the fields have spacetime indices, there could be additional terms like $N_{\alpha a} \psi_{\mu}^{\alpha} \partial^{\mu}\phi^{a}$, but they can be treated in the same matter, so we won't lose generality if we just ignore this issue here.)

First we would like to re-express this Lagrangian, using integration by parts (remember that the Lagrangian is integrated over spacetime to give the action of the system describing its dynamics), as follows:

$$\mathcal{L}_2 = \frac{1}{2} \phi^{a} \hat{Q}_{a b} \phi^{b},$$

where $\hat{Q}$ is the second-order linear differential operator acting on fields. It is called the Euler-Lagrange operator because it generates the classical equations of motion through

$$\hat{Q}_{ab} \phi^{b}_{\text{classical}} = 0.$$

For example, for the multiplet of Klein-Gordon fields it turns out to be

$$\hat{Q}_{ab} = \delta_{ab} \Box + M_{ab},$$

where $M_{ab}$ is called the mass matrix. The basis in which $M_{ab}$ is diagonal is a proper basis for expressing fields associated to elementary particles, the diagonal values being the masses squared of the elementary particles. The d'Alambert operator is $\Box = \partial_{\mu} \partial^{\mu}$.

In the quantum theory we want to calculate the propagator, or the time-ordered product of two field operators:

$$\Delta^{ab} (x, y) = \left< \phi^{a}(x) \phi^{b}(y) \right>.$$

It turns out that the propagator is equal to the Feynman Green's function of the differential operator $\hat{Q}$, which can be derived in the path integral formalism:

$$\hat{Q}_{ab}(x) \Delta^{bc} (x, y) = i \delta_a^c \delta^{(4)} (x - y).$$

This is what is meant by treating the quadratic term in the Lagrangian properly.

At this point it is worth mentioning that sometimes the operator $\hat{Q}_{ab}$ is singular, that is, doesn't have an inverse in the class of functions with radiation boundary conditions. This is because of gauge invariance. The simplest case where this shows up is the free Maxwell Lagrangian.

The modern way of dealing with this is through the formal manipulation with path integrals called the Faddeev-Popov procedure, which introduces additional terms in the Lagrangian (gauge-fixing term and maybe ghost fields). The resulting Lagrangian is still applicable to the same physical model (which is guaranteed by the Faddeev-Popov procedure), but its differential operator is not singular and the propagator can be calculated. This propagator turns out to be unphysical and depends on the unphysical gauge fixing parameter, but when used to calculate S-matrix elements between physical states, the dependence on the unphysical parameter disappears and gauge invariance is restored.

(In fact, gauge invariance is still present in the modified Lagrangian in the form of BRST supersymmetry. Do not confuse it with SUSY.)

Now consider a perturbation of the Lagrangian, i.e. a higher-order term. We deal with such perturbations using, rather unimaginatively, perturbation theory. In the path integral formalism it can be done by Taylor-expanding the exponential of the interaction Lagrangian and making it a part of the correlation functional, keeping the quadratic term as the effective action functional. Then we can apply Wick's theorem (which is only valid for quadratic actions, but hey, that's what is left after we expanded the interaction term) and that would lead us to Feynman rules.

This part is usually the same in all theories, and the final Feynman rules can be easily predicted by simply looking on the structure of the interaction term in the Lagrangian. That is what is meant by "reading off Feynman rules".

For example, consider a single Klein-Gordon field with a 4-th order interaction term

$$\mathcal{L}_4 = - \frac{\lambda}{4!} \phi^4.$$

We would like to Taylor-expand it in any expression for the correlation function of any functional $F$:

$$\left< F[\phi] \right> = \int D\phi \exp \left[ i \int (\mathcal{L}_2 + \mathcal{L}_4) \right] F[\phi] =$$

$$\left< F[\phi] \left( 1 + i \int \mathcal{L}_4 + \frac{i^2}{2} \intop_x \intop_y \mathcal{L}_4 (x) \mathcal{L}_4 (y) + \dots \right) \right>_0,$$

where the subscript $<>_0$ means that we use the free theory action, which is $\int \mathcal{L}_2$, for which the Wick's theorem is applicable.

Each integral in the series above then corresponds to the addition of an interaction vertex to the Feynman diagram. The expression for the vertex is easy to deduce: it is equal to

$$- \frac{i \lambda}{4!} \int d^4 x,$$

with the integral being over the position of the vertex. The factor of $4!$ also appears in the numerator because we have exactly $4!$ ways of contracting 4 operators at the same point with 4 other operators by propagators. Thus the factors nicely cancel out (in fact, it was the reason for choosing $\lambda$ such that $4!$ enters the denominator of $\mathcal{L}_4$ in the first place).

So, we could either keep $4!$ in the vertex expression and consider the $4!$ different contractions which appear after using Wick's theorem unequivalent, or we can consider them equivalent and cancel the factors of $4!$ which is what is usually done in the literature.