  Within the framework of generalized collinear factorization in perturbative
QCD (pQCD), we study the effect of initial multiple parton scattering and
induced parton energy loss in Drell-Yan (DY) process in proton-nucleus
collisions. We express the contribution from multiple parton scattering and
induced gluon radiation to the DY dilepton spectra in terms of nuclear modified
effective beam quark distribution functions. The modification depends on the
quark transport parameter in nuclear medium. This is similar to the final-state
multiple parton scattering in deeply inelastic scattering (DIS) off large
nuclei and leads to the suppression of the Drell-Yan cross section in $p+A$
relative to $p+p$ collisions. With the value of quark transport parameter
determined from the nuclear modification of single-inclusive DIS hadron spectra
as measured by the HERMES experiment, we calculate DY spectra in $p+A$
collisions and find the nuclear suppression due to beam parton energy loss
negligible at the Fermilab energy $E_{\rm lab}$=800 GeV in the kinematic region
as covered by the E866 experiment. Most of the observed nuclear suppression of
DY spectra in E866 experiment can be described well by parton shadowing in
target nuclei as given by the EPS08 parameterization. The effect of beam parton
energy loss, however, becomes significant for DY lepton pairs with large beam
parton momentum fraction $x^{\prime}$ or small target parton momentum fraction
$x$. We also predict the DY cross section in $p+A$ collisions at lower beam
proton energy $E_{\rm lab}$=120 GeV and show significant suppression due to
initial state parton energy loss at moderately large $x^{\prime}$ where the
effect of parton shadowing is very small.
