If your question is whether it is possible to explicitly evaluate the integral in the right-hand side of your equation when $u(x) = (1 + |x|^2)^{-p}$, then I do not know the answer, but I would not be surprised if it be affirmative. The closest things that comes to my mind are:
The explicit expression for $(-\Delta)^s u$ when $p = \tfrac{N + 1}{2}$ or $p = \tfrac{N - s}{2} + n$ ($n = 1, 2, \ldots$), found by S. Samko (see [1]).
The explicit expression for $(-\Delta)^s (1 - |x|^2)_+^{-p}$, found independently by B. Dyda [2] and P. Biller, C. Imbert, G. Karch [3].
The methods for doing integrals of this form, involving Kelvin transformation, were already developed by M. Riesz.
If, however, you are looking for an arbitrary explicit way to evaluate $\|(-\Delta)^{s/2} u\|_2 = \langle (-\Delta)^s u, u \rangle$, then there are other approaches. The most natural thing to do would be to find the Fourier transform of $u$, and use Plancherel's theorem. Since $u(x) = (1 + |x|^2)^{-p}$ is a radial function, its Fourier transform is given by the Hankel transform of the profile $(1 + r^2)^{-p}$. I did not attempt to search for the explicit expression, but I bet it is given in one of the standard tables of integrals.
Yet another approach is to use Mellin transform (rather than Fourier transform) to find an explicit expression for $(-\Delta)^s u$. For $u(x) = (1 + |x|^2)^{-p}$ this is indeed possible, and the result involves the hypergeometric function $_2F_1$; see Corollary 2 (or Corollary 1) in my joint paper with B. Dyda and A. Kuznetsov [4]. The last step would be to evaluate the inner product of $(-\Delta)^s u$ with $u$ (or the $L^2$ norm of $(-\Delta)^{s/2} u$, whichever turns out simpler); again, this is likely to be found in standard integral tables.
EDIT: One more thought: if one is able to evaluate the convolution of $u$ with the Gauss–Weierstrass kernel, or, even better, the value of $$q_t = \int_{\mathbb{R}^N} \int_{\mathbb{R}^N} (u(x) - u(y))^2 (4 \pi t)^{-N/2} \exp(-|x - y|^2 / (4 t)) dx dy,$$ then it can be convenient to use the subordination formula: $$\|(-\Delta)^{s/2} u\|_2^2 = \frac{1}{|\Gamma(-s)|} \int_0^\infty q_t t^{-1 - s} dt.$$
EDIT: As suggested by leo monsaingeon, if the convolution of $u$ with $(y^2 + |\cdot|^2)^{-(N + 2 s)/2}$ is known, then $(-\Delta)^s u$ and $\|(-\Delta)^{s/2} u\|_2$ can be evaluated using the Caffarelli–Silvestre extension technique; see [5].
(One could also consider, for example, Balakhrishnan's formula, but I doubt this is ever useful in calculations: the resolvent kernel for $\Delta$ is not the simplest convolution factor).
References:
[1] Samko, S.: Hypersingular Integrals and Their Applications, Analytical Methods and Special Functions, vol. 5. Taylor & Francis, Ltd., London (2002)
[2] Dyda, B.: Fractional calculus for power functions and eigenvalues of the fractional Laplacian. Fract. Calc. Appl. Anal. 15(4), 536–555 (2012)
[3] Biler, P., Imbert, C., Karch, G.: Barenblatt profiles for a nonlocal porous medium equation. C. R. Math. Acad. Sci. Paris 349(11–12), 641–645 (2011)
[4] Dyda, B., Kuznetsov, A., Kwaśnicki, M.: Fractional Laplace operator and Meijer G-function. Constructive Approx. 45(3), 427–448 (2017)
[5] Caffarelli, L., Silvestre, L.: An Extension Problem Related to the Fractional Laplacian. Comm. PDE 32(8), 1245–1260 (2007)
