Making use of functionals in Martin Siggia Rose formalism
There are two issues at stake here.
First, concerning the question itself :
By definition, $\Gamma$ is the Legendre transform of $W=\ln Z$, $$ \Gamma[\phi]=-W[J]+\phi.J,\\ \frac{\delta \Gamma}{\delta\phi_A}=J_A, $$ where I use condensed notations $J_A$, $\phi_A$, $J.\phi=J_A \phi_A$, etc. with $A$ collecting space-time coordinates, indices, and so on.
We also know that $\frac{\delta W}{\delta J_A}=\phi_A$ and $\frac{\delta^2 W}{\delta J_A\delta J_B}=\frac{\delta \phi_A}{\delta J_B}=G_{AB}$.
Thus, we find that $$ \frac{\delta^2 \Gamma}{\delta\phi_A\delta \phi_B}=\frac{\delta J_A}{\delta\phi_B}=\left(\frac{\delta \phi_B}{\delta J_A}\right)^{-1}=(G^{-1})_{AB}. $$
Up to the Fourier transform, this is the equation given in the OP. In fact, this equation, and the one one can obtain by taking more functional derivatives, just give us the relationship between the connected correlation functions, and the vertex function, nothing more.
This leads us to the second issue : using the formal definition of $\Gamma$ does not help to compute neither $\Gamma^{(2)}$ nor $G$ (unless we know $W$ of course). We thus need another way to compute $\Gamma$, and then use this to compute $G$.
One possibility is to compute $\Gamma$ as a loop expansion (it is better to compute $\Gamma$ than $W$ because there are less diagrams (only the 1PI) and resummations of the self-energies are already performed explicitly).
There are other kind of approximations, using Ansatz of $\Gamma$ and using some RG equation to find the coefficients (in the context of out-of-equilibrium dynamics, see for instance Phys. Rev. Lett. 92, 195703 (2004))