Note that the argument of the cosine is a dot product between two $n$-dimensional vectors. It depends on the angle between $x$ and $y$, so we cannot simply replace $x$ by the radial coordinate $r$ when switching to polar coordinates. However, for a given $y \in \mathbb{R}^n \setminus \{0\}$ we can always rotate the $x$-coordinate system in such a way that $y$ points along the $x_1$-axis. Then $x \cdot y = \lvert y \rvert x_1$ holds and we obtain
$$ I (y) = \int \limits_{\mathbb{R}^n} \frac{1 - \cos(\lvert y \rvert x_1)}{\lvert x \rvert^{n+q}} \, \mathrm{d} x \, ,$$
since the Lebesgue measure is invariant under rotations. Introducing $z = \lvert y \rvert x$, we find the desired result
$$ I(y) = \int \limits_{\mathbb{R}^n} \frac{1 - \cos(z_1)}{\lvert z \rvert^{n+q}} \, \mathrm{d} z \, \frac{\lvert y \rvert^{n+q}}{\lvert y \rvert^n} \equiv c \lvert y \rvert^q \, .$$
Clearly, we have
$$ c = \int \limits_{\mathbb{R}^n} \frac{1 - \cos(z_1)}{\lvert z \rvert^{n+q}} \, \mathrm{d} z > 0 \, ,$$
since the integrand is positive almost everywhere. $c < \infty$ follows as in GReyes' answer and we can even give an explicit estimate:
\begin{align}
c &= \int \limits_{B_1(0)} \frac{1 - \cos(z_1)}{\lvert z \rvert^{n+q}} \, \mathrm{d} z + \int \limits_{\mathbb{R}^n \setminus B_1(0)} \frac{1 - \cos(z_1)}{\lvert z \rvert^{n+q}} \, \mathrm{d} z \leq \int \limits_{B_1(0)} \frac{z_1^2}{2\lvert z \rvert^{n+q}} \, \mathrm{d} z + \int \limits_{\mathbb{R}^n \setminus B_1(0)} \frac{2}{\lvert z \rvert^{n+q}} \, \mathrm{d} z \\
&\leq \frac{1}{2} \int \limits_{B_1(0)} \frac{1}{\lvert z \rvert^{n+q-2}} \, \mathrm{d} z + 2 \int \limits_{\mathbb{R}^n \setminus B_1(0)} \frac{1}{\lvert z \rvert^{n+q}} \, \mathrm{d} z = \operatorname{vol}(\mathbb{S}^{n-1}) \left[\frac{1}{2} \int \limits_0^1 r^{1-q} \, \mathrm{d} r + 2 \int \limits_1^\infty r^{-(1+q)} \, \mathrm{d} r\right] \\
&= \operatorname{vol}(\mathbb{S}^{n-1}) \left[\frac{1}{2(2-q)} + \frac{2}{q}\right] = \operatorname{vol}(\mathbb{S}^{n-1}) \frac{8-3q}{2 q (2-q)} < \infty \, .
\end{align}
In order to find the exact value of $c$ we need to introduce polar coordinates. Since $z_1 = r \cos(\varphi_1)$, the integrand depends on the angle $\varphi_1$ and only the integrations over $\varphi_2, \dots, \varphi_{n-1}$ are trivial. Therefore, we need to divide the prefactor $\operatorname{vol}(\mathbb{S}^{n-1})$ by the contribution of the $\varphi_1$-integral:
$$ c = \frac{\operatorname{vol}(\mathbb{S}^{n-1})}{\int_0^\pi \sin^{n-2}(\varphi_1) \, \mathrm{d} \varphi_1} \int \limits_0^\pi \int \limits_0^\infty \frac{1 - \cos(r \cos(\varphi_1))}{r^{n+q}} r^{n-1} \sin^{n-2} (\varphi_1) \, \mathrm{d} r \, \mathrm{d} \varphi_1 \, .$$
The integral in the denominator is a beta function integral. The remaining double integral can be simplified by exploiting the symmetry of the integrand about $\varphi_1 = \frac{\pi}{2}$ and integrating by parts with respect to $r$:
$$ c = \frac{\operatorname{vol}(\mathbb{S}^{n-1})}{\operatorname{B}\left(\frac{n-1}{2},\frac{1}{2}\right)} \frac{2}{q} \int \limits_0^{\pi/2} \int \limits_0^\infty \frac{\sin(r \cos(\varphi_1))}{r^q} \cos(\varphi_1) \sin^{n-2}(\varphi_1) \, \mathrm{d} r \, \mathrm{d} \varphi_1 \, .$$
Finally, the change of variables $u = r \cos(\varphi_1)$ leads to two independent integrals. One of them is another beta function integral and the other has been solved here (we can simplify the result by applying the gamma function reflection formula). We find
\begin{align}
c &= \frac{\operatorname{vol}(\mathbb{S}^{n-1})}{\operatorname{B}\left(\frac{n-1}{2},\frac{1}{2}\right)} \frac{2}{q} \int \limits_0^{\pi/2} \sin^{n-2}(\varphi_1) \cos^{q}(\varphi_1) \, \mathrm{d} \varphi_1 \int \limits_0^\infty \frac{\sin(u)}{u^q} \, \mathrm{d} u \\
&= \frac{\operatorname{vol}(\mathbb{S}^{n-1})}{\operatorname{B}\left(\frac{n-1}{2},\frac{1}{2}\right)} \frac{1}{q} \operatorname{B}\left(\frac{n-1}{2},\frac{q+1}{2}\right) \frac{\pi}{2 \operatorname{\Gamma}(q) \sin\left(\frac{\pi}{2} q\right)} \, .
\end{align}
After plugging in the known values for the surface area and the beta function we arrive at
$$ c = \frac{\pi^{\frac{n+1}{2}} \operatorname{\Gamma}\left(\frac{q+1}{2}\right)}{\sin\left(\frac{\pi}{2} q\right) \operatorname{\Gamma}(q + 1) \operatorname{\Gamma}\left(\frac{q+n}{2}\right)} \, . $$
In particular, for $q=1$ we have
$$ c = \frac{\pi^{\frac{n+1}{2}}}{\operatorname{\Gamma}\left(\frac{n+1}{2}\right)} \, . $$
