I find anomalies an extremely unintuitive subject, because they're studied so indirectly. In the standard textbook presentation, one computes an abstract quantity that should be zero classically (say, some time-ordered correlation function) and finds that it isn't, and that's it.
I'm wondering how anomalies manifest at a basic level, say, at the level of the equations of motion. However, in this case I'm not sure how they can possibly appear.
Naively, upon canonical quantization, interacting quantum fields will obey equal-time canonical commutation relations, just as the corresponding classical fields have canonical Poisson brackets. Then the equations of motion for classical fields translate to operator equations for the fields, e.g. in $\phi^4$ theory we would have $$(\partial^2 + m^2) \hat{\phi} = \hat{\phi}^3$$ modulo formal issues involving multiplying operator-valued distributions.
Next, we may take the time-ordered vev of both sides. Pulling time derivatives past the time-ordering symbol produces the appropriate contact terms to arrive at the Schwinger-Dyson equations.
At the classical level, a conserved current obeys $\partial_\mu j^\mu = 0$. Hence by the same reasoning as for $\hat{\phi}$, this should hold as an operator equation, $\partial_\mu \hat{j}^\mu = 0$. Taking the time-ordered vev of both sides and pulling out the time derivative reproduces the Ward-Takahashi identity, with no anomaly terms.
I'm sure that almost everything I said above is extremely naive from the point of view of formal QFT, but morally speaking, where does the above logic go wrong? For the axial anomaly, the textbook arguments show that $$\partial_\mu \hat{j}^{\mu 5} \sim - \frac{g^2}{16 \pi^2} \epsilon^{\mu\nu\rho\sigma} \hat{F}_{\mu\nu} \hat{F}_{\rho \sigma}$$ where the $\sim$ indicates that only some matrix elements of the operators on both sides are shown to match. Is this in fact an operator equality? If so, where does the extra term slip in above?