Harmonic oscillator in spherical coordinates

Indeed, the supersymmetric operators do not factorize the Hamiltonian of the three-dimensional harmonic oscillator, there is an additional term. See Creation and annihilation operators, symmetry and supersymmetry of the 3D isotropic harmonic oscillator, equation 16.

Supplementing @CarloBeenakker's answer: a different sort of "factorization" of the Laplacian can be done in terms of the Dirac operator, which in some regards better preserves the nature of the situation, but/and has scalar coefficients in an appropriate Clifford algebra, which in general is not commutative. E.g., on $\mathbb R$, it is $\mathbb C$, which is convenient. But already in two dimensions the relevant Clifford algebra is the Hamiltonian quaternions...

Depending what you want to do, this might be a useful direction to consider, in addition to the "scalarization" of immediate reduction to the spherically symmetric case.