We perform theoretical and numerical studies of the full relativistic two-point galaxy correlation function, considering the linear-order scalar and tensor perturbation contributions and the wide-angle effects. Using the gauge-invariant relativistic description of galaxy clustering and accounting for the contributions at the observer position, we demonstrate that the complete theoretical expression is devoid of any long-mode contributions from scalar or tensor perturbations and it lacks the infrared divergences in agreement with the equivalence principle. By showing that the gravitational potential contribution to the correlation function converges in the infrared, our study justifies an IR cut-off (kIR ≤ H0) in computing the gravitational potential contribution. Using the full gauge-invariant expression, we numerically compute the galaxy two-point correlation function and study the individual contributions in the conformal Newtonian gauge. We find that the terms at the observer position such as the coordinate lapses and the observer velocity (missing in the standard formalism) dominate over the other relativistic contributions in the conformal Newtonian gauge such as the source velocity, the gravitational potential, the integrated Sachs-Wolf effect, the Shapiro time-delay and the lensing convergence. Compared to the standard Newtonian theoretical predictions that consider only the density fluctuation and redshift-space distortions, the relativistic effects in galaxy clustering result in a few percent-level systematic errors beyond the scale of the baryonic acoustic oscillation (~ 2% at 150 Mpc/h and redshift one). Our theoretical and numerical study provides a comprehensive understanding of the relativistic effects in the galaxy two-point correlation function, as it proves the validity of the theoretical prediction and accounts for effects that are often neglected in its numerical evaluation.