Understanding the behavior of large atomic nuclei (heavy ions) in high-energy collisions has been the focus of a concerted research effort over the past 10-15 years. Much of the latest progress in the field has centered around transverse momentum-dependent (or "unintegrated") parton distributions: in particular the prediction of the high-energy behavior of these distributions, in the form of the Balitsky-JIMWLK equations, and the development of the hybrid factorization framework, which connects the unintegrated parton distributions to predictions for experimentally measured cross sections. With the advent of high-energy proton-nucleus collisions at RHIC and the LHC, we are able to experimentally test these predictions for the first time. In this dissertation, I show two case studies of these predictions, to illustrate the use of the hybrid factorization at leading and next-to-leading order.First, as a simple example, I analyze the azimuthal angular correlation for a Drell-Yan process, the production of a lepton pair with an associated hadron. The correlation for back-to-back emission turns out to be determined by the low-momentum region of the unintegrated gluon distribution, and the correlation for parallel emission is determined by the high-momentum region. Accordingly, a proper prediction of the correlation at all angles requires a gluon distribution with physically realistic behavior at both high and low momenta. Furthermore, the properties of the central double peak that emerges in Drell-Yan production can provide some insight into the form of the gluon distribution.I'll then describe a numerical calculation of the cross section for inclusive hadron production, which incorporates all corrections up to next-to-leading order in the strong coupling. This calculation illustrates several obstacles presented by subleading terms, including the removal of divergences by renormalizing the integrated and unintegrated parton distributions. The results of the calculation are negative at high transverse momentum, which is surprising but may be mathematically reasonable, since the perturbative approximation to the cross section may break down under those kinematic conditions. However, it may be possible to make meaningful predictions for the nuclear modification ratio R_pA despite the negative cross section.Moving beyond next-to-leading order, it may be possible to cure the negativity of the inclusive hadron cross section by altering the formulas used. I'll show two possible methods of doing so: first, a straightforward resummation of selected higher-order terms corresponding to gluon loop diagrams is able to mitigate the negativity, though it requires some alterations of unclear theoretical origin. A more promising alternative seems to be use of exact kinematic definitions, incorporating terms which disappear in the infinite-energy limit; this constrains the kinematics to eliminate the region of phase space which most strongly contributes to the negativity. In this way, the calculation can be adapted to produce reasonable results at high transverse momentum.