The study of uncertainty propagation poses a great challenge to design high fidelity numerical methods. Based on the stochastic Galerkin formulation, this paper addresses the idea and implementation of the first flux reconstruction scheme for hyperbolic conservation laws with random inputs. High-order numerical approximation is adopted simultaneously in physical and random space, i.e., the modal representation of solutions is based on an orthogonal polynomial basis and the nodal representation is based on solution collocation points. Therefore, the numerical behaviors of the scheme in the (physical-random) phase space can be designed and understood uniformly. A family of filters is developed in multi-dimensional cases to mitigate the Gibbs phenomenon arising from discontinuities in both physical and random space. The filter function is switched on and off by the dynamic detection of discontinuous solutions, and a slope limiter is employed to preserve the positivity of physically realizable solutions. As a result, the proposed method is able to capture the stochastic flow evolution where resolved and unresolved regions coexist. Numerical experiments including a wave propagation, a Burgers’ shock, a one-dimensional Riemann problem, and a two-dimensional shock-vortex interaction problem are presented to validate the current scheme. The order of convergence of the high-order scheme is identified. The capability of the scheme for simulating smooth and discontinuous stochastic flow dynamics is demonstrated. The open-source codes to reproduce the numerical results are available under the MIT license (Xiao et al. in FRSG: stochastic Galerkin method with flux reconstruction. https://github.com/CSMMLab/FRSG, (2021). https://doi.org/10.5281/zenodo.5588317).