Efficient prediction of mode shape variation under uncertainties is important for design and control. While Monte Carlo simulation (MCS) is straightforward, it is computationally expensive and not feasible for complex structures with high dimensionalities. To address this issue, in this study we develop a multi-fidelity data fusion approach with an enhanced Gaussian process (GP) architecture to evaluate mode shape variation. Since the process to acquire high-fidelity data from full-scale physical model usually is costly, we involve an order-reduced model to rapidly generate a relatively large amount of low-fidelity data. Combining these with a small amount of high-fidelity data altogether, we can establish a Gaussian process meta-model and use it for efficient model shape prediction. This enhanced meta-model allows one to capture the intrinsic correlation of model shape amplitudes at different locations by incorporating a multi-response strategy. Comprehensive case studies are performed for methodology validation.