© 2022 authors. Published by the American Physical Society. Published by the American Physical Society under the terms of the ''https://creativecommons.org/licenses/by/4.0/''Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP3.We reconstruct ghost and gluon spectral functions in 2+1 flavor QCD with Gaussian process regression. This framework allows us to largely suppress spurious oscillations and other common reconstruction artifacts by specifying generic magnitude and length scale parameters in the kernel function. The Euclidean propagator data are taken from lattice simulations with domain wall fermions at the physical point. For the infrared and ultraviolet extensions of the lattice propagators as well as the low-frequency asymptotics of the ghost spectral function, we utilize results from functional computations in Yang-Mills theory and QCD. This further reduces the systematic error significantly. Our numerical results are compared against a direct real-time functional computation of the ghost and an earlier reconstruction of the gluon in Yang-Mills theory. The systematic approach presented in this work offers a promising route toward unveiling real-time properties of QCD.