We present a pipeline to infer the equation of state of neutron stars from observations based on deep neural networks. In particular, using the standard (deterministic), as well as Bayesian (probabilistic) deep networks, we explore how one can infer the interior speed of sound of the star given a set of mock observations of total stellar mass, stellar radius and tidal deformability. We discuss in detail the construction of our simulated data set of stellar observables starting from the solution of the post-Newtonian gravitational equations, as well as the relevant architectures for the deep networks, along with their performance and accuracy. We further explain how our pipeline is capable to detect a possible QCD phase transition in the stellar core. Our results show that deep networks offer a promising tool towards solving the inverse problem of neutron stars, and the accurate inference of their interior from future stellar observations.