We present an algorithm for constructing the Wilson operator product expansion (OPE) for perturbative interacting quantum field theory in general Lorentzian curved spacetimes, to arbitrary orders in perturbation theory. The remainder in this expansion is shown to go to zero at short distances in the sense of expectation values in arbitrary Hadamard states. We also establish a number of general properties of the OPE coefficients: (a) they only depend (locally and covariantly) upon the spacetime metric and coupling constants, (b) they satisfy an associativity property, (c) they satisfy a renormalization group equation, (d) they satisfy a certain microlocal wave front set condition, (e) they possess a "scaling expansion". The latter means that each OPE coefficient can be written as a sum of terms, each of which is the product of a curvature polynomial at a spacetime point, times a Lorentz invariant Minkowski distribution in the tangent space of that point. The algorithm is illustrated in an example.