In this paper we construct the Wilson short distance operator product expansion for the gluon, quark and ghost propagators in QCD, including operators of dimension two and three, namely, A^2, m^2, m A^2, \ovl{\psi} \psi and m^3. We compute analytically the coefficient functions of these operators at three loops for all three propagators in the general covariant gauge. Our results, taken in the Landau gauge, should help to improve the accuracy of extracting the vacuum expectation values of these operators from lattice simulation of the QCD propagators.