We present a novel strategy to renormalize lattice operators in QCD+QED, including first order QED corrections to the non-perturbative evaluation of QCD renormalization constants. Our procedure takes systematically into account the mixed non-factorizable QCD+QED effects which were neglected in previous calculations, thus significantly reducing the systematic uncertainty on renormalization corrections. The procedure is presented here in the RI'-MOM scheme, but it can be applied to other schemes (e.g. RI-SMOM) with appropriate changes. We discuss the application of this strategy to the calculation of the leading isospin breaking corrections to the leptonic decay rates $\Gamma(\pi_{\mu 2})$ and $\Gamma(K_{\mu 2})$, evaluated for the first time on the lattice. The precision in the matching to the $W$-regularization scheme is improved to $\mathcal{O}(\alpha_{em}\alpha_s(M_W))$ with respect to previous calculations. Finally, we show the updated precise result obtained for the Cabibbo-Kobayashi-Maskawa matrix element $|V_{us}|$.