Integrated Ward Identity Suppose you have the following ward identity : $$\int_{M} d^4x\ \epsilon(x)\ \partial_{\mu} \langle j_{\mu}(x)O(y)\rangle = - \ \langle\delta O(y)\rangle$$ where $\delta O(y)$ can be written in the following way to get the usual local form of the ward identity: $$\delta O(y) = \int_M d^4x \ \delta(x-y)\ \epsilon(x)\ \delta O(y)$$
Now suppose we choose $\epsilon$ (which is arbitrary) such that it's costant in M. Is it possible to use Stokes theorem and get :  $$\int_{\partial M} d\sigma_\mu\ \langle j_{\mu}(x)O(y)\rangle = - \ \langle\delta O(y)\rangle $$ or are we neglecting contact terms given by the $T$ product by doing so? Moreover is it possible to do the spectral decomposition of $\langle j_{\mu}(x)O(y)\rangle $ even if it's integrated?
 A: As it was noted in the comments, the derivatives $\partial_\mu$ is outside the correlator in the Ward identity so that the Stokes theorem trivially applies.
As for the second question, one can perform the spectral decomposition of $\langle j^\mu(x)\mathcal{O}(0)\rangle$ getting something along the lines of
$$
\langle j^\mu(x)\mathcal{O}(0)\rangle\sim \int \frac{d^4p}{(2\pi)^4} e^{ipx} p^\mu\int_0^\infty d\mu^2\frac{\rho(\mu^2)}{p^2-\mu^2+i\epsilon}
$$
that plugged into the integrated (say over all spacetime) Ward identity  gives
$$
\int d^4p \delta^4(p) p^2 \int d\mu^2 \frac{\rho(\mu^2)}{p^2-\mu^2+i\epsilon}
$$
which is clearly zero unless the integral over $d\mu^2$ doesn't generate a $1/p^2$ pole, that is a massless particle known as the Goldstone boson.
In fact, if the right-hand side isn't zero, the spectral function must contain such a Goldstone pole, 
$$
\rho(\mu^2)\sim \langle\mathcal{O}(0)\rangle\delta(\mu^2)+\sigma(\mu^2)
$$
so that 
$$
\int d\mu^2 \frac{\rho(\mu^2)}{p^2-\mu^2+i\epsilon}=-\frac{\langle\mathcal{O}(0)\rangle}{p^2+i\epsilon}\,.
$$
This is in fact another way to prove the Goldstone theorem.
PS: the $\sim$ symbol in the equations above means up to O(1) factors and signs which I haven't bothered to keep track of
A: You are describing precisely the route one follows to obtain the transformation laws of local operators from the Ward identity. I'll review the argument for you here in some detail, and I think that may resolve your confusion. 
Consider a quantum field theory with a continuous global symmetry under which local operators $\left\{\mathcal{O}_i \right\}$ transform as $\mathcal{O}_i \to \mathcal{O}_i + \epsilon\, \delta \mathcal{O}_i$, where $\epsilon$ is the parameter for the symmetry. By definition, there exists a local current operator $j_\mu(x)$ which satisfies the operator equation $\partial^\mu j_\mu (x) = 0$. This is the quantum version of Noether's theorem. By operator equation, I mean that $\partial^\mu j_\mu(x)$ vanishes inside of correlation functions with any other operators, provided that the other operators are not inserted at the same point as $j_\mu$. That is, $\left< \partial^\mu j_\mu(x) \mathcal{O}_1(x_1)\cdots \mathcal{O}_n(x_n) \right> = 0$ so long as $x_i \neq x$ for every $i\in\left\{1,\ldots,n\right\}$. 
More generally, when the operator insertions do approach the current insertion point, one picks up contact terms proportional to $\delta(x_i - x)$. Then Noether's theorem is generalized to the Ward identity (which is again part of the definition of a symmetry):
$$\left< \partial^\mu j_\mu(x) \mathcal{O}_1(x_1)\cdots \mathcal{O}_n(x_n) \right>=\sum_{i=1}^n \delta(x-x_i) \left<\mathcal{O}_1(x_1) \cdots \delta \mathcal{O}_i(x_i)\cdots \mathcal{O}_n(x_n)\right>.$$
Suppose we're studying a quantum field theory in $d$-dimensional flat spacetime. Let's integrate the Ward identity with respect to $x$ over a $d$-dimensional ball $B$ which contains the point $x_1$, but none of the other $x_i$:
$$\left< \left(\int_B \mathrm{d}\star j \right)\,\mathcal{O}_1(x_1)\cdots\right>=\left<\delta\mathcal{O}_1(x_1) \cdots \right>.$$
(I'm using differential form notation $\partial^\mu j_\mu \, \mathrm{d}^dx = \mathrm{d} \star j$, where $j = j_\mu \mathrm{d} x^\mu$ is the 1-form current). 
The ball $B$ is bounded by a $(d-1)$-dimensional surface $\Sigma = \partial B$. Applying Stokes' theorem, we find
$$\left< Q_\Sigma \,\mathcal{O}_1(x_1)\cdots\right>=\left<\delta\mathcal{O}_1(x_1) \cdots \right>,$$
where $Q_\Sigma := \int_\Sigma \star j$ is by definition the charge contained inside $\Sigma$.
Now, $Q_\Sigma$ is called a topological surface operator because I can freely deform the surface $\Sigma$ by moving it around or distorting it without affecting correlation functions of $Q_\Sigma$, provided that $\Sigma$ always surrounds the point $x_1$ and I don't drag it across any other operator insertion points $x_i$. To see why, consider the charge $Q_{\Sigma'}$ associated with a slightly deformed surface $\Sigma'$. The charges differ by
$$Q_\Sigma - Q_{\Sigma'} = \int_{\Sigma-\Sigma'} \star j.$$
$\Sigma$ and $\Sigma'$ are codimension 1 surfaces in spacetime, and their difference defines the boundary of a codimension 0 annulus $V$, with $\partial V = \Sigma-\Sigma'$. Thus,
$$Q_\Sigma - Q_{\Sigma'} = \int_{\partial V} \star j = \int_V \mathrm{d} \star j.$$
But as we just discussed, $\mathrm{d} \star j = \partial^\mu j_\mu \, \mathrm{d}^d x$ vanishes inside correlation functions up to contact terms. So as long as I don't cross an operator insertion in deforming $\Sigma$ to $\Sigma'$, the charge $Q_\Sigma$ is independent of deformations of $\Sigma$. 
Thus, integrating the Ward identity has led to the operator equation
$$Q_\Sigma \mathcal{O}(x) = \delta \mathcal{O}(x),$$
where $\Sigma$ is a codimension 1 surface surrounding the point $x$.
Now typically we evaluate charges by integrating over spatial slices. So let us exploit the above topological invariance to make the "ball" $B$ a thin box obtained by dragging a spatial slice $\Sigma^-$ at time $x^0 - \Delta$ through $x$ and ending with another spatial slice $\Sigma^+$ at time $x^0 + \Delta$, i.e. $\Sigma = \partial B = \Sigma^+ - \Sigma^-$.  Then the integrated Ward identity becomes
$$Q_{\Sigma^+} \mathcal{O}(x) - Q_{\Sigma^-}\mathcal{O}(x) = \delta \mathcal{O}(x).$$ 
Finally, we interpret this operator equation as a time-ordered correlation function:
$$\left< Q_{\Sigma^+} \mathcal{O}(x) - \mathcal{O}(x) Q_{\Sigma^-} \right> = \left< \delta \mathcal{O}(x) \right>.$$
And in the limit $\Delta \to 0$ we find
$$\left<[Q,\mathcal{O}(x)]\right> = \left< \delta\mathcal{O}(x)\right>,$$
where $Q$ is the charge associated to the spatial slice through $x^0$. This is of course the familiar form of the  transformation law for a local operator charged under a symmetry generator $Q$. 
