Exactly how you count symmetries depends on your normalizations, and in particular on whether you use divided (i.e. "exponential") power series or ordinary power series. I believe strongly that divided powers are the way to go. To establish notation, I will first review the preliminaries, that I'm sure you already know.
So let me consider Dyson series for integrals of the form:
$$\int_{x\in X} \exp\left( \sum_{n\geq 2} \frac{c_n x^n}{n!} \right) {\rm d}x$$
In actual physics examples, $X$ is an infinite-dimensional space and the measure ${\rm d}x$ does not exist, but no matter. Each coefficient $c_n$ is a symmetric tensor $X^{\otimes n} \to \mathbb R$ or $\mathbb C$, and we suppose that $c_2$ is invertible as a map $X \to X^*$ — in infinite-dimensional settings, this is a very problematic supposition, and leads to questions of renormalization, which I will not address here. (Incidentally, in all actual physical quantum field theories, $c_2$ is not a priori invertible, because of gauge symmetry; so I am assuming that you have fixed that however suits your fancy.) Finally, Dyson series / Feynman diagrams are by definition perturbative, meaning that you need some perturbation parameter. There are various choices for how to do this; ultimately what's important is that ratios $\lvert c_n\rvert / \lvert c_2 \rvert^{n/2}$ are infinitesimal for $n \geq 3$. In your example, $c_2$ is normalized to unity, and $c_4 = \lambda \ll 1$. Another option is to set all coefficients $\lvert c_n\rvert \sim \hbar^{-1}$, where $\hbar \ll 1$.
In any case, the Dyson series is correctly calculated by expanding
$$ \int_{x\in X} \exp\left( \sum_{n\geq 2} \frac{c_n x^n}{n!} \right) {\rm d}x = \int_{x\in X} \exp \left( \frac{c_2x^2}2\right) \sum_{m=0}^\infty \left( \sum_{n\geq 3} \frac{c_n x^n}{n!} \right)^m {\rm d}x, $$
using "Wick's theorem" for $b_\ell$ a totally symmetric tensor $X^{\otimes \ell} \to \mathbb R$:
$$ \int_{x\in X} \exp \left( \frac{c_2x^2}2\right) \frac{b_\ell x^\ell}{\ell!} = \begin{cases} 0, & \ell \text{ odd} \\ \text{normalization} \times \frac{b_\ell (c_2/2)^{\ell/2}}{(\ell/2)!}, & \ell \text{ even} \end{cases},$$
(where $\text{normalization}$ involves $\det c_2$, factors of $\sqrt{2\pi}$, and doesn't make really sense in infinite dimensions), and recognizing what all of these exponential generating functions are counting.
When all the dust has settled, what these are counting is (all: disconnected, empty, etc., are allowed) graphs with trivalent-and-higher vertices, mod symmetries. If you throw a $\log$ out in front of the whole integral, then the resulting exponential generating function counts connected graphs, still mod symmetries.
In any case, to actually answer your question for the Dyson series above: the value of each graph is computed by placing a $c_n$ at each vertex of valence $n$ and a $(c_2)^{-1}$ on each edge, and contracting these tensors according to the graph. The symmetry factor is precisely the number of automorphisms of the graph, defined as follows:
Definition: A Feynman graph is a collection $H$ of half edges along with two partitions: one (edges) into sets of size precisely $2$, and the other (vertices) into sets of size at least $3$. (Since each graph will be weighted by the coefficients $c_n$, if you want graphs for a theory with most $c_n = 0$, you can restrict to only graphs with the prescribed valences. If $X$ splits as a direct sum $X = X_1 \oplus X_2$, you can reasonably define graphs with half-edges colored by $X_1,X_2$. If you want to compute an integral with an "observable", you should consider graphs with prescribed "external half-edges".) This is the correct definition, because it gives the correct notion of iso/automorphism. In particular, an isomorphism of Feynman graphs is a bijection of half edges that induces bijections on the partitions. An automorphism is an isomorphism from a Feynman graph to itself. Then the symmetry factor is precisely the number of automorphism in this sense.
(A different notion of "graph" more commonly used in mathematics is that of an "adjacency matrix", which is a $\mathbb Z_{\geq 0}$-valued matrix indexed by the vertices of the graph. The natural notion of isomorphism of such things is a bijection of vertices that induces an equality of adjacency matrices, but this is not the right one for counting symmetries of Feynman diagrams, as, for example, it gives "$1$" for each of the figure-eight and tadpole diagrams.)
For small graphs, the automorphism group splits as a direct product. For example, for the figure eight, the automorphism group is $(\mathbb Z/2) \ltimes (\mathbb Z/2)^2$, where the $(\mathbb Z/2)^2$ acts as flips of the two lobes of the figure eight, and the left-hand $\mathbb Z/2$ switches the two lobes. In a "$\phi^3$" theory, the theta graph has automorphism group $(\mathbb Z/2) \ltimes S_3$, where $S_3$, the symmetric group on three objects, acts to permute the edges, and the $\mathbb Z/2$ switches the two vertices. So the figure eight has $2\times 2^2 = 8$ symmetries, and the theta has $2\times 3! = 12$ symmetries. For small graphs, you can really just read off these semidirect-product decompositions: self-loops contribute factors of $\mathbb Z/2$, each collection of $k$ parallel edges contributes a factor of $S_k$, if you permute vertices, that's another factor, etc. For large graphs, the computations become much harder, and I think I can find a graph whose symmetry group is any finite group you want (but don't quote me on that thought). A slow way to do the computation is to consider all possible permutations of half edges, and see if they induce automorphism. This is slow because there are $(\text{lots})!$ many such permutations. The whole point of working with Dyson series is to avoid this. Fortunately, no one ever computes large graphs, because even just contracting all the tensors is so hard for those.
Finally, I cannot help but mention one more fact you probably know. These Dyson series almost never have positive radius of convergence in the perturbation parameters. In particular, even after dividing by symmetry factors, there are still a lot of graphs, and so the coefficients grow as $n!$. A good exercise is to compute the Dyson series for the conditionally-convergent Riemann integral $\int_{\mathbb R} \exp \frac i \hbar \bigl( \frac{x^2}2 + \lambda \frac{x^3}6 \bigr){\rm d}x$. (The $i$ is there to make the integral conditionally converge; I tend to work with $\hbar$ as my perturbation parameter.) Then the Dyson series grows as $\sum \frac{(6n!)}{(2n)!\,(3n)!} \alpha^n$, where $\alpha$ is linear in $\lambda\hbar$ (something like $\alpha = \lambda\hbar / 6$). The point is that by Stirling's formula $\frac{(6n!)}{(2n)!\,(3n)!} \approx \beta^n n!$ for some positive $\beta$. In any case, this is not at all surprising. Integrals of the form $\int \exp \sum c_nx^n/n!$ are essentially never analytic at $0,\infty$ in the coefficients — for example, Gauss's formula $\int_{\mathbb R}\exp(-a^{-1}x^2/2){\rm d}x = \sqrt{2\pi a}$ has a ramified point at $a = 0,\infty$, and so the Riemann integral extends from its domain of convergence (${\rm Re}(a) > 0$) to a double-valued complex function in $a$.
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.
Best Answer
From the sound of it, you are reading Costello's book.
In point particle QFT the stable graphs you are referring to can be thought of as the paths of particles through spacetime the vertices are where the particles interact.
Now in the Deligne-Mumford case you can think of taking the stable graph and "thickening" it and replacing the particles with little loops of string to obtain a Riemann surface. The interactions then can be thought of as corresponding to the joining and splitting of these little loops of string.
It happens to be the case that in bosonic string theory one is in this second "thickened" case and the symmetries of the theory are such that in doing the path integral one ends up integrating over the Deligne-Mumford spaces of such Riemann surfaces. Thus, the Deligne-Mumford stratification is of use there. Furthermore, you can relate Deligne-Mumford spaces of such Riemann surfaces to point particle QFT by simply taking the limit of infinite string tension.