We perform a non-perturbative study of the scale-dependent renormalisation factors of a complete set of dimension-six four-fermion operators. The renormalisation-group (RG) running is determined in the continuum limit for a specific Schrdinger 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 $mathcal{O}(a)$-improved Wilson fermions. For one of these operators, the computation had been performed in ref. [1]; 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 $Delta 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 ref. [2] 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 (matrix-SSFs), 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 $mathcal{O}(Lambda_{rm QCD})$ to $mathcal{O}(M_W)$.