Indeed,
canonical quantization works just when it works.
It is in my view wrong and dangerous to think that this is the way to construct quantum theories even if it sometimes works: it produced astonishing results as the theoretical explanation of the hydrogen spectrum.
However, after all the world is quantum and classical physics is an approximation: the quantization procedures go along the wrong direction! There are in fact several no-go results against a naive validity of such procedures cumulatively known as Groenewold -Van Hove's theorem.
However, the question remains: why does that weird relation between Poisson brackets and commutators exist?
In fact, this relation motivates the naive quantization procedures.
In my view, the deepest answer relies upon the existence of some symmetry groups in common with classical and quantum theory.
These groups $G$ of transformations are Lie groups and they are therefore characterized by their so called Lie algebras $\mathfrak{g}$, which are vectors spaces equipped with a commutator structure $[a,b] \in \mathfrak{g}$ if $a,b\in \mathfrak{g}$. We can think of $a\in \mathfrak{g}$ as the generator of a one-parameter subgroup of $G$ usally denoted by $\mathbb{R} \ni t \mapsto \exp(ta) \in G$. If $a_1, \ldots, a_n \in \mathfrak{g}$ form a vector basis, it must hold $$[a_i,a_j] = \sum_k C^k_{ij}a_k\tag{1}\:,$$
for some real constants $C_k^{ij}$. These constants (almost) completely determine $G$.
For instance, if $G=SO(3)$ the group of 3D rotations, the one-parameters subgroups are rotations around fixed axes and it is always possible to choose $C_k^{ij}= \epsilon_{ijk}$ (the so-called Ricci symbol).
In classical physics, one represents the theory in the Hamiltonian formulation. States are points of a $2n$ smooth dimensional manifold $F$ called the space of phases, with prefereed classes of coordinates, said canonical, denoted by $q^1,\ldots, q^n, p_1,\ldots, p_n$.
If $G$ is a symmetry group of the system, then there is a faithful representation $G \ni g \mapsto \tau_g$ of it in terms of (canonical) transformations $\tau_g : F \to F$ which move the classical states according to the transformation $g$. The representation $G \ni g \mapsto \tau_g$ admits an infinitesimal description in terms of infinitesimal canonical transformations strictly analogous to the infinitesimal description of $G$ in terms of its Lie algebra $\mathfrak{g}$. In this case the corresponding of the Lie algebra is a linear space of smooth functions, $A \in C^\infty(F, \mathbb{R})$ representing classical observables, and the Poission bracket $\{A,B\} \in C^\infty(F, \mathbb{R})$.
An (actually central) isomorphism takes place between the Lie algebra $(\mathfrak{g}, [\:,\:])$ and the similar Lie algebra
$(C^\infty(F, \mathbb{R}), \{\:\:\})$ made of physical quantities where the commutator $\{\:\:\})$ is just the famous Poisson bracket.
If $a_k\in \mathfrak{g}$ corresponds to $A_k\in C^\infty(F, \mathbb{R})$ and (1) is valid for $G$, then
$$\{A_i,A_j\} = \sum_k C^k_{ij}A_k + c_{ij}1 \tag{2}$$
where the further constants $c_{ij}$, called central charges, depend on the representation.
$$a \mapsto A\tag{2'}$$ defines a (projective or central) isomorphism of Lie algebras.
When passing to the quantum description, if $G$ is still a symmetry group a similar mathematical structure exists.
Here, the space of (pure) states is a complex Hilbert space $H$ and the (pure) states are normalized vectors $\psi\in H$ up to phases.
If $G$ is a symmetry group there is a (projective/central) unitary representation $G \ni g \mapsto U_g$ in terms of unitary operators $U_g : H\to H$.
The one-parameter subgroups of $G$ are now represented by unitary groups of exponental form (I will systematically ignore a factor $1/\hbar$ in front of the exponent) $$\mathbb{R} \ni t \mapsto e^{-it \hat{A}}\:,$$
where $\hat{A}$ is a (uniquely determined) selfadjoint operator.
Again, if (1) is valid and $\hat{A}_k$ corresponds to $a_k\in \mathfrak{g}$, we have that
$$[-i\hat{A}_i,-i\hat{A}_j]= -i\sum_k C^k_{ij}\hat{A}_k -i c'_{ij}I \tag{3}$$
where $[\:,\:]$ is the commutator of operators.
In other words $$a \mapsto -i\hat{A} \tag{3'}$$ defines a (projective) isomorphism of Lie algebras.
I stress that the isomorphisms (2') and (3') independently exist and they are just due to the assumption that $G$ is a symmetry group of the system and the nature of the representation theory machinery.
Using these two isomprphisms, we can construct a third isomorphism (assuming $c_{ij}=c'_{ij}$) that interpolates between the classical and the quantum realm.
In this way, if $A \in C^\infty(F, \mathbb{R})$ corresponds to $\hat{A} : H \to H$ (actually one should restrict to a suitable dense domain), then
$$\{A,B\} \quad \mbox{corresponds to} \quad i[\hat{A},\hat{B}]\tag{4}$$
when comparing (2) and (3). (I again ignored a factor $\hbar$ since I have assumed $\hbar=1$ in the exponential expression of the one-parameter unitary groups.)
It is now clear that (4) is the reason of the correspondence principle of canonical quantization when the same symmetry group exists both in classical and in quantum physics.
In non relativistic physics, the relevant symmetry group is the Galileo group. This plays a crucial role both in classical and in non-relativistic quantum physics.
So we must have a (central) representation of its Lie algebra both in classical Hamiltonian and in Quantum physics.
Relying upon the above discussion, we conclude that
the isomorphism relating the isomorphic classical and quantum representations of the Galileo group -- the map associating classical quantities to corresponding operators preserving the commutation relations -- includes the so called canonical quantisation procedure
Let us illustrate this fact in details. The Lie algebra
$\mathfrak{g}$ includes a generator $p$ which, in classical Hamiltonian theory, describes the momentum (generator of the subgroups of translations) and another generator $k$ (generator of the subgroup of classical boost) corresponding to the position up to a constant corresponding to the mass of the system $m$.
Let us focus on the three levels.
Geometrically
$$[k,p]=0\:.$$
In the Hamiltonian formulation, a central charge show up
$$\{k,p\}= m 1$$
so that, defining $x:= k/m$, we have
$$\{x,p\}= 1\:.$$
In Quantum physics, in view of the discussion above, we should find for the corresponding generators/observables
$$[-i\hat{K},-i\hat{P}]= -im \hat{I}$$
hence, defining $\hat{X}:= \frac{1}{m}\hat{K}$,
$$[\hat{X},\hat{P}]= i \hat{I}$$
This correspondence, which preserves the commutation relation, can be next extended from the initial few observables describing the Lie algebra to a larger algebra of observables said the universal enveloping algebra. It is constructed out of the Lie algebra of the Galileo group. It includes for instance polynomials of observables.
Summing up: there are some fundamental symmetry groups in common with classical and quantum physics. These groups are the building blocks used to construct the theory, since they are deeply connected to basic notions as the concept of reference frame and basic physical principles as the relativity principle. The existence of these groups creates a link between classical and quantum physics. This link passes through the commutator structure of (projective) representations of the said group which is (projective) isomorphic to the Lie algebra of the symmetry group. Quantization procedures just reflect this fundamental relationship.
Next the two theories evolve along disjoint directions and, for instance, in quantum theory, further symmetry groups arise with no classical corresponding.
Best Answer
I do not know about deep questions. And people seem to give pretty deep answers here. My contribution is to show
$$ \lim_{\hbar \to \infty} \frac{1}{i\hbar}[ F(p,x) , G(p,x)] = \{F(p,x), G(p,x)\}_{P.B.} $$
where $ [ F, G] = F G - G F $ and
$$ \{ F(p,x), G(p,x) \} = \frac{\partial F}{\partial x} \frac{\partial G}{\partial p} - \frac{\partial G}{\partial x} \frac{\partial F}{\partial p}. $$
Preliminars.
With $[x, p] = i \hbar$, you can show the following two equalities:
$$ [x, f(p) ] = i \hbar \frac{\partial f}{\partial p} $$
and
$$ [p , g(x) ] = - i \hbar \frac{\partial g}{\partial x}. $$
I think this is almost mandatory for every QM course, so I will skip this derivation. In any case, the standard route is to consider the commutator of x with increasing powers of p; then use induction when developing $f(p)$ as a Taylor series.
A more illustrative example is the following:
$$ [x^{2} , f(p) ] = [x ,f(p) ] x + x [x, f(p)] \\ = i \hbar f'(p)\, x + i \hbar x \, f'(p) = 2 i\hbar x f'(p) - i\hbar [x , f'(p)]\\ = 2 i \hbar x f'(p) - (i \hbar)^{2} f''(p) $$
where I have introduced the pretty useful notation $ f'(p) = d f /dp $.
By now you can see the fun is in arbitrary powers of $x$. You should more or less be able to guess the result and prove it by induction.
Lemma.
$$ [x^{n} , f(p) ] = \sum_{j=1}^{n} (-)^{j+1} \binom{n}{k} \, (i \hbar)^{j} x^{n-j} \, f^{(j)}(p) $$
Proof: You do it. Use induction. It should be more or less straightforward. By the way, $\binom{n}{k}$ denotes the binomial coefficient.
Moment of truth.
The previous argument can be used to include an analytic function of $x$. Consider
$$ [ g(x) , f(p)] = \Biggl[ \, \sum_{k=0}^{\infty} \frac{1}{k!} g^{(k)} (0) x^{k}, \, f(p) \Biggr] = \sum_{k=0}^{\infty} \frac{1}{k!} g^{(k)} (0) \Biggl[ x^{k}, f(p) \Biggr] \\ = \sum_{k=0}^{\infty} \frac{1}{k!} g^{(k)} (0) \sum_{j=1}^{k} (-)^{j+1} C^{k}_{j} \, (i \hbar)^{j} x^{k-j} \, f^{(j)}(p) \\ = \sum_{j=1}^{\infty} (-)^{j+1} \, (i \hbar)^{j} g^{(j)}(x) \, f^{(j)}(p). $$
The trick in the fourth equality is to switch the sums (and then expand $C^{k}_{j}$... everything fits nicely).
It is interesting to notice that the double summations collapsed into one. This somehow makes sense by dimensional analysis, powers of x and p decrease together so that $\hbar$ appears.
The final part is the most subtle point. A general $f(x,p)$ is tricky because $x$ and $p$ do not commute. So you would have problems with "hermiticity" and ordering. I will choose every $p$ to be the left of every $x$. Once this is agreed, a general $F(p,x)$ can be written as
$$ F(p, x) = \sum_{n=0}^{\infty} \alpha_{n} (p) \,\, f_{n} (x). $$
Now, we can compute
$$ \Biggl[ F(p,x) , G(p,x) \Biggr] = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \Biggl[ \alpha_{n} (p) \,\,f_{n} (x), \, \beta_{m} (p) \, \, g_{m} (x) \Biggr] \\ = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \alpha_{n} (p) \biggl[ \, f_{n} (x) , \beta_{m} ( p) \biggr] g_{m} (x) + \beta_{m} (p) \, \biggl[ \, \alpha_{n} ( p) , g_{m} (x) \biggr] \, f_{n} (x) \\ = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} \, \alpha_{n} (p) \, \biggl( \sum_{j=1}^{\infty} (-)^{j+1} \, (i \hbar)^{j} f^{(j)}_{n} (x) \, \beta^{(j)}_{m} (p) \biggr) \, g_{m} (x) + \beta_{m} (p) \biggl( \sum_{j=1}^{\infty} (-)^{j} \, (i \hbar)^{j} g^{(j)}_{m}(x) \, \alpha^{(j)}_{n}(p) \biggr) \, f_{n} (x) $$
specially using
$$ \sum_{n=0}^{\infty} \alpha_{n} (p) \, f_{n}^{(j)} (x) = \frac{\partial^{j}}{\partial x^{j}} \biggl( \sum_{n=0}^{\infty} \alpha_{n} (p) \, f_{n} (x) \biggr) = \frac{\partial^{j}}{\partial x^{j}} F(p,x) $$
you see that you get the desired result (after changing the summations):
$$ \Biggl[ F(p,x), G(p,x) \Biggr] = \sum_{j=1}^{\infty} (-)^{j} \frac{(i \hbar)^{j}}{j!} \Biggl( \frac{\partial^{j} G}{\partial x^{j}} \frac{\partial^{j} F}{\partial p^{j}} - \frac{\partial^{j} F}{\partial x^{j}} \frac{\partial^{j} G}{\partial p^{j}} \Biggr) $$
because you see the only term that survives after dividing by (i \hbar) is the first one. This gives you the Poisson bracket. I didn't do any involved computations because they are long. It is more or less convincing.