  In this paper, we consider Muttalib-Borodin ensemble of Laguerre type, a
determinantal point process over $[0,\infty)$ which depends on the varying
weights $x^{\alpha}e^{-nV(x)}$, $\alpha>-1$, and a parameter $\theta$. For
$\theta$ being a positive integer, we derive asymptotics of the associated
biorthogonal polynomials near the origin for a large class of potential
functions $V$ as $n\to \infty$. This further allows us to establish the hard
edge scaling limit of the correlation kernel, which is previously only known in
the special cases and conjectured to be universal. Our proof is based on the
Deift/Zhou nonlinear steepest descent analysis of two $1 \times 2$
vector-valued Riemann-Hilbert problems that characterize the biorthogonal
polynomials and the explicit constructions of $(\theta+1)\times(\theta+1)$
local parametrices near the origin in terms of the Meijer G-functions.
