In a previous article [1] a method has been introduced to derive the all order Bose-Einstein distribution of the non interacting Bosons as the solution of the Wigner equation. The process was a perturbative one where the Bose-Einstein distribution was taken as the unperturbed solution. In this article it is shown that the same formalism is also applicable in the case of interacting Bosons. The formalism has been applied to calculate the quantum second virial coefficient of the Bosons interacting pairwise via Lenard-Jones potential and compared with the previous result.