Using the hypergeometric representation of the incomplete Beta function
\begin{equation}
B_x\left( c,d \right)=\frac{x^c}{c}{}_2F_1\left( c,1-d;1+c;x \right)
\end{equation}
the integral can be written as
\begin{align}
I\left( a,b,c,d \right)&=\int_{0}^1x^{a-1}(1-x)^{b-1}B_x(c,d)\,dx\\
&=\frac{1}{c}\int_{0}^1x^{a+c-1}(1-x)^{b-1}{}_2F_1\left( c,1-d;1+c;x \right)\,dx\\
&=\frac{1}{c}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \int_0^1x^{a+c+k-1}(1-x)^{b-1}\,dx\\
&=\frac{1}{c}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \frac{\Gamma(b)\Gamma(a+c+k)}{\Gamma(a+b+c+k)}\\
&=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\sum_{k=0}^\infty \frac{(c)_k(1-d)_k}{(1+c)_kk!} \frac{(a+c)_k}{(a+b+c)_k}\\
&=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\,{}_3F_2\left( 1-d,a+c,c;1+c,a+b+c ;1\right)
\end{align}
Using this identity for the generalized hypergeometric function:
\begin{equation}
{}_3F_2\left( a_1,a_2,a_3;b_1,b_2;1 \right)=\frac{\Gamma(b_1)\Gamma(b_1+b_2-a_1-a_2-a_3)}{\Gamma(b_1-a_1)\Gamma(b_1+b_2-a_2-a_3)}{}_3F_2\left( a_1,b_2-a_2,b_2-a_3;b_2,b_1+b_2-a_2-a_3;1 \right)
\end{equation}
Here we choose $a_1=1-d,a_2=a+c,a_3=c,b_1=a+b+c,b_2=1+c$ to obtain
\begin{align}
I\left( a,b,c,d \right)&=\frac{1}{c}\frac{\Gamma(b)\Gamma(a+c)}{\Gamma(a+b+c)}\frac{\Gamma(a+b+c)\Gamma(b+d)}{\Gamma(a+b+c+d-1)\Gamma(b+1)}\,{}_3F_2\left( 1-d,1-a,1;1+c,1+b ;1\right)\\
&=\frac{1}{bc}\frac{\Gamma(a+c)\Gamma(b+d)}{\Gamma(a+b+c+d-1)}\,{}_3F_2\left( 1,1-a,1-d;1+b,1+c ;1\right)\\
&=\frac{a+b+c+d-1}{bc}B(a+c,b+d)\,{}_3F_2\left( 1,1-a,1-d;1+b,1+c ;1\right)
\end{align}
which gives simple results if $a$ or $d$ are positive integer, as the hypergeometric series is finite: as $1-a\le0$ or $1-b\le0$, in the series definition of the hypergeometric function, there are $\operatorname{min}(a,b)$ terms (numerator of the coefficients cancel after that). One can also check the given result when $a=c$ and $b=d$.
Let $\lambda> 0$ and $\beta = \lambda (\alpha - 1)$. Then
$$
\Gamma (1 - \alpha ,\beta ) \sim \frac{{\beta ^{1 - \alpha } \mathrm{e}^{ - \beta } }}{{\alpha - 1}}\sum\limits_{n = 0}^\infty {\frac{{p_n (\lambda )}}{{(\lambda + 1)^{2n + 1} }}\frac{1}{{(\alpha - 1)^n }}}
$$
as $\alpha\to+\infty$ uniformly in $\lambda > 0$. Here $p_0(\lambda)=1$ and
$$
p_n (\lambda ) = \lambda (1 + \lambda )p'_{n - 1} (\lambda ) - (2n - 1)\lambda p_{n - 1} (\lambda )
$$
for $n\geq 1$. These are polynomials in $\lambda$ of degree $n$.
To achive the best numerical accuracy, stop the series after about $\left\lfloor {\alpha \sqrt {(\lambda + \log \lambda + 1)^2 + \pi ^2 } } \right\rfloor$ terms. See this paper for more details.
Best Answer
We may use the following integral representation for the incomplete gamma function (DLMF ref. see identity 8.6.4):
Given $n\in\mathbb{N}^{+}\land a,b,\epsilon\in\mathbb{R}\land1<a\land1\le b\le n\land0<\epsilon$, define $I_{n}{\left(a,b;\epsilon\right)}$ via the integral
$$I_{n}{\left(a,b;\epsilon\right)}=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)\left[a\left(-ix\right)^{a}\,\Gamma{\left(-a,-ix\right)}\right]^{n}}{x^{2}+\epsilon^{2}}.$$
Then,
$$\begin{align} I_{n}{\left(a,b;\epsilon\right)} &=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)\left[a\left(-ix\right)^{a}\,\Gamma{\left(-a,-ix\right)}\right]^{n}}{x^{2}+\epsilon^{2}}\\ &=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)}{x^{2}+\epsilon^{2}}\left[\frac{e^{ix}}{\Gamma{\left(a\right)}}\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\ &=\frac{1}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{2ie^{-ibx}\left(bx-i\right)}{x^{2}+\epsilon^{2}}e^{inx}\left[\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\ &=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\left[\int_{0}^{\infty}\frac{t^{a}e^{-t}}{t-ix}\,\mathrm{d}t\right]^{n}\\ &=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\prod_{k=1}^{n}\int_{0}^{\infty}\mathrm{d}t_{k}\,\frac{t_{k}^{a}e^{-t_{k}}}{t_{k}-ix}\\ &=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{x^{2}+\epsilon^{2}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,e^{-\sum_{k=1}^{n}t_{k}}\prod_{k=1}^{n}\left(\frac{t_{k}^{a}}{t_{k}-ix}\right)}\\ &=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\left(\prod_{k=1}^{n}t_{k}^{a}\right)e^{-\sum_{k=1}^{n}t_{k}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(bx-i\right)e^{i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}-ix\right)}}\\ &=\small{\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}}\\ &=:\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\,f_{n}{\left(b;\epsilon\right)},\\ \end{align}$$
where for $n\in\mathbb{N}^{+}\land b,\epsilon\in\mathbb{R}\land1\le b\le n\land0<\epsilon$ we've defined the auxiliary function denoting the innermost integration,
$$f_{n}{\left(b;\epsilon\right)}:=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}.$$
For the integration over $x$, I will appeal to a well-suited proposition from Gradshteyn's Tables. Now, I usually try to make my posts as self-contained as possible by proving any non-trivial lemmas I plan on using, rather than taking the lazy way out and simply citing the result from an outside source. But I can't resist doing so, partly out of laziness, but mostly just to show off how Gradshteyn can on occasion be eerily clairvoyant.
The following partial fraction decomposition is easily verified:
$$\frac{-bx-i}{x^{2}+\epsilon^{2}}=\frac{1}{2i\epsilon}\left[\frac{\left(1-b\epsilon\right)}{\epsilon-ix}+\frac{\left(1+b\epsilon\right)}{\epsilon+ix}\right].$$
Then, assuming $0<\Re{\left(t_{k}\right)}\land b<n$,
$$\begin{align} f_{n}{\left(b;\epsilon\right)} &=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{\left(-bx-i\right)e^{-i\left(n-b\right)x}}{\left(x^{2}+\epsilon^{2}\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\ &=\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{1}{2i\epsilon}\left[\frac{\left(1-b\epsilon\right)}{\epsilon-ix}+\frac{\left(1+b\epsilon\right)}{\epsilon+ix}\right]\frac{e^{-i\left(n-b\right)x}}{\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\ &=\frac{1-b\epsilon}{2i\epsilon}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{e^{-i\left(n-b\right)x}}{\left(\epsilon-ix\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\ &~~~~~+\frac{1+b\epsilon}{2i\epsilon}\int_{-\infty}^{\infty}\mathrm{d}x\,\frac{e^{-i\left(n-b\right)x}}{\left(\epsilon+ix\right)\prod_{k=1}^{n}\left(t_{k}+ix\right)}\\ &=\frac{1-b\epsilon}{2i\epsilon}\cdot\frac{2\pi e^{-\left(n-b\right)\epsilon}}{\prod_{k=1}^{n}\left(\epsilon+t_{k}\right)}.\\ \end{align}$$
Thus,
$$\begin{align} I_{n}{\left(a,b;\epsilon\right)} &=\frac{2i}{\left[\Gamma{\left(a\right)}\right]^{n}}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\,f_{n}{\left(b;\epsilon\right)}\\ &=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\frac{\left(\prod_{k=1}^{n}t_{k}^{a}\right)}{\exp{\left(\sum_{k=1}^{n}t_{k}\right)}}\cdot\frac{1}{\prod_{k=1}^{n}\left(\epsilon+t_{k}\right)}\\ &=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\int_{[0,\infty)^{n}}\mathrm{d}^{n}\mathbf{t}\,\prod_{k=1}^{n}\left(\frac{t_{k}^{a}e^{-t_{k}}}{\epsilon+t_{k}}\right)\\ &=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\left[\int_{0}^{\infty}\mathrm{d}t\,\left(\frac{t^{a}e^{-t}}{\epsilon+t}\right)\right]^{n}\\ &=\frac{2\pi\left(1-b\epsilon\right)e^{-\left(n-b\right)\epsilon}}{\left[\Gamma{\left(a\right)}\right]^{n}\epsilon}\left[\epsilon^{a}e^{\epsilon}\,\Gamma{\left(1+a\right)}\,\Gamma{\left(-a,\epsilon\right)}\right]^{n}\\ &=2\pi a^{n}\left(1-b\epsilon\right)\epsilon^{na-1}e^{b\epsilon}\left[\Gamma{\left(-a,\epsilon\right)}\right]^{n}.\\ \end{align}$$