Open quantum systems play an important role in developing quantum sciences, and therefore the study of corresponding numerical method is of great significance. For the open quantum systems, the quasi-adiabatic propagator path integral invented in 1990s is one of the few numerically exact methods. However, its computational complexity scales exponentially with system size and correlation length, and therefore its application is limited in practical calculation. In recent years, the study and application of tensor network have made rapid progress. Representing the path integral by tensor network makes the computational complexity increase polynomially, thus greatly improving the computational efficiency. Such a new method is called time-evolving matrix product operator. At the very beginning, the reduced density matrix is represented as a matrix product state. Then the time evolution of the system can be achieved by iteratively applying matrix product operators to the matrix product state. The iterative process is amenable to the standard matrix product states compression algorithm, which keeps the computational cost on a polynomial scale. The time-evolving matrix product operator is an efficient, numerically exact and fully non-Markovian method, which has a broad application prospect in the study of quantum open systems. For instance, it is already used in the study of the thermalization, heat statistic, heat transfer and optimal control of the quantum open systems, and conversely it can be also used to investigate the effect of the system on the environment. In addition, the TEMPO method is naturally related to the process tensor, and can be used to calculate the correlation function of the system efficiently. In this article we review this method and its applications. We give a brief introduction of the path integral formalism of Caldeira-Leggett model. According to the path integral formalism, we demonstrate the usage of quasi-adiabatic propagator path integral method. we give the basic idea of matrix product states, and we show how to recast quasi-adiabatic propagator path integral method into time-evolving matrix product operators method by employing the concept of matrix product states and matrix product operators, and give a review of its applications. In addition, we use the calculation results of physical quantities, correlation functions and heat currents in the spin-boson model to illustrate the applications of the time-evolving matrix product operator method.