That's a trivial simplification of the Johnson-Lindenstrauss lemma .
The matrix $X$ can be seen as a set of $q$ points in $p$ dimensions. Let us first prove the following result, than your question will follow directly with a union bound: (by the way, I'm going to assume that your matrix $G$ is normalized)
Theorem
Let $x \in \mathbb{R}^p$ and assume that the entries of the matrix $\tilde{G} \in \mathbb{R}^{s \times p}$ are drawn at random from a Gaussian distribution. Let $G = \frac{1}{\sqrt{s}}\tilde{G}$. Then $$ \mathbb{P}\left( |\|Gx\|_2^2 - \|x\|_2^2| \leq \varepsilon \|x\|_2^2 \right) \geq 1-2e^{-(\varepsilon^2-\varepsilon^3)s/4} $$
The proof uses a small lemma about the $\chi^2$ distribution which can be proved via Markov's inequality, after applying an exponential.
Proof
First of all, remember that $\mathbb{E}[\|Gx\|_2^2] = \|x\|_2^2$.
It is also worth mentioning that the individual components of the image vector $(Gx)_j/\|x\|_2$ for $1 \leq j \leq s$ are independent Normal variable with variance $1/\sqrt{s}$. As a consequence $\|Gx\|_2^2/\|x\|_2^2$ behaves as a $\chi^2/s$ random variable with $s$ degrees of freedom.
Thus it follows, using classical results on the $\chi^2$ distribution, that $$ \mathbb{P}\left(\|Gx\|_2^2 \geq (1+\varepsilon) \|x\|_2^2 \right) = \mathbb{P}(\sum_{j=1}^s \frac{(Gx)_j}{\|x\|_2} \geq 1 + \varepsilon) = \mathbb{P}(\frac{1}{s}\chi^2 \geq 1 + \varepsilon) = \mathbb{P}(\chi^2 \geq (1+\varepsilon)s) \leq e^{-(\varepsilon^2-\varepsilon^3)s/4}. $$
Similarly, the following bounds hold
$$ \mathbb{P}\left(\|Gx\|_2^2 \leq (1-\varepsilon) \|x\|_2^2 \right) = \mathbb{P}(\chi^2 \leq (1-\varepsilon)s) \leq e^{-(\varepsilon^2-\varepsilon^3)s/4}. $$
Combining both bounds together yields the expected result.
The result you asked for can be obtained by considering the norm $\|X\|_2^2 = \sum_{j = 1}^s\|x_j\|_2^2$ and using a union bound, playing with the constants
EDIT: adding some information relative to the estimations with the $\chi^2$ random variables.
Let $x_1,\cdots,x_s$ be independent random Gaussian variables.
$P=\mathbb{P}(\chi_s^2\geq (1+\varepsilon)s) = \mathbb{P}(\sum_{i=1}^sx_i^2 \geq (1+\varepsilon)s) = \mathbb{P}(\operatorname{exp}\left(\lambda \sum_{i=1}^sx_i^2\right) \geq \operatorname{exp}\left(\lambda (1+\varepsilon)s \right))$.
From here, we have Markov's inequality which yields $$P\leq \frac{\mathbb{E}\left(\operatorname{exp}\left(\lambda \sum_{i=1}^sx_i^2\right)\right)}{ \operatorname{exp}\left(\lambda (1+\varepsilon)s \right)} = \frac{\left(\mathbb{E}\left(\operatorname{exp}\left(\lambda x_i^2\right)\right)\right)^s}{ \operatorname{exp}\left(\lambda (1+\varepsilon)s \right)} = \frac{\left(\frac{1}{1-2\lambda}\right)^{k/2}}{\operatorname{exp}\left(\lambda (1+\varepsilon)s \right)}.$$
This expression is true for all $\lambda$, and in particular for $\lambda = \varepsilon/(2(1+\varepsilon))$, which achieves the minimum. We then only use the fact that $1+\varepsilon \leq \operatorname{exp}(\varepsilon-(\varepsilon^2-\varepsilon^3)/2)$. The other bound is obtained in a similar way.