In QCD the anomalous dimensions of gauge invariant operators of twist 2 play a key role, because they control the scale dependence of the parton distribution functions. Notably, the flavour singlet operators, such as those associated to the gluon distribution, mix under renormalisation with a set of unphysical operators, also known as aliens. Missing this effect leads to wrong results already for the two-loop anomalous dimensions. The correct renormalisation of gluonic operators is an important step towards the computation of the scale evolution of flavour singlet parton distributions, which is now required to 4 loops. Leveraging both the background field method and an enhanced BRST symmetry, we construct the required ghost and alien operator basis up to 4 loops for arbitrary mass dimensions. Furthermore, we extract anomalous dimensions at 4 loops, for the physical operators of mass-dimension 4 and 6, and at 3 loops for mass-dimension 8.