Suppose that $f$ is a weakly holomorphic modular form of weight $k$ on $\Gamma_0(N)$, and let $p$ be prime with $p\nmid N$. Further suppose that the $q$-expansion of $f$ at the cusp $\infty$ has $p$-integral coefficients. Is it true that $f|_kW_e$ also has $p$-integral coefficients where $$W_e = \begin{pmatrix} ae & b \\ Nc & de\end{pmatrix} \quad a,b,c,d\in \mathbb Z\;\text{ and }\det W_e = e$$ is the Atkin-Lehner involution for $e$, an exact divisor of $N$?
The reason I am hoping this is true is that I would like to be able to apply $W_e$ to both sides of a congruence of the form $f \equiv g \mod{p}$ where $f$ and $g$ are weakly holomorphic modular forms.
If we only wanted $f|_kW_e$ to have rational coefficients, then this seems to be a standard fact coming from the modular curve $\Gamma_0(N)\backslash \mathbb H^*$ being defined over $\mathbb Q$. I suspect that the corresponding fact for $\mathbb Z_{(p)}$ follows from similar reasoning, but I have not been able to find a reference for it.
This follows from the fact that the (open) modular curve $Y_0(N) = \Gamma_0(N) \backslash \mathbb{H}$ has a model as a scheme over $\mathbf{Z}[1/N]$, and the action of $W_q$ for $q \| N$ extends to this model.
The existence of the model you can find in all sorts of places, e.g. Katz's article in the Antwerp proceedings. I don't remember offhand whether the fact that $W_q$ extends to the model is written down explicitly in any of the standard references but it's an easy computation.