Yes, actually you made a calculation mistake. Indeed, careful calculation shows that
$\nabla \cdot \mathbf{F} = 0$
away from the origin. To be specific, let us rename the coordinate variables as $(x_1, x_2, x_3) = (x, y, z)$. Then
$\begin{align*}\nabla \cdot \mathbf{F} &= \sum_{i=1}^{3} \frac{\partial}{\partial x_i}\frac{x_i}{r^3} = \sum_{i=1}^{3} \frac{\partial x_i}{\partial x_i}\frac{1}{r^3} + x_i \frac{\partial}{\partial x_i}\frac{1}{r^3} \\ &= \sum_{i=1}^{3} \frac{1}{r^3} + x_i \left(-\frac{3}{r^4}\right)\frac{\partial}{\partial x_i}\sqrt{x_1^2+x_2^2+x_3^2} \\ &= \frac{3}{r^3} - \frac{3}{r^4} \sum_{i=1}^{3} x_i \cdot \frac{x_i}{r} \\ &= \frac{3}{r^3} - \frac{3}{r^3}=0. \end{align*}$
Now you may apply the divergence theorem to obtain
$\int_{S}\mathbf{F}\cdot\mathbf{n}\;da = \int_{V}\nabla \cdot \mathbf{F}\;dv = 0,$
where $V$ is the region inside the sphere. This calculation is valid since the origin does not lie on and inside the sphere $S$.
If you made a bold and brave resolution to attack this problem by brutal integration, then you may do as follows: First, make the substitution $(x',y',z')=(x-2,y,z)$. This is equivalent to shifting the whole situation so that the center of $S$ coincides with the origin. (This makes the domain of integration simpler in expense of the simplicity of the integrand.) Then make a variant of spherical coordinate change as follows:
$ \begin{align*} x' &= r \cos\phi \\ y' &= r \sin\phi\cos\theta\\ z' &= r \sin\phi\sin\theta. \end{align*}$
Of course, $0 \leq \phi \leq \pi$ and $0 \leq \theta < 2\pi$. Then on $S$, with this new coordinate system, $S$ is characterizes as the locus of $r = 1$, and on $S$ we have $da = r^2\sin\phi\,d\phi d\theta = \sin\phi\,d\phi d\theta$. Thus our integral becomes
$\begin{align*}\int_S \mathbf{F}\cdot\mathbf{n}\;da &=\int_{0}^{2\pi}\int_{0}^{\pi} \frac{(x'+2,y',z')}{((x'+2)^2+y'^2+z'^2)^{3/2}} \cdot \frac{(x',y',z')}{(x'^2+y'^2+z'^2)^{1/2}}\;r^2 \sin\phi\,d\phi d\theta\\ &=\int_{0}^{2\pi}\int_{0}^{\pi} \frac{(x'+2,y',z') \cdot (x',y',z')}{(r^2 + 4r\cos\phi + 4)^{3/2}}\;r \sin\phi\,d\phi d\theta\\ &=\int_{0}^{2\pi}\int_{0}^{\pi} \frac{r (r^2 + 2r\cos\phi) \sin\phi}{(r^2 + 4r\cos\phi + 4)^{3/2} }\;d\phi d\theta\\ &=\int_{0}^{2\pi}\int_{0}^{\pi} \frac{(1 + 2\cos\phi) \sin\phi}{(5 + 4\cos\phi)^{3/2} }\;d\phi d\theta\\ &=2\pi \int_{0}^{\pi} \frac{(1 + 2\cos\phi) \sin\phi}{(5 + 4\cos\phi)^{3/2} }\;d\phi \end{align*}$
Finally, let $t = 4\cos\phi$. Then $dt = -4\sin\phi \, d\phi$ and therefore
$\begin{align*} \int_S \mathbf{F}\cdot\mathbf{n}\;da &= \frac{\pi}{4} \int_{-4}^{4} \frac{2 + t}{(5 + t)^{3/2} }\;dt \\ &= \frac{\pi}{4} \left[ -2 \cdot \frac{2 + t}{\sqrt{5 + t}}\right]_{-4}^{4} + \frac{\pi}{2} \int_{-4}^{4} \frac{1}{\sqrt{5-t}}\;dt \\ &= -2\pi + \frac{\pi}{2} \left[2\sqrt{5+t}\right]_{-4}^{4} = 0. \end{align*}$