We perform a non-perturbative study of the scaledependent renormalisation factors of a complete set of dimension-six four-fermion operators without power subtractions. The renormalisation-group (RG) running is determined in the continuum limit for a specific Schrödinger Functional (SF) renormalisation scheme in the framework of lattice QCD with two dynamical flavours (N f = 2). The theory is regularised on a lattice with a plaquette Wilson action and O(a)-improved Wilson fermions. For one of these operators, the computation had been performed in Dimopoulos et al. (JHEP 0805, 065 (2008). arXiv:0712.2429); the present work completes the study for the rest of the operator basis, on the same simulations (configuration ensembles). The related weak matrix elements arise in several operator product expansions; in F = 2 transitions they contain the QCD long-distance effects, including contributions from beyond-Standard Model (BSM) processes. Some of these operators mix under renormalisation and their RG-running is governed by anomalous dimension matrices. In Papinutto et al. (Eur Phys J C 77(6), 376 (2017). arXiv:1612.06461) the RG formalism for the operator basis has been worked out in full generality and the anomalous dimension matrix has been calculated in NLO perturbation theory. Here the discussion is extended to the matrix step-scaling functions, which are used in finite-size recursive techniques. We rely on these matrix-SSFs to obtain non-perturbative estimates of the operator anomalous dimensions for scales ranging from O(QCD) to O(M W).