Let's take for granted the Gaussian integration formula, which holds for both bosonic and fermionic integrals, if they are properly interpreted:

__Theoreom (Gauss, Wick):__ Let $X$ be a vector space with a chosen volume form ${\rm d}x$, $f: X \to \mathbb C$ a polynomial, and $a: X^{\vee 2} \to \mathbb C$ a symmetric bilinear form with inverse $a^{-1}\in X^{\vee 2}$ such that the Gaussian measure $\exp(-\frac12 a\cdot x^{\otimes 2}){\rm d}x$ is defined.  (So for example $X$ can be bosonic finite-dimensional and $a$ can have positive-definite real part; or $a$ can be invertible pure-imaginary and all integrals can be taken as conditionally convergent; or $X$ can have even-dimensional fermionic parts and the integral can be defined a la Berezin.)  Then we have
$$ \int_X \sum f^{(n)} \cdot \frac{x^{\otimes n}}{n!} \exp \left(-\frac12 a\cdot x^{\otimes 2}\right){\rm d}x = \sqrt{\det(2\pi a)} \sum f^{(2k)} \cdot \frac{(a^{-1})^{\otimes 2k}}{2^kk!} $$
Or, anyway, when $X$ is finite-dimensional Bosonic this is correct.  In the fermionic case, it's off by some $\sqrt{-2\pi}$s, and $\det$ is Berezin's superdeterminant.  Here $f^{(n)} : X^{\vee n} \to \mathbb C$ is the $n$th Taylor coefficient of $f$ at $0$; it's a symmetric tensor.  In the fermionic case, some care must be taken with words like "symmetric", but I will ignore this subtlety.

__Proof:__ Integrate by parts.  $\Box$

Now the trick is to interpret the RHS combinatorially: you should draw each summand as a graph with one vertex, labeled $f^{(2k)}$, and $k$ self-loops, labeled $a^{-1}$; then if you count automorphisms correctly, the denominator of the summand is the number of automorphisms of the graph, and the numerator is the "evaluation" of the graph as a picture of tensor contractions.  You can also draw the left hand summands combinatorially: a vertex with $n$ incoming strands corresponds to $f^{(n)}$, and the $n!$ counts the symmetries.

What would have happened if you had not just a single polynomial but a product?  You can draw $\frac1{m!} f^{(m)} \otimes \frac1{n!} g^{(n)}$ as two vertices, one labeled $f$ with $m$ incoming strands and the other labeled $g$ with $n$ incoming strands, and the symmetries are correct, and using $\frac1{m!} f^{(m)} \otimes \frac1{n!} g^{(n)} = \binom{m+n}{m,n} (f^{(n)}\otimes g^{(m)})$ and the above theorem, the RHS now can be taken as a sum over all graphs (possibly disconnected) with two vertices, one labeled $f$ and the other labeled $g$, where the edges are still valued $a^{-1}$ and a labeled graph is weighted by its automorphism group.

Ok, so now let's try to calculate asymptotics of integrals, and I will henceforth ignore the "determinant" prefactors.  My domain of integration will always be a "small" neighborhood of $0$.  I'm interested in the $\hbar \to 0$ asymptotics of
$$ \int_X \exp \frac1\hbar\left( -  a\cdot \frac{x^{\otimes 2}}2 + \sum_{n\geq 3} b^{(n)}\cdot \frac{x^{\otimes n}}{n!} \right) {\rm d}x $$
First rescale $x \mapsto \sqrt\hbar x$; this just rescales the overall integral by $\hbar^{\dim X/2}$, and I'm dropping those terms.  Then the integral is $\exp\bigl( -  a\cdot \frac{x^{\otimes 2}}2 + O(\hbar) \bigr)$, so keep the $O(1)$ part in the exponent but expand the $O(\hbar)$ part out:
$$ = \int_X {\rm d}x \exp \left(-\frac12 a\cdot x^{\otimes 2}\right) \times \sum_{m\geq 0} \frac1{m!} \left( \sum_{n\geq 3} b^{(n)}\cdot \frac{x^{\otimes n}}{n!} \right)^m$$
Expanding the sum further, the $b$s still look like vertices with $n\geq 3$ incoming strands (and $n!$ symmetries), but now I get to have $m$ many of them, weighted by the $m!$ symmetries from permuting the vertices.  So the sum is over all collections of trivalent-and-higher vertices, counted with symmetry, and an $n$-valent vertex is valued $\hbar^{\frac n 2 - 1}b^{(n)}.

Then we integrate by connecting them up.  All together, we get:
$$ = \sum_{\text{graphs } \Gamma } \frac{ \hbar^{\#} \operatorname{ev} (\Gamma) }{ \lvert \operatorname{Aut} \Gamma \rvert } $$
Graphs can be disconnected, empty, have parallel edges and self loops, etc.  We compute $\operatorname{ev}(\Gamma)$ by assigning $a^{-1}$ to each edge, $b^{(n)}$ to an $n$-valent vertex, and contracting tensors.  The power on $\hbar$ is $-1$ for each vertex, and $\frac12$ for each half-edge (each part of an edge that arrives at a vertex), i.e. it's $-1$ for each vertex, $+1$ for each edge, i.e. it's the negative of the Euler characteristic of the graph.

Finally, let me get to the version from Getzler-Kapranov.  Above, I used the fact that if $\star$ is a sum of connected things (counted with symmetry), then $\exp(\star)$ is a sum (counted with symmetry) over all possible collections of disjoint copies of $\star$.  We've ended up with a sum-with-symmetry over disjoint things.  Taking $\log$ gives the sum of connected things.  For a connect graph, the negative Euler characteristic is precisely one less than the first Betti number.

__Exercise:__ Redo the above calculations with $O(\hbar)$ corrections to the exponent of the initial integral, to end up with the precise Getzler-Kapranov result.

-----

Finally, I should say one more thing about Feynman diagrams.  Feynman and Dyson disagreed about the meaning of Feynman's diagrams: Feynman thought of them as pictures of particles interacting, and Dyson thought of them they way I do above: as graphs computing the asymptotics of integrals.

The point is the following.  The most important integrals like those above that physicists care about come from particularly nice quantum field theories, where $X$ is a space of sections of some vector bundle (with some extra structure) on a Riemannian manifold, and $a$ is the Laplacian, and the $b^{(n)}$ are all "local".  In this situation, $a^{-1}$ computes the "heat flow" for the Riemannian manifold, which is nothing more nor less than an integral over all paths connecting two endpoints of some "amplitude for a free particle to travel along this path".  Then $\operatorname{ev}(\Gamma)$ can be interpreted as an integral over all embeddings of $\Gamma$ into your manifold of: the amplitude for your particle to travel along the edges, times an amplitude for an "interaction" at the vertices.

One good reference is the first chapter or two of Costello's recent book on QFT.