Lattice simulations show that the running coupling constants obtained through QCD vertices differ in the infrared even though their bare values are the same. For instance, at low momenta the strength of the quark-gluon coupling is about twice the ghost-gluon coupling. In particular, none of them diverge at a finite non-zero infrared momentum scale contrarily to what is predicted by standard perturbation theory. Moreover, the ghost-gluon coupling remains moderate even in the infrared. This observation motivates the use of perturbation theory in the ghost-gluon sector in the frame of a massive deformation of QCD Lagrangian in Landau gauge. However, perturbation theory in the quark sector within this massive Lagrangian does not bring as good results as in the pure Yang-Mills case, which is consistent with the larger value of the quark-gluon coupling constant. We propose, then, a controlled systematic expansion in full QCD based in two small parameters: first the Yang-Mills sector couplings and second the inverse of the number of colors (large-Nc limit). This systematic expansion allows us to properly introduced the use of the renormalization group for the rainbow resummation. At leading order, this double expansion leads to the well-known rainbow approximation for the quark propagator whose solution shows spontaneous chiral symmetry breaking for sufficiently large quark-gluon coupling constant. We present here some of the results discussed in \cite{Pelaez:2017bhh,Pelaez:2020ups} and work in progress.