We investigate the nonequilibrium stationary states of systems consisting of chemical reactions among molecules of several chemical species. To this end we introduce and develop a stochastic formulation of nonequilibrium thermodynamics of chemical reaction systems based on a master equation defined on the space of microscopic chemical states, and on appropriate definitions of entropy and entropy production, The system is in contact with a heat reservoir, and is placed out of equilibrium by the contact with particle reservoirs. In our approach, the fluxes of various types, such as the heat and particle fluxes, play a fundamental role in characterizing the nonequilibrium chemical state. We show that the rate of entropy production in the stationary nonequilibrium state is a bilinear form in the affinities and the fluxes of reaction, which are expressed in terms of rate constants and transition rates, respectively. We also show how the description in terms of microscopic states can be reduced to a description in terms of the numbers of particles of each species, from which follows the chemical master equation. As an example, we calculate the rate of entropy production of the first and second Schlogl reaction models.