For simplicity of notation I'll leave out the index $k$. You've already noted that the error is
$$E=R\cdot G-S=S(HG-1)+NG\tag{1}$$
The MSE can be computed by integrating the power spectrum of the error. Since the power spectrum is real-valued and non-negative, its integral (i.e., the MSE) can be minimized by making the power spectrum as small as possible. Assuming white data and white noise with powers $\sigma_s^2$ and $\sigma_n^2$, respectively, which are independent of each other, we can compute the error power spectrum from (1) as
$$S_E=\sigma_s^2|HG-1|^2+\sigma_n^2|G|^2\tag{2}$$
The expression for the equalizer response $G$ that minimizes (2) can be found in an elegant way by completing the square in (2):
$$S_E=(\sigma_s^2|H|^2+\sigma_n^2)\left|G-\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\right|^2+\frac{\sigma_n^2}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{3}$$
The expression (3) can be minimized by choosing $G$ such that the left-hand term vanishes:
$$G_{opt}=\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{4}$$
which is the optimal MMSE solution.
It is not immediately obvious that (2) and (3) are identical, but proving their equivalence is simple (yet tedious) by multiplying both expressions out and comparing them.