We propose and analyze a time-stepping discontinuous Petrov- Galerkin method combined with the continuous conforming finite element method in space for the numerical solution of time- fractional subdiffusion problems. We prove the existence, uniqueness, and stability of approximate solutions and derive error estimates. To achieve high order convergence rates from the time discretizations, the time mesh is graded appropriately near t = 0 to compensate for the singular ( temporal) behavior of the exact solution near t = 0 caused by the weakly singular kernel, but the spatial mesh is quasi uniform. In the L infinity((0, T); L-2( Omega))- norm, ((0, T) is the time domain and Omega is the spatial domain); for sufficiently graded time meshes, a global convergence of order k(m)+(alpha/2)+h(r+1) is shown, where 0 < alpha < 1 is the fractional exponent, k is the maximum time step, h is the maximum diameter of the elements of the spatial mesh, and m and r are the degrees of approximate solutions in time and spatial variables, respectively. Numerical experiments indicate that our theoretical error bound is pessimistic. We observe that the error is of order k(m)+(alpha/2)+h(r+1), that is, optimal in both variables.