For a general vector bundle, I there is no "$E$-valued integration" as you put it. You are trying to add up elements in the fibres of $E$, but since the fibres over different points are not the same vector space you can't add their elements.
For the trivial bundle $M \times \mathbb{R}^k$ - and with a fixed choice of trivialisation! - you can carry out the integral component by component, But if you change the trivialisation you will get a different answer. Moreover, you can change the trivialisation in a way that varies over the manifold, so there is no hope that the integral will just change by a linear map of $\mathbb{R}^k$. You can see this behaviour even with ordinary functions. A function is a section of the trivialised rank 1 bundle. If you change the trivialisation, but insist on regarding ordinray $p$-forms as bundle valued, you multiply all your forms by a fixed nowhere vanishing function. This can change the integral over a $p$-cycle in a more-or-less arbitrary way.
What you can integrate $E$-valued forms against is $E^*$-valued ones. Given $a \in \Omega^p(E)$ and $b \in \Omega^q(E^*)$ their wedge product is an ordinary $(p+q)$-form which you can then integrate over a $(p+q)$-cycle. Now you have a version of Stokes theorem. If you have a connection $A$ in $E$ then you can check that
$$
d(a \wedge b) = d_A(a) \wedge b \pm a \wedge d_A(b).
$$
So Stokes theorem gives
$$
\int_{M}d_A(a) \wedge b \pm a \wedge d_A(b) = \int_{\partial M} a \wedge b.
$$
In the case when $b$ is a covariant constant section of $E$ and $M$ has dimension one more than the degree of $a$, we get
$$\int_M \langle d_A(a) , b \rangle=\int_{\partial M} \langle a, b \rangle$$
This is just the usual Stokes theorem, for the component of $a$ in the direction $b$. Since $b$ is covariant-constant, $\langle d_A(a), b \rangle = d \langle a, b \rangle$.
The Euler operator $d+d^*$ of a spin manifold is the spin Dirac operator
twisted by the $\mathbb Z/2$-graded spinor bundle $W=S^+\ominus S^-$, as explained in Heat Kernels and Dirac Operators by Berline, Getzler, Vergne [BGV, chapter 4.1]. The Chern character $\mathrm{ch}(W)$ equals $e(TM)\hat A(TM)^{-1}$, see [BGV] or Nicolaescu's notes. The cohomological index theorem therefore gives
$$\mathrm{ind}(d+d^*)^+=(\hat A(TM)\wedge\mathrm{ch}(W))[M]=e(TM)[M]$$
The de Rham representative of $e(TM)$ is given by the Pfaffian, and you get your second formula [BGV, Theorem 4.6].
The formulas above are pairings of "characteristic cohomology classes" with "fundamental homology classes". If you regard the Dirac operator as a fundamental $K$-homology class for the $K$-orientation of the tangent bundle given by the chosen spin structure, then the left hand side also constitutes a pairing of the $K$-theory class of the spinor bundle with the corresponding fundamental class.
To avoid the spin condition, you define the ``twist Chern character'' $\mathrm{tr}(\mathrm{exp}(-F^{W/S}))$ as explained in [BGV, chapter 4.1].
The formulas above still hold.
You can even get rid of orientability by regarding $e(TM)\in H^{\dim M}(M,o(TM))$ as a cohomology class twisted by the orientation bundle. Similarly, the Pfaffian is in fact a differential form twisted by the orientation bundle. Then all formulas above make sense without an orientation of $TM$.
On the other hand, the Euler number is the "trace" of id$_M$ in a sense explained nicely here (Def 2.2) - no matter if you count cells or dimensions of (co-) homology groups.
To relate it to the pairings above, you may use either Hodge theory, or the Gauss-Bonnet-Chern theorem, or deduce $e(TM)[M]=\chi(M)$ from the Poincar'e-Hopf theorem, e.g. using Morse theory. It seems that because $\chi(M)$ and $e(TM)[M]$ are different kinds of objects, some work is needed for this step.
Best Answer
When $E\to M$ is an oriented vector bundle of rank $2n$ over a compact manifold $M$, it has a well-defined de Rham Euler class $e(E)$ in $H^{2n}_{dR}(M)$, and a representative $2n$-form for $e(E)$ can be computed as follows:
Fix a positive definite inner product $\langle,\rangle$ on $E$. (Since any two such inner products are equivalent under automorphisms of $E$, it doesn't matter which one.) Let $\nabla$ be a $\langle,\rangle$-orthogonal connection on $E$, and let $K^\nabla$ denote the curvature of $\nabla$, regarded as a $2$-form with values in ${\frak{so}}(E)$. Let $e(\nabla) = \mathrm{Pf}(K^\nabla/2\pi)$, which is a well-defined $2n$-form on $M$. (N.B.: The definition of $\mathrm{Pf}$ requires both the inner product and the orientation of $E$.) Then $e(E)=\bigl[e(\nabla)\bigr]\in H^{2n}_{dR}(M)$. (That $\bigl[e(\nabla)\bigr]$ is independent of the choice of $\nabla$ is one of the first results proved in Chern-Weil theory.)
If $M$ is a compact, oriented $2n$-manifold, then the value of $e(E)$ on $[M]$, the fundamental class of $M$, can be computed as follows: Let $Y$ be a section of $E$ that has only a finite number of zeroes. (By Whitney transversality, the generic section of $E$ satisfies this condition.) Using the orientations of $E$ and $M$, one defines the index of an isolated zero $z$ of $Y$, which is an integer $\iota_Y(z)$. Then $$ e(E)\bigl([M]\bigr) = \int_M e(\nabla) = \sum_{z\in Z}\ \iota_Y(z). $$ By the Poincaré-Hopf Theorem, when $E = TM$ (as oriented bundles), this sum is equal to the Euler characteristic of $M$, which explains why $e(E)$ is called the Euler class of $E$.
To prove this result, one chooses an orthogonal connection $\nabla$ on $E$ that depends on $Y$ and whose Euler form $e(\nabla)$ is easy to evaluate explicitly. This can be done as follows:
Let $Z\subset M$ be the (finite) zero set of $Y$ and let $U\subset M$ be an open neighborhood of $Z$ that consists of disjoint, smoothly embedded $2n$-balls, one around each element of $Z$. Let $\phi$ be a smooth function on $M$ that is identically equal to $1$ on an open neighborhood of $Z$ and whose support $K\subset U$ is a disjoint union of closed, smoothly embedded balls, one around each element of $Z$.
Now, $E$ is trivial over $U$, so choose a positively oriented $\langle,\rangle$-orthonormal basis of sections $s_1,\ldots, s_{2n}$ of $E$ over $U$ and let $\nabla_1$ be the (flat) connection on $E_U$ for which these sections are parallel. Let $\bar Y = Y/|Y|$ be the normalized unit section defined on $M\setminus Z$, and define a second connection on $U$ by $$ \nabla_2 s = \nabla_1 s + (1{-}\phi)\bigl( \bar Y\otimes \langle s,\nabla_1 \bar Y\ \rangle - \langle s,\bar Y\ \rangle\ \nabla_1\bar Y\ \bigr). $$ It is easily verified that this formula does define a connection on $U$, that $\nabla_2$ is $\langle,\rangle$-orthogonal, and that, outside of $K$ (i.e., where $\phi\equiv0$), the vector field $\bar Y$ is $\nabla_2$-parallel. (Note that, because $\phi\equiv1$ near $Z$, where $Y$ is not defined, this formula does extend smoothly across $Z$, agreeing with $\nabla_1$ on a neighborhood of $Z$.)
Meanwhile, on $M\setminus K$, write $E$ as a $\langle,\rangle$-orthogonal direct sum $E = \mathbb{R}\cdot Y \oplus E'$ and choose a $\langle,\rangle$-compatible connection $\nabla_3$ on $E$ over $M\setminus K$ that preserves this splitting. Using a partition of unity subordinate to the open cover of $M$ defined by $U$ and $M\setminus K$, construct a $\langle,\rangle$-orthogonal connection $\nabla$ that agrees with $\nabla_2$ on $K$ and with $\nabla_3$ on $M\setminus U$, and for which the line bundle $\mathbb{R}\cdot Y\subset E$ is parallel on $M\setminus K$.
Now, on $M\setminus K$, the curvature of $\nabla$ takes values in ${\frak{so}}(E')\subset{\frak{so}}(E)$, so $e(\nabla) = \textrm{Pf}\bigl(K^{\nabla}/(2\pi)\bigr)$ vanishes identically outside of $K$. It suffices now to evaluate the integral of $e(\nabla)$ over a single component $B$ of $K$, which may be assumed to be the unit ball in $\mathbb{R}^{2n}$, so restrict attention to this case. Write $\bar Y = s_1 u_1 +\cdots + s_{2n}\ u_{2n}$, and note that, on $B$, one has, by definition, $$ \nabla s_i = \nabla_2 s_i = \sum_{j=1}^{2n} s_j\otimes (1{-}\phi)(u_j\ du_i - u_i\ du_j), $$ i.e., the connection $1$-forms of $\nabla$ in this basis are $\omega_{ij} = (1{-}\phi)(u_i\ du_j - u_j\ du_i)$.
Using the identity ${u_1}^2 +\cdots + {u_{2n}}^2 = 1$, the curvature forms are easily computed to satisfy $$ \Omega_{ij} = d\omega_{ij} + \omega_{ik}\wedge\omega_{kj} = (1{-}\phi^2)\ du_i\wedge du_j - d\phi\wedge(u_i\ du_j - u_j\ du_i). $$ At this point, you have to know the definition of the Pfaffian. (I'll wait while you look it up.) Using the fact that the result has to be invariant under $SO(2n)$-rotations, it is easy to show that, on $B$, one has $$ e(\nabla) = \textrm{Pf}\left(\frac{\Omega}{2\pi}\right) = c_n (1{-}\phi^2)^{n-1} d\phi\wedge u^*(\Upsilon) $$ for some universal constant $c_n$ and where $\Upsilon$ is the $SO(2n)$-invariant $2n$-form of unit volume on $S^{2n-1}\subset\mathbb{R}^{2n}$ and $u: B\setminus z\to S^{2n-1}$ is $u = (u_1,\ldots, u_{2n})$. (I'll let you evaluate the constant $c_n$. This can be done in a number of ways, but it's essentially a combinatorial exercise.) In particular, $$ e(\nabla) = d\bigl(P_n(1{-}\phi)\ u^*(\Upsilon)\bigr) $$ where $P_n(t)$ is the polynomial in $t$ (of degree $2n{-}1)$) that vanishes at $t=0$ and satisfies $P'_n(1{-}t) = -c_n(1{-}t^2)^{n-1}$. By Stokes' Theorem, one has $$ \int_B e(\nabla) = P_n(1) \int_{\partial B}u^*(\Upsilon) = P_n(1)\ \textrm{deg}(u:\partial B\to S^{2n-1}) = P_n(1)\ \iota_Y(z) $$ Thus, it follows that there is a universal constant $C_n = P_n(1)$ such that $$ e(E)\bigl([M]\bigr) = \int_M e(\nabla) = C_n \ \sum_{z\in Z}\ \iota_Y(z). $$ Now, $C_1 = 1$, as one can show by computing with the constant curvature metric on the $2$-sphere and noting that, for $Y$ the gradient of the height function on $S^2$, the sum of the indices is $2$. Finally, by properties of the Pfaffian and the index, one sees that, if such a formula as above is to hold, then one must have $C_{n+m}=C_nC_m$. Thus, $C_n=1$ for all $n$, and the formula is proved.
(Of course, one can avoid these final tricks for evaluating the constants by carefully doing the combinatorial exercise, but I'll leave that to the curious.)