  Muttalib--Borodin ensembles are characterised by the pair interaction term in
the eigenvalue probability density function being of the form $\prod_{1 \le j <
k \le N}(\lambda_k - \lambda_j) (\lambda_k^\theta - \lambda_j^\theta)$. We
study the Laguerre and Jacobi versions of this model --- so named by the form
of the one-body interaction terms --- and show that for $\theta \in \mathbb
Z^+$ they can be realised as the eigenvalue PDF of certain random matrices with
Gaussian entries. For general $\theta > 0$, realisations in terms of the
eigenvalue PDF of ensembles involving triangular matrices are given. In the
Laguerre case this is a recent result due to Cheliotis, although our derivation
is different. We make use of a generalisation of a double contour integral
formula for the correlation functions contained in a paper by Adler, van
Moerbeke and Wang to analyse the global density (which we also analyse by
studying characteristic polynomials), and the hard edge scaled correlation
functions. For the global density functional equations for the corresponding
resolvents are obtained; solving this gives the moments in terms of
Fuss--Catalan numbers (Laguerre case --- a known result) and particular
binomial coefficients (Jacobi case). For $\theta \in \mathbb Z^+$ the Laguerre
and Jacobi cases are closely related to the squared singular values for
products of $\theta$ standard Gaussian random matrices, and truncations of
unitary matrices, respectively. At the hard edge the double contour integral
formulas provide a double contour integral form of the scaled correlation
kernel obtained by Borodin in terms of Wright's Bessel function.
