Mode shape information plays the essential role in deciding the spatial pattern of vibratory response of a structure. The uncertainty quantification of mode shape, i.e., predicting mode shape variation when the structure is subjected to uncertainty, can provide guidance for robust design and control. Nevertheless, computational efficiency is a challenging issue. Direct Monte Carlo simulation is unlikely to be feasible especially for a complex structure with a large number of degrees-of-freedom. In this research, we develop a new probabilistic framework built upon the Gaussian process meta-modeling architecture to analyze mode shape variation. To expedite the generation of input data set for meta-model establishment, a multi-level strategy is adopted which can blend a large amount of low-fidelity data acquired from order-reduced analysis with a small amount of high-fidelity data produced by high-dimensional full finite element analysis. To take advantage of the intrinsic relation of spatial distribution of mode shape, a multi-response strategy is incorporated to predict mode shape variation at different locations simultaneously. These yield a multi-level, multi-response Gaussian process that can efficiently and accurately quantify the effect of structural uncertainty to mode shape variation. Comprehensive case studies are carried out for demonstration and validation.