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}|$.