We perform a state-of-the-art Bayesian statistical inference to extract the equation of state (EoS) and heavy quark transport coefficient of quark-gluon plasma (QGP) created in relativistic heavy-ion collisions. Based on our quasi-particle linear Boltzmann transport (QLBT) model combined with (3+1)-dimensional viscous hydrodynamic simulation of the QGP and a hybrid fragmentation-coalescence approach for heavy flavor hadronization, the QGP EoS and the heavy quark transport coefficient can be simultaneously constrained from the D meson $R_{AA}$ and $v_2$ data at RHIC and the LHC, both consistent with the lattice QCD results. Our study provides a novel perspective on constraining the QGP EoS from hard probe observables.