You can follow large deviation-type estimates from random matrix theory, starting with Ben Arous & Guionnet's paper: https://link.springer.com/article/10.1007/s004400050119

You will eventually obtain that your Selberg-type integral behaves like
$$
\frac{1}{(n!)^2}e^{-n^2 E_*}
$$
where $E_*$ is the minimum of the functional
\begin{multline}
E(\mu,\nu)=\iint \log\frac1{|x-y|}d\mu(x)d\mu(y)\\+\iint \log\frac1{|x-y|}d\nu(x)d\nu(y)\\-2\iint\log\frac1{|1-xy|}d\mu(x)d\nu(y)
\end{multline}
where $(\mu,\nu)$ ranges over pairs of Borel probability measures supported in $[0,1]$. I guess you can work out the existence and unicity of a minimizer for $E$ and obtain Euler-Lagrange equations to characterize this minimiser and, hopefully, compute explicitly the minimal value $E_*$.

There is quite some technical work to do to fill the gaps, and sorry for the self-advertisement, but if you follow this approach let me tell you that I had to struggle with a similar two type particle problem in the paper http://www.worldscientific.com/doi/abs/10.1142/S2010326312500165 written with Arno Kuijlaars. I hope this can help.

By the way, minimizing $E$ is referred as to a vector equilibrium problem in potential theory.

ADDENDUM: After thinking further, it seems more natural to make the changes of variables $s_i=1/x_i$, which leads to the same asymptotics but with $E_*$ the minimum of
\begin{multline}
\tilde E(\mu,\nu)=\iint \log\frac1{|x-y|}d(\mu-\nu)(x)d(\mu-\nu)(y)+\int_1^\infty \log(x)d\mu(x)
\end{multline}
where $\mu,\nu$ are Borel probability measures living on $[1,\infty)$ and $[0,1]$ respectively.