Given a hyperspectral image, unmixing tries to estimate the spectral responses of the latent constituent materials and their corresponding fractions. Recently, Nonnegative Matrix Factorization (NMF) has been widely applied to solve the hyper-spectral unmixing problem because of its plausible physical interpretation. In this paper, we propose a novel method, Multilayer Manifold and Sparsity constrained Nonnegative Matrix Factorization (MMSNMF), for hyperspectral unmixing. In this approach, Multilayer NMF decomposes a hyperspectral image iteratively at several layers. In order to consider both the manifold structure of hyperspectral image and the sparsity of abundance matrix, we impose a graph regularization term and a sparsity regularization term on both the spectral signature matrix and the abundance matrix. Experimental results on both synthetic and real data validate the effectiveness of the proposed method in hyperspectral unmixing.