With breakthroughs in scientific computing, computer simulations are quickly replacing physical experiments in modern scientific problems. Such simulations, however, are very time-consuming, and predictive models are used to emulate the expensive computer code. In many problems, these simulations are often performed in multiple stages, with the accuracy (or fidelities) at each stage controlled by a tuning parameter. This provides a flexible multi-stage, multi-fidelity framework for efficiently simulating lower-fidelity training data. We propose a new Multi-stage Multi-Fidelity Gaussian process (M$^2$GP) model, which leverages this multi-stage, multi-fidelity simulation data to efficiently train an emulator for the high-fidelity computer code. The M$^2$GP uses a generalized Brownian sheet model to embed prior information on error convergence from numerical analysis. We propose an efficient algorithmic framework for M$^2$GP prediction, and provide novel design methods for allocating simulation runs given a computational budget. The improvement of M$^2$GP over existing methods is then demonstrated in numerical experiments and an application to emulation of heavy-ion collisions.