Numerical Stochastic Perturbation Theory (NSPT) enables very high order computations in Lattice Gauge Theories. We report on the determination of the gluon condensate from lattice QCD measurements of the basic plaquette. This is a long standing problem, which was eventually solved a few years ago in pure gauge. In this context NSPT is crucial: it is actually the only tool enabling the subtraction of the power divergent contribution associated to the identity operator in the OPE for the plaquette. This subtraction is actually a delicate issue, since the perturbative expansion of the plaquette is on general ground expected to be an asymptotic one, due to renormalons. This in turn results in ambiguities and the separation of scales in the OPE does not correspond to a separation of perturbative and non-perturbative contributions. All in all, one needs to absorb the ambiguities attached to the perturbative series into the definition of the condensate itself, i.e. one needs a prescription. A possible one amounts to summing the perturbative series up to its minimal term, which means computing up to orders which only NSPT can aim at. Our computation is the first one in QCD, with massless staggered fermions. In order to remove the zero-mode of the gauge field, twisted boundary conditions are adopted for the latter, consistently coupled to fermions in the fundamental representation supplemented with smell degrees of freedom.