My formula in the comment for $\|z\|=1$, was obtained as follows.

First note that the integral is well defined in $[0,\infty]$, and by monotone convergence, is given by
$$
\int_{S^{2n-1}}\frac{1}{|1-\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
=\lim_{r\rightarrow 1^{-}}
\int_{S^{2n-1}}\frac{1}{|1-r\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
$$
Now for $r\in(0,1)$, one can use the convergent representation
$$
\frac{1}{1-r\langle z,\zeta\rangle}=\sum_{k\ge 0} r^k \langle z,\zeta\rangle^k
$$
and its conjugate, in order to write
$$
\int_{S^{2n-1}}\frac{1}{|1-r\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
=\int_{S^{2n-1}}
\sum_{k,\ell\ge 0} r^{k+\ell}
\langle z,\zeta\rangle^k \langle \zeta,z\rangle^{\ell}
\ d\sigma(\zeta)
$$
$$
=\sum_{k,\ell\ge 0} r^{k+\ell}
\int_{S^{2n-1}}
\langle z,\zeta\rangle^k \langle \zeta,z\rangle^{\ell}
\ d\sigma(\zeta)
$$
The analogue of the Isserlis-Wick formula for integrating monomials on the complex unit sphere is
$$
\int_{S^{2n-1}} \zeta_{i_1}\cdots \zeta_{i_k}\ \overline{\zeta_{j_1}}\cdots\overline{\zeta_{j_{\ell}}}\ \ d\sigma(\zeta)
=\delta_{k,\ell}\ \frac{k!}{n(n+1)\cdots(n+k-1)}\ S_{i_1,\ldots,i_k}^{j_1,\ldots,j_k}
$$
where
$$
S_{i_1,\ldots,i_k}^{j_1,\ldots,j_k}:=\frac{1}{k!}\ \sum_{\sigma\in\mathfrak{S}_k}\delta_{i_1,j_{\sigma(1)}}\cdots\delta_{i_k,j_{\sigma(k)}}
$$
is a symmetrizer. The formula holds for all assignments of the $i$ and $j$ indices in $\{1,2,\ldots,n\}$.
This is easily proved from the Isserlis-Wick theorem for the complex Gaussian
$$
\frac{1}{\pi^n}\int_{\mathbb{C}^n}
e^{-|\langle \zeta,\zeta\rangle|^2}
\ \zeta_{i_1}\cdots \zeta_{i_k}\ \overline{\zeta_{j_1}}\cdots\overline{\zeta_{j_{\ell}}}\ \ \prod_{a=1}^{n}d({\rm Re}\ \zeta_a)\ d({\rm Im}\ \zeta_a)
$$
$$
=\delta_{k,\ell}\sum_{\sigma\in\mathfrak{S}_k}\delta_{i_1,j_{\sigma(1)}}\cdots\delta_{i_k,j_{\sigma(k)}}\ .
$$
by separating the radial part from the spherical part, using the homogeneity of the integrand.

As a result,
$$
\int_{S^{2n-1}}
\langle z,\zeta\rangle^k \langle \zeta,z\rangle^{\ell}
\ d\sigma(\zeta)
$$
vanishes when $k\neq \ell$, and is equal to
$$
\frac{k!}{n(n+1)\cdots(n+k-1)}=\frac{1}{\binom{n+k-1}{k}}=\frac{(n-1)!}{(k+1)(k+2)\cdots(k+n-1)}
$$
when $k=\ell$.
We then get the convergent series representation
$$
\int_{S^{2n-1}}\frac{1}{|1-r\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
=(n-1)!\ \sum_{k=0}^{\infty}\frac{r^{2k}}{(k+1)(k+2)\cdots(k+n-1)}
$$
After taking the limit $r\rightarrow 1$, we get
$$
\int_{S^{2n-1}}\frac{1}{|1-\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
=(n-1)!\ \sum_{k=0}^{\infty}\frac{1}{(k+1)(k+2)\cdots(k+n-1)}
$$
which diverges for $n=1,2$, and otherwise is a telescopic sum, given the identity
$$
\frac{1}{(k+1)(k+2)\cdots(k+n-1)}=\frac{1}{(n-2)}
\left[
\frac{1}{(k+1)(k+2)\cdots(k+n-2)}
-\frac{1}{(k+2)(k+3)\cdots(k+n-1)}
\right]\ .
$$
Finally, for $n\ge 3$,
$$
\int_{S^{2n-1}}\frac{1}{|1-\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)=\frac{n-1}{n-2}\ .
$$

Now going back to the original question of computing the integral when $z$ is in the open unit ball, rather than the sphere (as I misunderstood), one has, by the above computation,
$$
\int_{S^{2n-1}}\frac{1}{|1-\langle z,\zeta\rangle|^2}\ d\sigma(\zeta)
=(n-1)!\
\sum_{k=0}^{\infty}\frac{r^{2k}}{(k+1)(k+2)\cdots(k+n-1)}\ ,
$$
when $\|z\|=r\in (0,1)$.
With the telescoping sum identity above, one can then put this in the form mentioned by Peter. I don't know if one has completely explicit formulas for the result which essentially is a hypergeomtric series
$$
{}_2 F_1\left[\begin{array}{c}1,1\\ n\end{array}; r^2\right]\ .
$$
Note that Peter's answer is of course simpler than mine, but I wanted to advertise the Isserlis-Wick approach which is much more powerful than say Funk-Hecke, because it can also be used for integrating over several vectors $\zeta$. It also shows that what looks like a problem of multivariate calculus is really a problem of combinatorics. See these two other MO questions for similar examples:

Integration of a function over 7-sphere

Moments of Plücker coordinates on complex Grassmannian