We present preliminary results of a calculation of the QCD renormalization constants (RCs) for quark bilinear operators, evaluated non-perturbatively on the lattice in the RI'-MOM scheme. The calculation is performed by using dedicated ensembles with $N_f=4$ degenerate dynamical twisted mass (clover) fermions and the Iwasaki gauge action. A detailed analysis is reported, with emphasis on the control or subtraction of the hadronic contaminations occurring in the lattice estimators of RCs and a check of proper scaling with $a^2$ of the final results. Such a careful study of systematic errors is the counterpart of the high statistical precision reached by current calculations of RCs in the RI'-MOM scheme and is important in order to quote accurate results in phenomenological applications, such as the computation of quark masses.