wsc's answer (i.e., Onsager's computation of the free energy) provides one alternative road to a proof of a phase transition in the Ising model. It implies the existence of a phase transition in dimension 2 (for the nearest-neighbor model). Combined with correlation inequalities, this implies existence of a phase transition in any dimension d≥2, and interactions of any range (provided that nearest-neighbors also interact).
An alternative approach is via reflection positivity, see, for example, these lecture notes. In particular Theorem 3.5 therein implies the existence of a phase transition in all $O(n)$ models in dimensions $3$ and more. This approach (via the so-called infrared bound) cannot yield the result in dimension $2$, however (but notice that the $O(n)$ model do not display spontaneous magnetization in dimension 2 when $n\geq 2$).
You can also use slight variants of the Peierls approach, using percolation arguments. For example, you can use a comparison (relying on the FKG inequality available in the random-cluster representation of the Ising model) with the percolation model, in order to show that the existence of a phase transition in the percolation model implies the existence of a phase transition in the Ising model. However, to prove the former, you use again a Peierls-type argument...