Splitting functions are the non-local anomalous dimensions of parton distribution functions, and as such important ingredients for as good as all LHC calculations. It has been known for a long time that the Mellin moments of these splitting functions correspond to anomalous dimensions of local leading-twist operators which appear in the lightcone operator product expansion. A long-standing problem which prevented the use of this framework beyond two loops was due to mixing with gauge-variant operators, which were difficult to classify. I will discuss a method based on generalised BRST and anti-BRST symmetry to derive the general form of this basis. I will then present a number 4-loop results obtained with this method.