It is well known that a basis of the modular forms of weight $k$ over $\text{SL}_2(\mathbb{Z})$ is $\{E_{4}^iE_{6}^j: 4i + 6j = k\}$. Moreover, let
$$d=\dim\mathcal M_k(\text{SL}_2(\mathbb{Z})),$$
then we have a basis of the form
$$\{\varepsilon_1=1+O(q),\ \varepsilon_2=q+O(q^2),\ \cdots\ ,\ \varepsilon_d=q^{d-1}+O(q^d)\}.$$
For example, $\varepsilon_i=\Delta^{i-1}E_{k-12(i-1)}$, where $\Delta$ is the normalized cusp form of weight 12. It is easy to obtain that we have a basis of the form
$$\{\eta_1=1+O(q^d),\ \eta_2=q+O(q^d),\ \cdots\ ,\ \eta_d=q^{d-1}+O(q^d)\}.$$
My question is, can we write down each $\eta_i$ explicitly. Thanks in advance.
I think you're asking the completely practical question --- if we fix some $k$, can we explicitly write down a basis of modular forms over $\mathrm{SL}(2, \mathbb{Z})$ of weight $k$ using the Eisenstein series relations. The answer is yes, as long as $k$ is far too large.
We know the exact coefficients of the Eisenstein series. They are essentially $\sigma_4(n)$ and $\sigma_6(n)$ divisor functions. We also know the dimensions of the relevant spaces of Eisenstein series and modular forms through valence formulas.
One idea to write down modular forms is then to multiply out the various Eisenstein series and explicitly compute the first several ($k$ will do, less probably works too; see the Sturm bound for more) coefficients of each of the relevant forms $E_4^i E_6^j$. We can think of these as length $k$ vectors than span the space of the first $k$ coefficients of every modular form in $M_k(\mathrm{SL}(2, \mathbb{Z}))$.
The first $k$ coefficients of the elements that you call $\eta_i$ are what you get from reducing this matrix to row echelon form. If you track the row operations performed on the matrix, you can actually write each $\eta_i$ directly as linear combinations of the initial Eisenstein series products $E_4^i E_6^j$. Thus obtaining formulas for $\eta_i$ is an exercise in linear algebra.
I'll note that one can go further. We know how Hecke operators act on the coefficients. One can explicitly simultaneously diagonalize the space of modular forms and find a basis of Hecke eigenforms by studying the actions of the Hecke operators on the matrix formed from the first $k$ Fourier coefficients described above. More linear algebra!
In practice, this works fine for $\mathrm{SL}(2, \mathbb{Z})$, but other methods can be used to tabulate modular forms for more general congruence subgroups. See for example the LMFDB crowd (including me) paper on Computing classical modular forms for a description of other methods.