Proof that $\pi =\lim_{n\to\infty}\frac{2^{4n}n!^4}{n(2n)!^2}$

You are right that not all operations with asymptotically equivalent expressions are correct. For example, $f\sim g$ doesn't necessarily imply $e^f\sim e^g$. However, $f\sim g$ implies $f^k\sim g^k$ where $k$ is an integer. This follows directly from properties of limits: $$\lim_{x\to a}\frac{f(x)^k}{g(x)^k}=\left(\lim_{x\to a}\frac{f(x)}{g(x)}\right)^k=1^k=1 $$ (This can be generalized to $k\in\mathbb{R}$ provided $f(x)^k$ is defined). Since the only thing we're doing here is squaring, i.e. $k=2$, your manipulations are justified.

You can avoid the $\sim$ by writing the Stirling formula in the form $$n!=\sqrt{2\pi n}\left({n\over e}\right)^n\>g(n)\quad(n\geq1),\qquad \lim_{n\to\infty}g(n)=1\ ,$$ and similarly in the other equations, maybe with new $g$s.

The asymptotic equality $e^{f(n)} \sim e^{g(n)}$ means $$ \frac{{e^{f(n)} }}{{e^{g(n)} }}\to 1 \Leftrightarrow e^{f(n) - g(n)} \to 1 \Leftrightarrow f(n) - g(n) \to 0. $$ But $f(n) \sim g(n)$ does not necessarily imply that the difference of the two sides tend to $0$ (or even that it is bounded). Thus $f(n) \sim g(n) \not\Rightarrow e^{f(n)} \sim e^{g(n)}$. If, for example, $g(n)$ is evetually non-zero, then the reverse implication is true.

In your example $$ \log n! - (n\log n - n) = \frac{1}{2}\log (2\pi n) + \mathcal{O}\!\left( {\frac{1}{n}} \right) $$ which is not bounded.