Fractional Gaussian noise models the time series with long-range dependence; when the Hurst index H > 1/2, it has positive correlation reflecting a persistent autocorrelation structure. This paper studies the numerical method for solving stochastic fractional diffusion equation driven by fractional Gaussian noise. Using the operator theoretical approach, we present the regularity estimate of the mild solution and the fully discrete scheme with finite element approximation in space and backward Euler convolution quadrature in time. The O(τ H−ρα ) convergence rate in time and O(h min(2,2−2ρ, H α ) ) in space are obtained, showing the relationship between the regularity of noise and convergence rates, where ρ is a parameter to measure the regularity of noise and α ∈ (0, 1). Finally, numerical experiments are performed to support the theoretical results.