We study open heavy flavor meson production in proton–nucleus (pA) collisions at RHIC and LHC energies within the Color Glass Condensate framework. We use the unintegrated gluon distribution at small Bjorkenʼs x in the proton obtained by solving the Balitsky–Kovchegov equation with running coupling correction and constrained by global fitting of HERA data. We change the initial saturation scale of the gluon distribution for the heavy nucleus. The gluon distribution with McLerran–Venugopalan model initial condition is also used for comparison. We present transverse momentum spectra of single D and B productions in pA collisions, and the so-called nuclear modification factor. The azimuthal angle correlation of open heavy flavor meson pair is also computed to study the modification due to the gluon saturation in the heavy nucleus at the LHC.