OLDER EDIT. (Elementary derivation) I realized that my original answer was actually overkill for the question. The said integral in question follows from the definition of the multivariate Gamma function
\begin{equation*}
\Gamma_p(a) := \int_{A > 0} \exp(-\mbox{tr}(A))\det(A)^{a-(p+1)/2}(dA),
\end{equation*}
where $\Re(a)>(p-1)/2$.
From this it follows (by a change of variables) that for a positive definite matrix $S$,
\begin{equation*}
\int_{A > 0} \exp(-\mbox{tr}(S^{-1}A))\det(A)^{a-(p+1)/2}(dA) = \Gamma_p(a)\det(S)^a,
\end{equation*}
so that with $S=(I-U)^{-1}$ we obtain the integral in question.
Of course, to complete the picture it may be helpful to express $\Gamma_p(a)$ in more elementary terms. Chapter 2 of Muirhead's book provides these details. I cite the result that provides this expression.
Theorem (Muirhead (1982), Thm 2.1.2) Let $\Re(a) > (p-1)/2$. Then,
$$
\Gamma_p(a) = \pi^{p(p-1)/4}\prod_{j=1}^p \Gamma(a - (j-1)/2)
$$
(Hint: To prove the above, write the Cholesky decomposition $A=T'T$, with these change of variables, the original Gamma integral factorizes into a product of Gaussian and Gamma integrals.)
The part that I recalled below provides yet another representation that expresses the multiplicative determinantal lhs in terms of an infinite sum.
OLDER STUFF
This is actually somewhat classical knowledge. Here are two related pointers.
A partition $\tau=(t_1,\ldots,t_m)$ is a vector of nonnegative integers listed in increasing order, and $|\tau|$ denotes $t_1+\cdots+t_m$. The generalized Pochhammer symbol $(a)_\tau$ is defined as
\begin{equation*}
\newcommand{\risingf}[2]{{{#1}}^{\overline{{#2}}}}
(a)_\tau := \frac{\Gamma_d(a+\tau)}{\Gamma_d(a)} = \prod_{l=1}^m \risingf{\bigl(a - \tfrac{1}{2}(l-1)\bigr)}{t_l}
\end{equation*}
Let $C_\tau(X)$ be the Zonal Polynomial with signature partition $\tau$. Then, the following representation exists
For a matrix $U$ satisfying $\|U\| < 1$, we have the following "binomial-theorem"
\begin{equation}
\frac{1}{|I-U|^a} = \sum_{k\ge 0}\sum_{|\tau| = k} \frac{(a)_\tau C_\tau(U)}{k!}.
\end{equation}
Using representations for these Zonal polynomials, one can obtain the integral representation mentioned in the original post.
More directly, you can look at Chapter 7 of R. Muirhead, "Aspects of Multivariate Statistical Theory", where you'll see that actually, $|I-U|^{-a}={}_1F_0(a;U)$, a matrix argument hypergeometric function. I've to run now, if I get a chance I'll clean up my answer and fill in the details.