Measurements of open charm and beauty production cross sections in deep inelastic $ep$ scattering at HERA from the H1 and ZEUS Collaborations are combined. Reduced cross sections are obtained in the kinematic range of negative four-momentum transfer squared of the photon $2.5$ GeV$^2<Q^2<2000$ GeV$^2$ and Bjorken scaling variable $3cdot10^{-5}<x_{text{Bj}}<5cdot10^{-2}$. The combination method accounts for the correlations of the statistical and systematic uncertainties among the different datasets. Perturbative QCD calculations are compared to the combined data. A next-to-leading order QCD analysis is performed using these data together with the combined inclusive deep inelastic scattering cross sections from HERA. The running charm- and beauty-quark masses are determined as $m_c(m_c) = 1.290^{+0.046}_{-0.041}text{(exp/fit)}^{+0.062}_{-0.014}text{(model)}^{+0.003}_{-0.031}text{(parameterisation)}$ GeV and $m_b(m_b) = 4.049^{+0.104}_{-0.109}text{(exp/fit)}^{+0.090}_{-0.032}text{(model)}^{+0.001}_{-0.031} text{(parameterisation)}$~GeV.