The physics of gluon saturation or color glass condensate (CGC) has been one of the main driving forces for the future Electron Ion Collider. Significant progress has been made in the theory and phenomenology for computing physics observables measured at RHIC and LHC in the past decades. However, the higher-order perturbative calculations in the CGC formalism still remains a bit elusive, especially in comparison with the conventional QCD collinear factorization in the dilute regime. In this talk, using single hadron production in proton-nucleus collisions as an example, I point out some of the interesting difficulties and demonstrate how this can be solved. For example, in the standard collinear factorization, the natural hard scale will automatically arise from a higher-order calculation, while choosing the natural rapidity scale for small-x remains quite tricky even with explicit higher-order calculation. I further show how to perform threshold resummation and demonstrate how this would solve the well-known negative cross section problem for this process.