*Not a complete answer but too long for a comment:*

Let $I=\left]0,1\right[$ be the unit interval. Let $$F:I^2\to M, (x,y)\mapsto \left(x,\frac{x}{x+y}\right),$$

where $$M=\{(a,b)\in I\times I\mid \frac{a}{a+1}\le b\}=\{(a,b)\in I\times I\mid a\le\frac{b}{1-b}\}.$$

(Here I have the convention that $\frac{b}{b-1}=\infty$ for $b=1$.)

Note that $F$ is a diffeomorphism with inverse

$$F^{-1}(a,b)=\left(a, \frac{a\cdot(1-b)}b\right).$$

So, where $\operatorname{Jac}$ denotes the Jacobian determinant,

$$|(\operatorname{Jac}F^{-1})(a,b)|=\frac{a}{b^2}.$$

Use the following result from probability Theory (here, $f_Z$ denotes the density of the random vector $Z$):

$$\bbox[15px,border:1px groove navy]{f_{F\circ(X,Y)}(a,b)=f_{(X,Y)}(F^{-1}(a,b))\cdot|(\operatorname{Jac}F^{-1})(a,b)|\cdot1_{\operatorname{Im}(F)}(a,b).}$$

In our case, $$f_{(X,\frac{X}{X+Y})}(a,b)=f_{(X,Y)}\left(a,\frac{a\cdot(1-b)}b\right)\cdot\frac{a}{b^2}\cdot 1_{M}(a,b).$$

Note that, since $X$ and $Y$ are independent, the joint density is obtained simply as a product of the individual densities of a Beta distribution. Say $X\sim \operatorname{Beta}(\alpha,\beta)$ and $Y\sim\operatorname{Beta}(\alpha,\beta)$. Let $$c=\operatorname{B}(\alpha,\beta)\cdot \operatorname{B}(\gamma, \delta),$$ where $\operatorname{B}$ is the Beta function.

Then

\begin{split}
f_{\frac{X}{X+Y}}(b)&=\int_\mathbb R f_{(X,\frac{X}{X+Y})}(a,b)\,\mathrm da\\
&=1_{I}(b)\cdot c^{-1} b^{\gamma-3} (1-b)^{\delta-1}\cdot\int_0^{\min(1,\frac{b}{1-b})} a^\alpha(1-a)^{\beta-1}\cdot\left(1+\frac{a\cdot(b-1)}b\right)^{\delta-1}\,\mathrm da.
\end{split}