For local observables $\{O_i(x_i)\}^n_{i = 1}$, one defines the Ward identity as
$$\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O(x_i)\rangle = \sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O(x_n)\rangle\tag{1}$$
Under the local variation of the field $$\Psi(x)\mapsto \Psi(x)+\epsilon\delta\Psi(x)\tag{2}$$
Whenever $\delta_{\epsilon}\mathcal{D}[\Psi(x)] = 0$, $\delta_{\epsilon} Z[\Psi(x)] = 0$, one obtains the following variation
$$\begin{align}\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle &= \delta_{\epsilon}\biggr(\frac{\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}}{Z[\Psi(x)]}\biggr)\\&= \frac{\delta_{\epsilon}\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)Z[\Psi(x)] - \delta_{\epsilon} Z[\Psi(x)]\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)}{Z^2[\Psi(x)]}\\&= \frac{1}{Z[\Psi(x)]}\delta_{\epsilon}\biggr(\int \mathcal{D}[\Psi(x)]\prod^n_{i = 1}O_i(x_i)e^{\frac{i}{\hbar}S[\Psi(x)]}\biggr)-\underbrace{\frac{1}{Z[\Psi(x)]}\delta_{\epsilon} Z[\Psi(x)]\langle \prod^n_{i = 1}O_i(x_i)\rangle}_{ = 0} \\ &= \sum^n_{i = 1}\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\langle\delta_{\epsilon} S[\Psi(x)]\prod^n_{i = 1}O_i(x_i)\rangle\\&= \sum^n_{i = 1}\langle O_1(x_1)\cdots\delta O_i(x_i)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\int_{M} d^n x\epsilon(x)\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle\end{align}\tag{3}$$
Where we assumed that $$Z[\Psi(x)] = \displaystyle\int D[\Psi(x)]e^{\frac{i}{\hbar}S[\Psi(x)]}\tag{4}$$
Further recalling that $$\int_{V}O_i(x)\delta(x-x_i)d^nx = \begin{cases}O_i(x_i), x_i\in V\\ 0, x_i\notin V\end{cases}\tag{5}$$
We observe that under arbitrary spacetime-dependent parameter $\epsilon : M\rightarrow \mathbb{R}$,
$$\begin{align}\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle &= \int_{M}d^nx\sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O_n(x_n)\rangle+\frac{i}{\hbar}\int_{M} d^nx\epsilon(x) \partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O_i(x_i)\rangle\end{align}\tag{6}$$
From which we conclude that whenever $\delta_{\epsilon}\langle \prod^n_{i = 1}O_i(x_i)\rangle = 0$,
$$-i\hbar\partial_{\mu}\langle j^{\mu}(x)\prod^n_{i = 1}O(x_i)\rangle = \sum^n_{i = 1}\delta(x-x_i)\langle O_1(x_1)\cdots\delta O_i(x)\cdots O(x_n)\rangle\tag{7}$$
However, I am not sure how the Dirac delta function $\delta(x-x_i)$ arises in the Ward identity in $(1)$. Can you provide clarification?