Using a model based on the Color Glass Condensate framework and the dilute-dense factorization, we systematically study the azimuthal angular correlations between a heavy flavor meson and a light reference particle in proton-nucleus collisions. The obtained second harmonic coefficients (also known as the elliptic flows) for $J/psi$ and $D^0$ agree with recent experimental data from the LHC. We also provide predictions for the elliptic flows of $Upsilon$ and $B$ meson, which can be measured in the near future at the LHC. This work can shed light on the physics origin of the collectivity phenomenon in the collisions of small systems.