Biogas is a pervade renewable energy source in rural areas, which can be used to generate power, heating and gas to satisfy the multi-energy needs of users. However, the energy efficiency and thermal comfort of heat customers after the connection of biogas digester is a challenging modelling problem. This paper establishes a bi-level optimization model of integrated biogas energy system considering the thermal comfort of heat customers and the price fluctuation of natural gas. In detail, using the penalty factor to describe the thermal comfort and using the DBN (Dynamic Bayesian network) model to predict natural gas fluctuations. The upper optimization target is the largest profit of the energy agency, and the lower objective is the minimum heating cost considering the penalty factor. The uncertainty of weights and biases in DBN prediction is addressed using the Sparrow search algorithm, and the uncertainty in probability results is tackled using an improved Latin hypercube sampling method. The optimization results show the proposed bi-level model can greatly increase the profit of energy agency by 18.5%, and the DBN prediction model has 89.3% prediction accuracy. Furthermore, the influence of diverse biogas scales and operation modes of the system, fluctuated natural gas price effects and various penalty factors of the user economy and comfort during heating are also analyzed.