We construct a framework for the study of fluctuations in the nonequilibrium relaxation of glassy systems with and without quenched disorder. We study two types of two-time local correlators with the aim of characterizing the heterogeneous evolution in these systems: in one case we average the local correlators over histories of the thermal noise, in the other case we simply coarse-grain the local correlators obtained for a given noise realization. We explain why the noise-averaged correlators describe the fingerprint of quenched disorder when it exists, while the coarse-grained correlators are linked to noise-induced mesoscopic fluctuations. We predict constraints on the distribution of the fluctuations of the coarse-grained quantities. In particular, we show that locally defined correlations and responses are connected by a generalized local out-of-equilibrium fluctuation-dissipation relation. We argue that large-size heterogeneities in the age of the system survive in the long-time limit. A symmetry of the underlying theory, namely invariance under reparametrizations of the time coordinates, underlies these results. We establish a connection between the probabilities of spatial distributions of local coarse-grained quantities and the theory of dynamic random manifolds. We define, and discuss the behavior of, a two-time dependent correlation length from the spatial decay of the fluctuations in the two-time local functions. We characterize the fluctuations in the system in terms of their fractal properties. For concreteness, we present numerical tests performed on disordered spin models in finite and infinite dimensions. Finally, we explain how these ideas can be applied to the analysis of the dynamics of other glassy systems that can be either spin models without disorder or atomic and molecular glassy systems.