About the Majorana spinors, think of them as 'real' spinors, while Dirac spinors are 'complex'. For a real scalar field, $\phi = \phi^*$, while for a complex scalar field $\phi$ and $\phi^*$ are independent degrees of freedom. Similarly, for a Dirac spinor $\psi$ and $\psi^*$ are independent degrees of freedom, while for a Majorana spinor $\psi^* = \psi$. The anti-commutation relations for the gamma matrices will be the same as before.
To find the Green's function or propagator of the Majorana spinor $\chi$, you have to - identical to the case for a Dirac spinor - find the inverse of this operator
$$
\bigg[ \gamma^0(i \cancel{\partial}-m_\chi) \bigg]G(x-y) = i\delta(x-y)
$$
which will give you
$$
G(x-y) = \int \frac{d^4 p}{(2\pi)^4} \frac{i(\cancel{p} + m_\chi)\gamma^0}{p^2 - m^2 + i\epsilon} e^{-ip\cdot(x-y)}
$$
Note that $\gamma^0\cancel{p}\gamma^0\cancel{k} = \gamma^0 \gamma^\mu \partial_\mu \gamma^0 \gamma^\nu \partial_\nu = - (\gamma^0)^2 \gamma^\mu \partial_\mu \gamma^\nu \partial_\nu$. You therefore have to put the $\gamma^0$ in the above propagator on the right, to ensure the minus sign will cancel (as every time you pull the $\gamma^0$ 'through' a $\cancel{p}$ it obtains a minus sign).
Because the Majorana spinor is 'real', $\chi$ and $\chi^*$ are the same. $\psi$ is a Dirac spinor and therefore has different chiral components (that are obtained from the projection operators $P_L$ and $P_R$). The interaction terms tell you that this Majorana spinor $\chi$ interacts with the difference chiral parts (the left and right handed components) of the Dirac spinor $\psi$ in different ways
$$
\sqrt{2} \lambda \chi^\dagger\gamma^0 \bigg[ \phi_1 P_L - \phi_2 P_R \bigg] \psi+\mbox{h.c.}
$$
If you want to do Feynman diagrams this means that we will have two distinguish everywhere between the diagrams for $\phi_1$ and $\phi_2$, and $P_L \psi$ and $P_R \psi$, since these different components of the fields will interact in different ways.
For the first interaction term (involving $\chi^\dagger$, $\phi_1$ and $P_L \psi$) the interaction term will indeed be what you said it would be: $\sqrt{2}\lambda\gamma^0$.
Feynman rules by functional derivatives
It is not in general, it just coincides that way for polynomials of fields, without any derivatives or other complications. In truth, one takes functional derivatives until no field is left.
For example, schematically,
$$-i\frac{\delta^4}{\delta \phi^4} \frac{\lambda}{4!} \phi^4 \to -i\lambda$$
which is the entire reason for the factor of $4!$ - it is a convenient convention, but not a necessary coefficient and the physics remains the same without it. More complexly, we could have,
$$-i\frac{\delta^2}{\delta \phi^2} \frac{\delta}{\delta A_\mu}g\phi A^\nu \partial_\nu \phi \to -ig(p_1^\mu + p_2^\mu)$$
which describes a vector coupling to a scalar, where $p_1$ and $p_2$ would be momenta labelling two of the legs attached to the $A\phi\phi$ vertex.
Quadratic terms
The kinetic and mass term is,
$$\mathcal L = \frac12 (\partial \phi)^2 - \frac12 m^2 \phi^2.$$
In Fourier space, $\partial_\mu \phi \to ip_\mu \phi$ and so we have $(\partial \phi)^2$ must go as $p^2 \phi^2$. Interpreting the inverse in Fourier space as the multiplicative inverse, we then have that the propagator goes as,
$$\Delta \sim \frac{1}{p^2+m^2}.$$
Note that for the counterterm Lagrangian (for when you move on to renormalisation), the kinetic and mass counterterms are typically interpreted as interactions, and thus functional derivatives are taken, but no inversion is performed. However, this is a matter of choice, one could absorb coefficients into the propagator instead - they lead to the same physics.
Green's function
Note that a propagator - other than being the inverse of the quadratic terms - can also be interpreted as the Green's function of the equations of motion. This is a function which can be used to solve the equations, via,
$$\phi(x) = \int \mathrm dx' \, G(x,x') f(x')$$
where $(\square + m^2)\phi(x) = f(x)$. Conceptually, in the same way we can think of a function as being built of delta functions, we can think of a solution built up as Green's functions since,
$$(\square + m^2)G(x) \sim \delta^{(n)}(x).$$
Best Answer
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.
I hope this answers your question.