How are Feynman rules derived, in general? A few questions about the derivation of Feynman rules in specific instances have already been asked here 
(e.g.
Sign of Feynman rules with derivative couplings, Feynman rules for coupled systems, How can we derive the Feynman rule for the ordinary QED 3-vertex?, Recipe for computing vertex factors in Feynman diagrams).
However, a more general discussion seems to be missing, hence this question.


*

*How are Feynman rules for a generic theory derived from the Lagrangian density? How are the various methodologies (e.g. second quantization, functional quantization) related to each other? When (or if) is one preferable to the others?

*When and why are additional complications (like Faddeev-Popov procedures) involved in the derivation of the rules?
 A: *

*The simplest way is to look at arbitrary process amplitude (S-matrix), to expand it in a series of some constant, then - to use Wick theorem and, finally, to get that n-th amplitude consists of sum over multiplications of all possible numbers of propagators and normal ordering field operators. 
Sometimes it's convenient to use nonperturbative methods. Formally it means that we change interaction representation to Heisenberg representation and introduce n-point functions. It can be shown that n-point function in Heisenberg picture
$$
\tag 1 \langle | \hat{T}\left(\hat{\varphi}^{H}_{1}...\hat{\varphi}^{H}_{n}\right)|\rangle 
$$
consists all Feynman diagramms with $n$ external lines for a given theory. By using this idea you may "prove" Feynman rules even in nonperturbative formalism, but with modification of fields, mass and bounding constants (LSZ-theorem).
The connection between operator formalism and path integration formalism is in the different methods of calculating of n-point functions. For example, in terms of path integration $(1)$ looks like
$$
\langle | \hat{T}\left(\hat{\varphi}^{H}_{1}...\hat{\varphi}^{H}_{n} \right)|\rangle = \left[\frac{\delta }{\delta J_{1}}...\frac{\delta }{\delta J_{n}}\int D\varphi_{1}...D\varphi_{n}e^{i S[\varphi] + i\int d^{4}x J^{i}\varphi_{i}} \right]_{J = 0}.
$$
Path integration formalism is very convenient because it includes full action of theory (not only interaction summands), so contains all information about symmetries of theory. It helps to derive relations for vertex functions (like Slavnov-Taylor identities, Ward identities etc.). You may also "derive" Feynman rules by using only path integration formalism. 

*Complications of description mostly can be appeared only when you use nonperturbative description for processes. There are three examples. 


*

*When using standard perturbation theory (expanding S-operator in a series) you automatically don't take into account bounded states (like pion) and topological configurations (like instantons) which also appear in interaction theory. But non-perturbative methods allow to take into account these complications. 

*When you assume some theory with connections of the first kind between canonical coordinates and momentum (like gauge theories) you need to eliminate this problem, because it doesn't permit to determine independent set of dynamical variables (i.e. to quantize the theory). But when you use some gauge condition, it helps to eliminate the problem (in the language of path integration it means that you reduce the number of integrations in the fields lying on the gauge orbit), but the payment for this is appearance of fictive field (ghosts).

*Sometimes some symmetries are broken in nonperturbative description. For example, there is well-known the violation of CP-symmetry (when you assume topological constructions in path integral formalism) and chiral symmetry.
