A generalization of the exponential variational expansion in relative coordinates to four-body systems is presented. It is shown that all required matrix elements can be expressed as derivatives upon the nonlinear parameters of the so-called auxiliary or basic four-body (six-dimensional) integral. If such an integral can be computed analytically, then it is possible to use the proposed variational approach for highly accurate variational calculations for various four-body systems with arbitrary particle masses without adiabatic or similar approximations. For instance, this approach can be used for highly accurate, bound state, variational calculations of the positronium molecule Ps(2), four-body muonic molecules, H-2, D-2, and DT molecules and many other systems. In fact, the present work opens a new avenue in the study of four-body systems.