We describe the calculation of the two-loop massive operator matrix elements with massive external fermions in QED. We investigate the factorization of the O(α 2 ) initial state corrections to e + e − annihilation into a virtual boson for large cms energies s ≫ m 2 e into massive operator matrix elements and the massless Wilson coefficients of the Drell-Yan process adapting the color coefficients to the case of QED, as proposed by Berends et. al. in Ref. [1]. Our calculations show explicitly that the representation proposed in Ref.[1] works at one-loop order and up to terms linear in ln(s/m 2 e ) at two-loop order. However, the two-loop constant part contains a few structural terms, which have not been obtained in previous direct calculations.