Integral involving incomplete beta function I have the following integral,
$$\int_{0}^1x^{a-1}(1-x)^{b-1}B_x(c,d)dx$$
where $B_x(c,d) = \int_{0}^xt^{c-1}(1-t)^{d-1}dt$ is the incomplete beta function, and $a,b,c,d>0$. 
Question: Does this have a closed form?

My attempt:


*

*First, playing around in Wolfram Alpha makes me think that there may be a (simple?) closed form: example 1 and example 2. The second example can also be written as $\int x^2(1-x)B_x(12,2)dx = \frac{1}{4}(B_x(16,2) - x^4B_x(12,2))+\frac{1}{3}(x^3B_x(12,2)-B_x(15,2))$.

*It seems there is a reduction when $a=c$, $b=d$ as
$$\int_0^1 x^{a-1}(1-x)^{b-1}B_x(a,b)dx=\frac{1}{2}\left(\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\right)^2$$
where $\Gamma(a)$ is the gamma function. However, for the case when $a\neq c$, $b\neq d$, things are not quite as clear for me.

*I found a similar question here but it has not been answered.

*I also found this but I'm not sure if it is useful.

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$.
