  We demonstrate the QCD factorization for inclusive hadron production in $pA$
collisions in the saturation formalism at one-loop order, with explicit
calculation of both real and virtual gluon radiation diagrams. The collinear
divergences associated with the incoming parton distribution of the nucleon and
the outgoing fragmentation function of the final state hadron, as well as the
rapidity divergence with small-$x$ dipole gluon distribution of the nucleus are
factorized into the splittings of the associated parton distribution and
fragmentation functions and the energy evolution of the dipole gluon
distribution function. The hard coefficient function is evaluated at one-loop
order, and contains no divergence.
