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 MS and different magnetic guide fields B0. The dominant contribution to the energy cascade rate comes from the compressible flux, which depends weakly on the magnetic guide field B0, unlike the other terms whose modulus increase significantly with MS and B0. In particular, for strong B0 the source and hybrid terms are dominant at small scales with almost the same amplitude but with a different sign. A statistical analysis made with an isotropic decomposition based on the SO(3) rotation group is shown to generate spurious results in presence of B0, 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.