Neutron stars harbour extremely powerful magnetic fields, leading to their shape being deformed. Their magnetic deformation depends both on the geometry -and strength -of their internal magnetic field and on their composition, encoded by the equation of state. However, both the details of the internal magnetic structure and the equation of state of the innermost part of neutron stars are mostly unkown. We performed a study of numerical models of magnetised, static, axisymmetric neutron stars in general relativity and in one of its most promising extensions, scalar-tensor theories. We did so by using several realistic equations of state currently allowed by observational and nuclear physics constraints, considering also those for strange quark stars. We show that it is possible to find simple relations among the magnetic deformation of a neutron star, its Komar mass, and its circumferential radius in the case of purely poloidal and purely toroidal magnetic configurations satisfying the equilibrium criterion in the Bernoulli formalism. These relations are quasi-universal, in the sense that they mostly do not depend on the equation of state. Our results, being formulated in terms of potentially observable quantities, could help to understand the magnetic properties of neutron stars interiors and the detectability of continuous gravitational waves by isolated neutron stars, independently of their equation of state. In the case of scalar-tensor theories, these relations depend also on the scalar charge of the neutron stars, thus potentially providing a new way to set constraints on the theory of gravitation.