The large scale structure bispectrum couples large scales (where general relativity is relevant) and short scales (where non-linearities are important) such that it contains information about the dynamics of the early universe in the squeezed limit. We start by performing a calculation which is non-linear and is based in General relativity under the weak field approximation with the purpose to solve the dynamics for dark matter fluid. The solutions for dark matter perturbations are used to compute relativistic corrections to correlation functions. Specifically, we compute the three-point matter correlation function (matter bispectrum) at one loop, we find its relativistic corrections to be as large as the Newtonian result in the squeeze limit. In this limit we also find the relativistic correction bispectrum to be degenerated with a primordial non-Gaussianity signal. We also compute matter power spectrum which has a small relativistic correction. Our next step is to connect dark matter perturbation with the number density of galaxies through a bias expansion. To begin with we write the Lagrangian bias expansion in terms of operators built on the curvature of early time hypersurface of comoving observer. We take this as initial condition and evolve it to obtain the Eulerian bias description in the general relativistic framework which is also renormalized properly. These are the first steps towards to determine the gauge invariant galaxy bispectrum.