Let's assume I have Euclidean vector space ${\rm R}^n$ over real numbers ${\rm R}$ with the inner product $\left< \cdot | \cdot \right>$. Then, for a symmetric positive definite real $n\times n$ matrix $A$, we have
$$ \int d^n x~ e^{-\frac{1}{2} \left< x | A | x \right> + \left< J | x \right>} = \left( \frac{(2\pi)^n}{\det A} \right)^{1/2} e^{\frac{1}{2} \left< J | A^{-1} | J \right> }, $$ where $x$ and $J$ are Euclidean vectors.
How to prove from above that $$ \int d^n x~ e^{-\frac{1}{2} \left< x | A | x \right> + i \left< J | x \right>} = \left( \frac{(2\pi)^n}{\det A} \right)^{1/2} e^{-\frac{1}{2} \left< J | A^{-1} | J \right> } ~~ ?$$
I guess one should not be able to use simple transformation $J \to -i J$ since inner product of $\left< -i J | x \right>$ is not well defined?
Well, it amounts to what you said, although you still use the usual real inner product with the $i$ in there, and that gives you a factor of $-1$. To understand what's going on, you have to look at (the beginning of) the derivation of the first formula — you complete the square, substituting $x-A^{-1}J$ for $x$. Now you do the same thing with $iJ$. To finish the argument, you need to use a little bit of complex analysis to argue that the integral $\displaystyle\int_{-\infty}^\infty e^{-x^2}\,dx$ does not change if you push in the imaginary direction and consider $\displaystyle\int_{-\infty+ic}^{\infty+ic} e^{-z^2}\,dz$.