In your derivation you've implicitly assumed that the wavefunction does not change its values when you go to the Galilean boosted frame. In other words, you've assumed $\psi'(x', t') = \psi \bigl(x(x',t'), t(x',t') \bigr)$. However, this isn't right.
The wavefunction encodes information about a particle's momentum, so when you go to a different frame the wavefunction must change to represent the momentum the particle has in the new frame. For example, in the case of a plane wave, $\psi(x, t) = e^{i(kx - \omega t)}$. When you boost by velocity $u$, the wave's $k, \omega$ must change to match the new momentum and energy, like $k' = k + mu/\hbar$ and $\omega' = \hbar k'^2 / 2m$. So it's not the same function anymore; it has different wavelength and frequency, above and beyond the simple change of coordinates.
In other words, the Schrödinger equation is Galilean-invariant not in the sense that the same solution works after a boost, but in the sense that there exist solutions representing waves traveling at all different speeds. A boost maps a solution with $k, \omega$ to another solution with $k', \omega'$. (And for solutions that are not plane waves, we can use the Fourier transform to decompose them into plane waves, boost each one, and recompose them.)
This argument might be a bit unsatisfying since it relies on physical intuition about the meaning of wavelength and frequency, rather than being a purely mathematical derivation. I wonder if you might be able to derive it more rigorously from some operator-algebraic considerations, or some such.
The problem is also solved if $\vec\nabla A$ is non-zero. Your last line equals
$$ \mathrm e^{\mathrm i\alpha} \left( -\Delta \psi + \mathrm i q\, (\vec\nabla \cdot \vec A)\, \psi + 2 \mathrm i q\, \vec A \cdot \vec\nabla \psi + (q\vec A)^2\, \psi\right) \;. $$
Note that this is the same as
$$ \mathrm e^{\mathrm i\alpha}\, ( \vec p - q\vec A ) ( \vec p - q\vec A )\, \psi \;. $$
Edit to respond to comment.
You already did everything correctly in your question, $\vec p - q\vec A$ acts as follows on a wave function $\psi(x)$:
$$ (( \vec p - q\vec A ) \psi)(x) = -\mathrm i\hbar\, \vec \nabla \psi(x) - q\vec A(x)\, \psi(x) $$
If you act again from the left with $\vec p - q\vec A$, you just take the inner product.
Thinking about this in terms of gradients, divergences etc is possible but not too helpful, it's better to just think in components:
$$ (( \vec p - q\vec A )^2 \psi)(x) = \sum_{i=1}^3 \bigl( -\mathrm i\hbar\, \partial_i - q A_i(x) \bigr) \bigl( -\mathrm i\hbar\, \partial_i \psi(x) - q A_i(x)\, \psi(x) \bigr) $$
Best Answer
Note that the probability current in the presence of a EM field is given by
$$ \mathbf{j}=\frac{1}{2m}\left(\psi^{*}\mathbf{p}\psi-\psi\mathbf{p}\psi^{*}-2q\mathbf{A}\psi^{*}\psi\right) $$
As you note a local phase shift
$$ \psi^{\prime}=e^{iq\chi(\mathbf{r},t)/\hbar}\psi $$
leads to a gauge transformation of the vector potential
$$ \mathbf{A}^{\prime}=\mathbf{A}+\nabla\chi $$
Substituting these into the expression for the probability current gives
$$\begin{multline} \mathbf{j}^{\prime}=\frac{1}{2m}\left(e^{-iq\chi(\mathbf{r},t)/\hbar}\psi^{*}\mathbf{p}e^{iq\chi(\mathbf{r},t)/\hbar}\psi-e^{iq\chi(\mathbf{r},t)/\hbar}\psi\mathbf{p}e^{-iq\chi(\mathbf{r},t)/\hbar}\psi^{*}\right.\\\left.-2q\mathbf{A}e^{-iq\chi(\mathbf{r},t)/\hbar}\psi^{*}e^{iq\chi(\mathbf{r},t)/\hbar}\psi-2q\nabla\chi(\mathbf{r},t) e^{-iq\chi(\mathbf{r},t)/\hbar}\psi^{*}e^{iq\chi(\mathbf{r},t)/\hbar}\psi\right) \end{multline}$$
Operating with $\mathbf{p}\rightarrow-i\hbar\nabla$ one obtains
$$\begin{multline} \mathbf{j}^{\prime}=\frac{1}{2m}\left(\psi^{*}\left(-i\hbar\right)\nabla\psi+\psi^{*}\psi\frac{iq}{\hbar}(-i\hbar)\nabla\chi(\mathbf{r},t)-\psi\left(-i\hbar\right)\nabla\psi^{*}\right.\\\left.-\psi^{*}\psi\left(-\frac{iq}{\hbar}\right)(-i\hbar)\nabla\chi(\mathbf{r},t)-2q\mathbf{A}\psi^{*}\psi-2q\nabla\chi(\mathbf{r},t)\psi^{*}\psi\right) \end{multline}$$
Sorting it out one obtains
$$ \mathbf{j}^{\prime}=\frac{1}{2m}\left(\psi^{*}\mathbf{p}\psi-\psi\mathbf{p}\psi^{*}-2q\mathbf{A}\psi^{*}\psi\right)=\mathbf{j} $$ thus the probability current is gauge invariant.