Nonlinear normal modes (NNMs) are instrumental in spectral analysis, reduction, and identification of nonlinear vibrating systems, as they extend the notion of linear modes to nonlinear systems. Also called spectral submanifolds, they are (two-dimensional) invariant manifolds tangent to planar modes of the linearized dynamics. Moreover, NNMs have been studied under the lenses of Koopman operator theory, which aims at ``linearizing'' nonlinear dynamics through an infinite-dimensional lifting. In this context, it has been showed that NNMs can be computed as intersections of zero level sets of Koopman eigenfunctions. These results yielded a global parametrization of NNMs based on analyticity properties of the Koopman eigenfunctions and the so-called Koopman modes. In the present work, we develop a data-driven method to compute a global analytic parametrization of NNMs for nonlinear vibrating systems. To this aim, we combine the Koopman operator-based technique with the recent analytic extended dynamic mode decomposition (analytic EDMD) method which provides Taylor series of Koopman eigenfunctions and Koopman modes from data. While analytic EDMD is suited to dynamics with impulsive forcing (i.e. datasets containing several trajectories associated with different initial conditions), we also extend the proposed approach to dynamics with a general forcing term.