We study two-scale parabolic partial differential equations whose coefficient is stochastic and depends linearly on a sequence of pairwise independent random variables which are uniformly distributed in a compact interval. We cast the problem into a deterministic two-scale parabolic problem which depends on a sequence of real parameters in a compact interval. Passing to the limit when the microscale tends to zero, using two-scale homogenization, we obtain the parametric two-scale homogenized problem. This problem contains the solution to the homogenized equation which describes the solution to the original two-scale parabolic problem macroscopically, and the corrector which encodes the microscopic information. The solution of this two-scale homogenized equation is represented as a generalized polynomial chaos (gpc) expansion according to a polynomial basis of the $L^2$ space of the parameters. We use a semidiscrete Galerkin approximation which projects the solution into a space of parametric functions which contain a finite number of pre-chosen gpc modes. Analyticity of the solution of the two-scale homogenized equation with respect to the parameters is established. Under mild assumptions, we show the summability of the coefficients of the solution’s gpc expansion. From this, an explicit error estimate for the semidiscrete Galerkin approximation in terms of the number of the chosen gpc modes is derived when these gpc modes are chosen as the $N$ best ones according to the norms of the gpc coefficients. Regularity of the gpc coefficients and summability of their norms in the regularity spaces are also established. Using the solution of the best $N$ term semidiscrete Galerkin approximation, we derive an approximation for the solution of the original two-scale problem, with an explicit convergence rate in terms of the microscopic scale and the number of gpc modes in the Galerkin approximation.