In FEM, to discrete the weak form, elements are needed to build the FE approximation basis space. FE method involves integration of the weak form. Due to discretization, the solution is an approximation close to the exact solution. I have a doubt that if FEM discretization error is mainly the error when us integrating with elements the bi-linear term in the weak form? If not the case, what is the other error source for the discretization error?
No, it is not equivalent. The main discretization error comes from restricting the space in which the solution is sought from a space of (weakly) differentiable functions to a subspace of piecewise polynomials -- this is Céa's lemma, which says that the discretization error is bounded by the best approximation error. Typically, these can then be integrated (and differentiated) exactly, so that there is no additional integration error.
If you do integrate the weak form inexactly (e.g., if you have spatially varying, non-polynomial, coefficients), the (first) Strang lemma can be used to estimate the additional error (which is usually of the same order as the approximation error, if the quadrature rule is exact for constant coefficients).
A possible other source of discretization error arises in the case that the space for the discrete solution is not a subspace of the space for the continuous solution (so-called non-conforming approaches such as, e.g., discontinuous Galerkin methods). The error from this "gap" can be estimated by the second Strang lemma.