This paper presents a next-to-leading order (NLO) global DGLAP analysis of nuclear parton distribution functions (nPDFs) and their uncertainties. The analysis includes three types of experimental data: deep inelastic lepton-nucleus scattering, Drell-Yan dilepton production in p+A collisions, and inclusive pion production in d+Au and p+p collisions at RHIC. The results show that these data can be well described in a conventional collinear factorization framework. The analysis also finds that the nuclear modification factor R_{dAu} of the pion p_T-spectrum at midrapidity retains sensitivity to the gluon distributions, providing evidence for shadowing and EMC-effect in the nuclear gluons. The Hessian method is used to quantify the nPDF uncertainties, which originate from the uncertainties in the data. The obtained NLO and LO nPDFs and the corresponding error sets are collected in our new release called EPS09. These results should find applications in precision analyses of the signatures and properties of QCD matter at the LHC and RHIC. The paper also discusses the comparison of the results with the data and with earlier global analyses, and presents a concrete application of the obtained nPDFs and shows how to best estimate the PDF-related uncertainties for nuclear cross-sections. The paper concludes with a summary and a brief outlook.This paper presents a next-to-leading order (NLO) global DGLAP analysis of nuclear parton distribution functions (nPDFs) and their uncertainties. The analysis includes three types of experimental data: deep inelastic lepton-nucleus scattering, Drell-Yan dilepton production in p+A collisions, and inclusive pion production in d+Au and p+p collisions at RHIC. The results show that these data can be well described in a conventional collinear factorization framework. The analysis also finds that the nuclear modification factor R_{dAu} of the pion p_T-spectrum at midrapidity retains sensitivity to the gluon distributions, providing evidence for shadowing and EMC-effect in the nuclear gluons. The Hessian method is used to quantify the nPDF uncertainties, which originate from the uncertainties in the data. The obtained NLO and LO nPDFs and the corresponding error sets are collected in our new release called EPS09. These results should find applications in precision analyses of the signatures and properties of QCD matter at the LHC and RHIC. The paper also discusses the comparison of the results with the data and with earlier global analyses, and presents a concrete application of the obtained nPDFs and shows how to best estimate the PDF-related uncertainties for nuclear cross-sections. The paper concludes with a summary and a brief outlook.