The problem is that you're understanding the meaning of the coordinates in the hydrogen atom wave function incorrectly. The hydrogen atom begins with two coordinates, the positions of the electron and proton, $\mathbf{r}_e$ and $\mathbf{r}_p$, respectively. The very first step in solving the hydrogen atom is to produce the center of mass position, $\mathbf{R} \equiv \frac{m_p \mathbf{r}_p + m_e \mathbf{r}_e}{m_p + m_e}$. This coordinate is completely free, and so it, and its wave functions, are those of a free particle.
The next step is to produce a translation invariant coordinate, so that all changes from spatial translation are isolated in the center of mass coordinate. Since we have three vectors, there are a lot of different choices we could make, but the standard one is to use the difference of the electron and proton position:
$$\mathbf{r} \equiv \mathbf{r}_e - \mathbf{r}_p.$$
Under these two pairs of coordinates, the Hamiltonian takes the following forms:
\begin{align}
H & = \frac{\mathbf{p}_e^2}{2 m_e} + \frac{\mathbf{p}_p^2}{2 m_p} + V(|\mathbf{r}_e - \mathbf{r}_p|) \\
& = \frac{\mathbf{P}^2}{2 (m_e + m_p)} + \frac{\mathbf{p}^2}{2 \mu} + V(|\mathbf{r}|),
\end{align}
where $\mu$ is the reduced mass of the electron/proton system,
$$\mu \equiv \frac{1}{\frac{1}{m_e} + \frac{1}{m_p}},$$
$\mathbf{P}$ is the momentum conjugate to $\mathbf{R}$, and $\mathbf{p}$ is conjugate to $\mathbf{r}$.
You can see why we chose to define $\mathbf{r}$ the way we did - it's the coordinate of the potential, and this choice is the only one that makes the time independent Schrödinger equation (the eigenvalue equation for the Hamiltonian) a separable partial differential equation.
I stress again, $\mathbf{r}$ is completely invariant under translations, so it doesn't transform when you shift the origin of your coordinate system. It's also invariant under Galilean boosts, $\mathbf{r}_{e/p} \rightarrow \mathbf{r}_{e/p} + \mathbf{v} t$, which are just time parameterized translations. It is not invariant under Lorentz boosts, but we didn't expect it to be.
In the hydrogen ground and excited states, $r$ is the length of $\mathbf{r}$, so it is always $r^2 = r_e^2 + r_p^2 - 2r_e r_p \cos \alpha$, but in changing back to that basis you have to include the wave function for the position of the center of mass. Otherwise you have a mismatch between what is being described by the wave function (in this case, a probability amplitude for the relative positions) and the variables in it.
To be concrete, suppose a hydrogen atom is in its ground state and its center of mass is spread over a Gaussian wave packet, so:
\begin{align}
\Psi(\mathbf{R},\, \mathbf{r}) & = \psi_{\mathrm{com}}(\mathbf{R}) \, \psi_{\mathrm{internal}}(\mathbf{r}) \\
& = \frac{1}{\sqrt{\sigma_c^3 (2\pi)^{3/2}}}\mathrm{e}^{-R^2 / (4 \sigma_c^2)} \frac{2}{\sqrt{a_0^{3}}} \mathrm{e}^{-r / (2 a_0)}.
\end{align}
Changing variables to $\mathbf{r}_p$ and $\mathbf{r}_e$ requires multiplication by the square root of the determinant of the Jacobian matrix ($\Psi^2$ is part of a differential form):
\begin{align}
\Psi(\mathbf{R},\, \mathbf{r}) & = \Psi(\mathbf{r}_p,\, \mathbf{r}_e) \sqrt{\left|\frac{\partial (\mathbf{r}_p,\,\mathbf{r}_e)}{\partial (\mathbf{R},\, \mathbf{r}) }\right|} \\
& = \frac{2}{\sqrt{\sigma_c^3 (2\pi)^{3/2} a_0^3}} \exp\left(\frac{m_e^2 r_e^2 + m_p^2 r_p^2 + 2 m_e m_p r_e r_p \cos \alpha}{4 \sigma_c^2 (m_e+m_p)^2} - \frac{\sqrt{r_e^2 + r_p^2 - 2 r_e r_p \cos\alpha}}{2 a_0}\right).
\end{align}
Notice that since $\mathbf{r}_p \cdot \mathbf{r}_e = r_e r_p \cos\alpha$ the angle $\alpha$ depends on where in space both the proton and electron are. There is no special coordinate system where $\alpha = 0$, always.
