We describe the irrotational dust component of the universe in terms of a relativistic gradient expansion and transform the resulting synchronous metric to a Newtonian coordinate system. The two metrics are connected via a space-like displacement field and a time-like perturbation, providing a relativistic generalization of the transformation from Lagrangian to Eulerian coordinates. The relativistic part of the displacement field generates already at initial time a non-local density perturbation at second order. This is a purely relativistic effect since it originates from space-time mixing. We give two options, the passive and the active approach, on how to include the relativistic corrections for example in N-body simulations. In the passive approach we treat the corrections as a non-Gaussian modification of the initial Gaussian field (primordial non-Gaussianity could be incorporated as well). The induced non-Gaussianity depends on scale and the redshift at which initial conditions are set, with f_NL ~ few for small enough scales and redshifts. In the active approach we show how to use the relativistic trajectory to obtain the initial displacement and velocity of particles for N-body simulations without modifying the initial Gaussian field.