Why does the divergence in my integral has a change of type?

I am looking at the following integral in Euclidean space:

$$I = \int_{\mathbb{R}} d\tau_3 \int_{\mathbb{R}} d\tau_4 \int_{\mathbb{R}^4} d^4 x_5 \frac{1}{x_{15}^2 x_{25}^2 x_{35}^2 x_{45}^2} \tag{1}$$

with $$x_{ij} := x_i - x_j$$, $$x_1=(1,0,0,0)$$, $$x_2=(x_2^1,x_2^2,0,0)$$, $$x_3=(0,0,0,\tau_3)$$ and finally $$x_4=(0,0,0,\tau_4)$$. This integral is divergent at $$x_5 \sim \tau_3 \sim \tau_4$$.

Using dimensional regularization $$\left(\int d^4 x_5 \to \int d^{2\omega} x_5 \right)$$, I could extract the divergence of $$I$$ by setting $$x_{15} x_{25} \sim x_{13} x_{24}$$, integrating the remaining integrals and expanding the resulting $$\Gamma$$-function. The result is:

$$I_\text{div} = \pi^4 \frac{1}{\left| x_1 \right| \left| x_2 \right|} \frac{1}{\epsilon} \tag{2}$$

with $$\epsilon \to 0$$ as $$2\omega \to 4$$.

But I am interested in the full result of the integral, i.e. divergence + finite part. My idea now is to compute the integral given in $$(1)$$ as far as I can in $$4$$d, and to try to "push" the divergence into the last one-dimensional integral. For example, I can start by performing the integrals over the $$\tau$$'s:

$$\int_{\mathbb{R}} d\tau_3 \frac{1}{x_{35}^2} = \int_{\mathbb{R}} d\tau_3 \frac{1}{\left[ (\tau_3 - \tau_5)^2 + \vec{x}_5^2 \right]} = \frac{\pi}{\left| \vec{x}_5 \right|}, \tag{3}$$

where I used the fact that $$\vec{x}_3 = (0,0,0)$$. Regarding the notation, I used here $$\vec{x}^2 := \vec{x} \cdot \vec{x}$$. After that, I can also perform the integral over $$\tau_5$$:

\begin{align}I &= \int_{\mathbb{R}^3} d^3 x_5 \frac{1}{\vec{x}_5^2} \int_{\mathbb{R}} d\tau_5 \frac{1}{\left(\tau_5^2 + \vec{x}_{15}^2\right)\left(\tau_5^2 + \vec{x}_{15}^2\right)} \\ &= \int_{\mathbb{R}^3} d^3 x_5 \frac{1}{\vec{x}_5^2} \frac{1}{\left| \vec{x}_{15} \right| \vec{x}_{25}^2 + \left| \vec{x}_{25} \right| \vec{x}_{15}^2}. \tag{4} \end{align}

Going further integrating dimension by dimension, it seems that I can actually go pretty far, maybe all the way up to the last dimension. However the computations become quite involved, and I would like to have a way to check my result in between. For example, since I know that my divergence is logarithmic, I expect to be able to extract the divergence again and get something like:

$$I_\text{div} \sim \pi^4 \left| x_1 \right| \left| x_2 \right| \int_0^\infty \frac{dr}{r}. \tag{5}$$

So I tried to do that with eq. ($$5$$). The divergence seems to sit at $$\vec{x}_5 = (0,0,0)$$, hence I try to extract it by setting $$\left| \vec{x}_{15} \right| \vec{x}_{25}^2 + \left| \vec{x}_{25} \right| \vec{x}_{15}^2 \sim \left| x_1 \right| x_2^2 + \left| x_2 \right| x_1^2$$, and we have:

\begin{align} I_\text{div} &= \pi^3 \frac{1}{\left| x_1 \right| x_2^2 + \left| x_2 \right| x_1^2} \int_{\mathbb{R}^3} d^3 x \frac{1}{\vec{x}^2} \\ &= 4 \pi^4 \tag{6} \frac{1}{\left| x_1 \right| x_2^2 + \left| x_2 \right| x_1^2} \int_0^\infty dr \end{align}

But I am confused: the prefactor with the $$x_1$$, $$x_2$$ is not the same as in eq. ($$2$$), and moreover the divergence is now linear! Am I not allowed to do that? And if not, why?

I have heard that dimensional regularization "sees" only the logarithmic divergences, maybe that could be the reason?