Energy dependence in density functional theory

There are several things at play here. Kohn-Sham theory, as you know, involves creating an artificial system of non-interacting electron with the same density as that of the interacting system. Hohenberg-Kohn theorems, on the other hand, tell us that the energy of this system is a functional of the density. Your expression is correct, but in practice energy is computed by taking the sum of the eigenvalues $sum_j varepsilon_j,$ substracting the double counting of the electrostatic part (hartree and xc) and adding the electron-electron interaction energy. No need to compute all those integrals and the gradients of the eigenstates. In any case, yes, $E$ is a functional of the density because, after all, from HK theorems you know that the $varepsilon_j$ themselves are functions of the density. This does not mean that it is convenient to write explicitly $T_s$ in terms of the density! That would be very messy, and trying to do so will lead you to some approximations that might not be that good. To give a perhaps silly analogy: you know that $int_{-infty}^{x}e^{-y^2}dy$ is a function of x, but any attempt to write it explicitly in terms of $x$ (with elementary functions) will take you to approximations, never to the exact result. However, three comments are in order at this point:

1. Thomas-Fermi theory, which predates by 30 years DFT, did exactly that: they tried to write the kinetic energy in terms of the density, reasoning first for a model of uniform density, "LDA style", and then integrating.
2. An important improvement of the former, which accounts for spatial variation of the density, is the so-called Weizsacker's gradient correction. Actually, Kohn and Sham's original article mentions this, and they give a second-order expansion for the $T_s$!
3. Orbital-free DFT is a modern version of Thomas-Fermi model (with all the technicalities and formal aspects of DFT) which tries to do what you say: to write explicitly the dependency of $T_s$ in terms of $n.$ Generally, this implies losing some accuracy (the "real"/fully-correct explicit dependency of $T_s$ with $n$ could be something extremely complicated or not even expresable with elementary functions), but it might also speed up a lot the calculations, so that's why it is interesting, and depending on the system under study it might be very interesting to consider it as an alternative.

Although all of the above already (I hope) answers your question, I think there are some other interesting remarks that I can add: I would not say that the fact that E can written as a functional of the density follows from Kohn-Sham theory. Actually, the very expresion "functional of the density" is meaningless if one does not specify the functional space where these "densities" live. In principle, if you look at the most basic theory, the energy functional is defined only for $v-$representable densities; that is, densities given by (allow me to omit spin in this brief discussion): $$n(mathbf{r}) = Nint vert Psi(mathbf{r},mathbf{r}_2,cdots,mathbf{r}_N) vert^2 dmathbf{r}_2 cdots dmathbf{r}_N,$$ where $vert Psi rangle$ is the ground state of a hamiltonian $$hat{H}_v = hat{T}_{el}+hat{V}_{ee}+int v(mathbf{r})hat{n}(mathbf{r})dmathbf{r},$$ with $hat{T}_{el}$ the kinectic energy operator of the electrons, $hat{V}_{ee}$ the electron-electron (coloumbian) interaction operator and $v(mathbf{r})$ the function describing the external potential. Here, the density operator is defined as $hat{n}(mathbf{r})=sum_{j=1}^N delta(mathbf{r}-mathbf{r}_j),$ with $N$ the number of electrons in the system and the $mathbf{r}_j$ the position vector of each electronic degree of freedom, so that the last term could be written equivalently as $sum_{j=1}^N hat{v}(mathbf{r}_j).$

We say that such a density is $v-$representable and, actually, we know from the first Hohenberg-Kohn theorem that (up to a constant) they are in one-to-one correspondence with the set of admissible potentials. This set of $v-$representable densities is basically garbage. It does not even make sense to take functional derivatives in such a space, so we need to extend the domain of the functional. The way I see it, you can do this most comfortably in two stages: First, use the Levy-Lieb constrained search formulation of DFT to define an extension of the functional for the so called $Psi-$representable densities; that is, densities that can be written as: $$n(mathbf{r}) = Nint vert Psi(mathbf{r},mathbf{r}_2,cdots,mathbf{r}_N) vert^2 dmathbf{r}_2 cdots dmathbf{r}_N,$$ for antisymmetric $Psi,$ but without any reference to whether $Psi$ is a ground state of hamiltonians like $H_{v}$. The second stage involves proving that (not trivial!), actually, any "reasonable" density can be expressed in this way (and you can actually look for $Psi$ in the space of Slater determinants!). "Reasonable" here means some technical conditions, like being integrable and having integrable gradients. This takes you to the Sobolev spaces, which are Banach spaces where, finally, it makes sense to take functional derivatives. The fact that the Ground state density of $hat{H}_v$ is the minimizer of such extended functional follows then from the variational principle. I don't know many great detailed references about this (particularly the second stage; the constrained search formalism you will find in plenty of places) , but probably the best one to my knowledge, where one can read about DFT from the point of view of functional and convex analysis, is Eschring's book.