A new implementation of the general PQCD formalism of Collins, including heavy quark mass effects, is described. Important features that contribute to the accuracy and efficiency of the calculation of both neutral current (NC) and charged current (CC) processess are explicitly discussed. This new implementation is applied to the global analysis of the full HERA I data sets on NC and CC cross sections, with correlated systematic errors, in conjunction with the usual fixed-target and hadron collider data sets. By using a variety of parametrizations to explore the parton parameter space, robust new parton distribution function (PDF) sets (CTEQ6.5) are obtained. The new quark distributions are consistently higher in the region x ~ 10−3 than previous ones, with important implications on hadron collider phenomenology, especially at the LHC. The uncertainties of the parton distributions are reassessed and are compared to the previous ones. A new set of CTEQ6.5 eigenvector PDFs that encapsulates these uncertainties is also presented.