Renewable energy sources will anticipate significantly in the future energy system paradigm due to their low cost of operation and low pollution. Considering the renewable generation (e.g., wind) intermittency, flexible gas-fired power plants will continue to play their essential role as the main linkage of natural gas and electricity networks, and hence coordinated operation of these networks is beneficial. Furthermore, uncertainty is always found in gas demand prediction, electricity demand prediction, and output power of wind generation. Therefore, in this paper, a two-stage stochastic model for operation of natural gas and electricity networks is implemented. In order to model uncertainty in these networks, Monte Carlo simulation is applied to generate scenarios representing the uncertain parameters. Afterwards, a scenario reduction algorithm based on distances between the scenarios is applied. Stochastic and deterministic models for natural gas and electricity networks are optimized and compared considering integrated and iterative operation strategies. Furthermore, the value of flexibility options (i.e., electricity storage systems) in dealing with uncertainty is quantified. A case study is presented based on a high pressure 15-node gas system and the IEEE 24-bus reliability test system to validate the applicability of the proposed approach. The results demonstrate that applying the stochastic model of gas and electricity networks as well as considering integrated operation strategy in the presence of flexibility provides different benefits (e.g., 14% cost savings) and enhances the system reliability in the case of contingency.