If the Shimura variety is attached to the Shimura data $h:\mathbb S \to G_{/\mathbb R}$,
and if as usual $K$ denotes the centralizer in $G_{/\mathbb R}$ of $h$,
then the automorphic vector bundles are $G(\mathbb R)$-equivariant bundles on $X:= G(\mathbb R)/K(\mathbb R)$ attached to representations of the algebraic
group $K$.
Specifically, if $V$ is the representation (over $\mathbb C$, say), then the associated
vector bundle $\mathcal V$ is given by $(G(\mathbb R)\times V)/K(\mathbb R).$
Note that if $G = GL_2$ and we are in the modular curve case, then $K$ is equal to
$\mathbb C^{\times}$ (thought of as an algebraic group over $\mathbb R$), its complexification is $\mathbb G_m\times \mathbb G_m$, and so there are two integral
parameters describing its irreps (which are just one-dimensional in this case). One of them is the usual $k$ of $\omega^k$; the other can be chosen so as to correspond to twisting by a power of the determinant, which does not change the underlying
line bundle, but changes the equivariant structure (which we don't think about so explicitly in the classical language).
Now $K$ is the Levi in a parabolic $P$ (defined over $\mathbb R$), and there is an open embedding (of complex analytic spaces)
$X := G(\mathbb R)/K(\mathbb R) \hookrightarrow G(\mathbb C)/P(\mathbb C) =: D$.
(If you look at Deligne's Corvalis article and unravel things, you'll see that this is
how the complex structure on $X$ is defined: in the Hodge-theoretic formulation given there, $P$ is the parabolic preserving the Hodge filtration on the Hodge-structure corresponding to the base point of $X$.)
The automorphic vector bundles are naturally defined on the flag variety $D$
(as $\mathcal V := (G(\mathbb C)\times V)/P(\mathbb C)$), and then restricted to $X$.
(So in the modular curve case, $D$ is the projective line, and $\omega^k$ comes from
$\mathcal O(k)$.)
Now the automorphic bundles, being $G$-equivariant, descend to bundles on each
Shimura variety quotient $Sh(X,K_f)$ of $X \times G(\mathbb A_f)$, and have canonical models defined over the reflex field. In the PEL case, one will be able to construct these canonical models using data from the abelian varieties (think about how the abelian varieties give rise to Hodge structures with Mumford--Tate group equal to $G$ in the first place). In general, this result is due to Milne (see his answer).
Added: A colleague pointed out to me, regarding Kevin's initial question,
that it can be understood much more generally. Namely, all automorphic
vector bundles are locally free coherent sheaves, equivariant under the
action of the finite adeles, but not all coherent sheaves necessarily
have these properties! Practically nothing is known about the $K$-theory of
Shimura varieties, although the question is obviously of fundamental
interest.
The very deep current work on cycles on Shimura varieties, from various
points of view, must be the beginning of a substantial theory whose ultimate
shape no one is in a position to imagine. So Kevin's question should be
seen
as a program for future collaboration between number theorists and algebraic
geometers (at least).
Also, let me remark that I'm told that the term "automorphic vector bundle" was invented in
a conversation between Michael Harris and Jim Milne in the second half of the 80s.
Let $S = S_{\mathbf{Q}} = M_{13}(\Gamma_1(N),\mathbf{Q})$, and $S_{\mathbf{C}} = S \otimes \mathbf{C}$ denote the corresponding space of modular forms over $\mathbf{C}$.
Let $V \subset S \times S$ be the subspace cut out by pairs of forms $(A,B)$ satisfying the following equation:
$$A \cdot E_{12} = B \cdot \Delta $$
As equations in the Fourier coefficients of $A$ and $B$ these are linear equations with coefficients in $\mathbf{Q}$. Since, by the $q$-expansion principle, a modular form can be recovered from some finite number of Fourier coefficients, $V$ is determined by the null space of some finite matrix with coefficients in $\mathbf{Q}$. Since a linear system over $\mathbf{Q}$ has the same rank over $\mathbf{C}$, it follows that $V_{\mathbf{C}} = V \otimes \mathbf{C}$, where $V_{\mathbf{C}}$ is the set of solutions in $S_{\mathbf{C}} \times S_{\mathbf{C}}$ of the same equations.
On the other hand, there is an isomorphism:
$$V_{\mathbf{C}} \rightarrow M_{1}(\Gamma_1(N),\mathbf{C})$$
given by
$$(A,B) \mapsto \frac{A}{\Delta} = \frac{B}{E_{12}}$$
The point is that $E_{12}$ and $\Delta$ do not have any common zeros, so the image of this map clearly consists of holomorphic forms. Hence the map is well defined. However, if $F$ has weight one, then $(A,B) = (F \cdot \Delta, F \cdot E_{12})$ maps to $F$, so the map is surjective. It is clearly injective, so it is an isomorphism.
It follows that the image of $V$ under this map gives a rational basis for $M_1(\Gamma_1(N))$. Since $V$ is then preserved by Hecke operators (as is obvious on $q$-expansions), the result follows.
Best Answer
Here is one construction:
We have the exact sequence $$0 \to H^0(\omega^{\otimes k}(-\text{cusps})) \to H^0(\omega^{\otimes k}) \to H^0(\text{cusps}, \omega^{\otimes k}_{| \text{cusps}}).$$ (Here I am using $\omega$ for what you called $E$; this is the traditional notation for modular forms people.) It is easy to define a Hecke action on the third $H^0$ so that this exact sequence is Hecke equivariant.
The right hand map is surjective if $k > 2$, and its image has codimension one when $k = 2$. In any event, write $\mathcal I$ to denote the image, so that $$0 \to H^0(\omega^{\otimes k}(-\text{cusps})) \to H^0(\omega^{\otimes k}) \to \mathcal I \to 0$$ is short exact. One then shows that this short exact sequence has a unique Hecke equivariant splitting; i.e. there is a uniquely determined Hecke equivariant subspace $\mathcal E \subset H^0(\omega^{\otimes k})$ such that $\mathcal E$ projects isomorphically onto $\mathcal I$. This space $\mathcal E$ is the space of weight $k$ Eisenstein series (for whatever level we are working at).