The scale evolution of parton distribution functions is governed by the splitting functions. One of the most efficient methods is to extract splitting functions via the computations of off-shell operator matrix elements, where the physical twist-two operators mix with some unknown gauge-variant operators. We proposed a new method to systematically extract the counterterm Feynman rules resulting from those gauge-variant operators. As a first application of the new method, we rederived the three-loop singlet unpolarized splitting functions. As another application, we obtained a new result for the four-loop pure-singlet splitting functions with two closed fermionic loops. We also presented a new result for the $N_f\,C_F^3$ contribution to the four-loop non-singlet splitting function, where the gauge-variant operators are absent.