Three-dimensional direct numerical simulations are used to study the energy cascade rate in isothermal compressible magnetohydrodynamic turbulence. Our analysis is guided by a two-point exact law derived recently for this problem in which flux, source, hybrid and mixed terms are present. The relative importance of each term is studied for different initial subsonic Mach numbers $M_{S}$ and different magnetic guide fields $\boldsymbol{B}_{0}$ . The dominant contribution to the energy cascade rate comes from the compressible flux, which depends weakly on the magnetic guide field $\boldsymbol{B}_{0}$ , unlike the other terms whose moduli increase significantly with $M_{S}$ and $\boldsymbol{B}_{0}$ . In particular, for strong $\boldsymbol{B}_{0}$ the source and hybrid terms are dominant at small scales with almost the same amplitude but with a different sign. A statistical analysis undertaken with an isotropic decomposition based on the SO(3) rotation group is shown to generate spurious results in the presence of $\boldsymbol{B}_{0}$ , when compared with an axisymmetric decomposition better suited to the geometry of the problem. Our numerical results are compared with previous analyses made with in situ measurements in the solar wind and the terrestrial magnetosheath.