With the advantages of fast calculating speed and high precision, the physics-informed neural network method opens up a new approach for numerically solving nonlinear partial differential equations. Based on conserved quantities, we devise a two-stage PINN method which is tailored to the nature of equations by introducing features of physical systems into neural networks. Its remarkable advantage lies in that it can impose physical constraints from a global perspective. In stage one, the original PINN is applied. In stage two, we additionally introduce the measurement of conserved quantities into mean squared error loss to train neural networks. This two-stage PINN method is utilized to simulate abundant localized wave solutions of integrable equations. We mainly study the Sawada-Kotera equation as well as the coupled equations: the classical Boussinesq-Burgers equations and acquire the data-driven soliton molecule, M-shape double-peak soliton, plateau soliton, interaction solution, etc. Numerical results illustrate that abundant dynamic behaviors of these solutions can be well reproduced and the two-stage PINN method can remarkably improve prediction accuracy and enhance the ability of generalization compared to the original PINN method.