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}(λ_k - λ_j) (λ_k^θ- λ_j^θ)$. 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 $θ\in \mathbb Z^+$ they can be realised as the eigenvalue PDF of certain random matrices with Gaussian entries. For general $θ> 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 $θ\in \mathbb Z^+$ the Laguerre and Jacobi cases are closely related to the squared singular values for products of $θ$ 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.