  We present a comprehensive study of the quark sector of $2+1$ flavour QCD,
based on a self-consistent treatment of the coupled system of Schwinger-Dyson
equations for the quark propagator and the full quark-gluon vertex. The
individual form factors of the quark-gluon vertex are expressed in a special
tensor basis obtained from a set of gauge-invariant operators. The sole
external ingredient used as input to our equations is the Landau gauge gluon
propagator with $2+1$ dynamical quark flavours, obtained from studies with
Schwinger-Dyson equations, the functional renormalisation group approach, and
large volume lattice simulations. The appropriate renormalisation procedure
required in order to self-consistently accommodate external inputs stemming
from other functional approaches or the lattice is discussed in detail, and the
value of the gauge coupling is accurately determined at two vastly separated
renormalisation group scales.
  Our analysis establishes a clear hierarchy among the vertex form factors. We
identify only three dominant ones, in agreement with previous results. The
components of the quark propagator obtained from our approach are in excellent
agreement with the results from Schwinger-Dyson equations, the functional
renormalisation group, and lattice QCD simulation, a simple benchmark
observable being the chiral condensate in the chiral limit, which is computed
as $(245\,\textrm{MeV})^3$. The present approach has a wide range of
applications, including the self-consistent computation of bound-state
properties and finite temperature and density physics, which are briefly
discussed.
