  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(\alpha^2)$ initial state corrections to $e^+e^-$ annihilation into a virtual
boson for large cms energies $s \gg m_e^2$ 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.
