For continuous-time systems, you obtain the frequency response by evaluating the transfer function $H(s)$ on the imaginary axis $s=j\omega$ (assuming stability). In discrete-time, you get the $\mathcal{Z}$-transform instead of the Laplace transform, and the imaginary axis is replaced by the unit circle. So the frequency response is obtained by evaluating the transfer function $H(z)$ for $z=e^{j\omega}$. Note that for discrete-time systems, $\omega$ is a normalized frequency (in radians):
$$\omega=\frac{2\pi f}{f_s}\tag{1}$$
where $f_s$ is the sampling frequency.
Coming to you example, a frequency of $8000$ Hz with a sampling frequency of $4000$ Hz is the same as DC (Eq. $(1)$ should make that clear). Computing the frequency response at any frequency means to determine the desired $\omega$ from Eq. $(1)$, and then evaluate $H(z)$ with $z=e^{j\omega}$. In the case of DC ($\omega=0$) and Nyquist ($\omega=\pi$) this becomes especially easy: just evaluate $H(1)$ for DC, and $H(-1)$ for Nyquist.
For most other values of $\omega$, the frequency response is most likely complex-valued, so after computing $H(e^{j\omega})$ you just compute the magnitude of that complex number to obtain the gain at the given frequency. You don't need to find a general expression for $|H(z)|$.
