Dirac equation from a vierbein operator? Klein-Gordon equation can be derived straightforwardly by getting the mass-energy relation from special relativity in tensorial form, $$\eta^{\mu\nu}p_{\mu}p_{\nu} = m^2c^2$$
and promoting the variables to operators.
Dirac's motivation to find its equation was to find a first order equation instead of second order, so my idea was:
Can't we write the mass-energy relation using vierbeins? We get $$e^{\mu}p_{\mu} = \pm mc$$ and then we promote the variables to operators and get a first order equation for the wavefunction! If there isn't nothing wrong with my logic, what would be the vierbein operator? How is this connected to the Dirac equation?
 A: Short answer. A vierbein is a choice of a preferred inertial frame at a particular spacetime point. That is, it is a collection of four four-vectors. I insist on the two fours in the previous sentence. In fully expanded index notation, you’d write it as $e^\mu_\alpha$: sixteen components! So, on the surface, your formula doesn’t make sense because it has vectors on the left and scalars on the right.
Elaboration 1. What is a vierbein, really? In a Minkowski spacetime, there is a built-in way of identifying neighbourhoods of any two points (think about it literally: a point and its neighbours). It is consistent in the sense that going from $A$ to $B$ and then from $B$ to $C$ is the same thing as going directly from $A$ to $C$. In most other situations, there is no such natural prescription, and general relativity exploits this fact in a useful way.
How do you describe such a prescription, though, natural or not? A useful (mathematical) worldview tells you to keep things local if you can, so you consider tangent vectors—intuitively, “infinitely small” displacements—instead of finite neighbourhoods. Denote the (vector space of) tangent vectors at point $x \in M$ by $T_xM$; then our problem is to describe a consistent identification $T_xM \leftrightarrow T_yM$ for $x\not= y$. 
Choose a basis of tangent vectors at any point in spacetime, $\mathbf e_\alpha(x) \in T_xM, \alpha = 0,\ldots,3$, and you have an algorithm for moving vectors around: given a vector $\mathbf v \in T_xM$ at a point $x$, expand it in the local basis as $\mathbf v = v^\alpha\mathbf e_\alpha(x)$ and then use these coordinates to build a tangent vector at any other point: $v^\alpha\mathbf e_\alpha(y) \in T_yM$. Convince yourself that this prescription is consistent, and that any reasonable prescription can be described this way by choosing a basis at a single point and moving it to every other point. (What is reasonable?) To do physics, it is useful to limit ourselves to orthonormal (a.k.a. local inertial) frames with $\mathbf e_\alpha\cdot\mathbf e_\beta = \eta_{\alpha\beta}$.
Elaboration 2. Where are the sixteen components? Choose now a coordinate system on the spacetime, that is, an assignment to every1 point $x$ of a set of four numbers $x^\mu, \mu = 0,\ldots,3$. There is then a distinguished set of tangent vectors at every point of spacetime, the infinitesimal displacements along a particular coordinate conventionally denoted by $\partial_\mu\equiv\partial_\mu(x)$. These vectors form a basis, but they are not necessarily orthonormal: indeed, $\partial_r(r,\phi)$ has length one in polar coordinates, but $\partial_\phi(r,\phi)$ has length $r$. Any other vierbein we might have had can now be expanded in these:
$$ \mathbf e_\alpha(x) = e^\mu_\alpha(x)\partial_\mu(x). $$
Here you go, sixteen scalars $(e^\mu_\alpha)$. Their dependence on $x$ reflects the fact that the two prescriptions we now have can be different. And the indices $\alpha$ and $\mu$ are of completely different nature, too: they label different sets of basis vectors, so (e. g.) contracting them is as geometrically meaningful as $v^x\partial_r$. However, you can package the sixteen functions $e^\mu_\alpha$ into four sets of four both as $(e_\alpha)^\mu$, like we have done originally, and as $(e^\mu)_\alpha$, like in your formula (how can this be useful?).
Elaboration 3. What is a formula that does work? For Dirac’s (ultimately misguided) attempt at relativistic quantum mechanics, not only do you want a first-order equation instead of the energy-momentum relation: you also need the (second-order) relation to follow from it. A good way to ensure something second-order and linear follows from something first-order and linear is to demand that the latter squared equals the former2. Let’s try any inhomogeneous first-order operator, calling its coefficients $\gamma^\mu$ and $\delta$:
$$
  \left(p_\mu\cdot\gamma^\mu \pm m\cdot\delta\right)\psi = 0.
$$
The thing in the parentheses is an operator, so squaring it means applying it two times. You can’t multiply a scalar by a vector two times, so your original method can’t work. You can, however, pretend that these symbols multiply in some reasonable3 way and find out what that way must be so that $(\cdots)^2 = (p^2 - m^2)\cdot 1$. We need (although rarely write explicitly) the explicit $1$ in the last formula because we need to embed the scalar Klein-Gordon operator into the space those weird symbols live in. Work it out (you may assume $\delta = 1$ using the above definition of $1$) and you’ll get the correct result. Take care to avoid any prejudices about the $\gamma$s.

1 I pretend we only need one chart for the sake of simplicity: every argument here is local, so there’s little use concerning ourselves with global issues.
2 I don’t know if it’s the only one (modulo constant factors); edits welcome.
3 Associative, distributive and commuting with numbers and differential operators.

