The forthcoming p+Pb run at the LHC will provide crucial in formation on the initial state effects of heavy ion collisions and on the gluon saturation phenomena. In turn, most of the saturation inspired phenomenology in heavy ion collisions borrows substantial empiric information from the analysis of e+p data, where abundant high quality data on the small-x kinematic region is available. Indeed, the very precise combined HERA data provides a testing ground in which the relevance of novel QCD regimes, other than the successful linear DGLAP evolution, in small-x inclusive DIS data can be ascertained. We present a study of the dependence of the AAMQS fits, based on the running coupling BK non-linear evolution equations (rcBK), on the fitted dataset. This allows for the identification of the kinematical region where rcBK accurately describes the data, and thus for the determination of its applicability boundary. It also set important constraints to the saturation models used to model the early stages of heavy ion collisions. Finally we compare the rcBK results with NNLO DGLAP fits, obtained with the NNPDF methodology with analogous kinematical cuts. Further, we explore the impact on LHC phenomenology of applying stringent kinematical cuts to the low-x HERA data in a DGLAP fit.