We study the effect of the non-Gaussianity induced by gravitational evolution upon the statistical properties of absorption in quasar (quasi-stellar object) spectra. Using the generic hierarchical ansatz and the lognormal approximation, we derive the analytical expressions for the one-point probability distribution function (PDF) as well as for the joint two-point PDF of transmitted fluxes in two neighbouring quasi-stellar objects. These flux PDFs are constructed in three dimensions as well as in projection (i.e. in two dimensions). The PDFs are constructed by relating the lower-order moments (i.e. cumulants and cumulant correlators) of the fluxes to the three-dimensional neutral hydrogen distribution, which is, in turn, expressed as a function of the underlying dark matter distribution. Next, the lower-order moments are modelled using a generating function formalism in the context of a minimal tree-model for the higher-order correlation hierarchy. These different approximations give nearly identical results for the range of redshifts probed, and we also find very good agreement between our predictions and the outputs of hydrodynamical simulations. The formalism developed here for the joint statistics of flux-decrements concerning two lines of sight can be extended to multiple lines of sight, which could be particularly important for the three-dimensional reconstruction of the cosmic web from the spectra of quasi-stellar objects (e.g. in the Baryon Oscillation Spectroscopic Survey). These statistics probe the underlying projected neutral hydrogen field and are thus linked to hotspots of absorption. The results for the PDF and the bias presented here use the same functional forms of scaling functions that have previously been employed for the modelling of other cosmological observations, such as the Sunyaev-Zel'dovich effect.