A mathematical model of interaction of weakly non-spherical gas bubbles in liquid is proposed. It is a system оf the second order ODEs in the radii of bubbles, the position vectors of their centers and the amplitudes of their small deformations in the form of spherical harmonics. Compared to the available analogs, the proposed model equations are more accurate in terms of the ratio of the radii of bubbles to the distance between their centers. This allows one to study bubble-bubble interaction at smaller distances between the bubbles. The model equations are quite compact, which simplifies their theoretical analysis and construction of algorithms of their solution. Five problems (dynamics of a single bubble under harmonic excitation, bubble collapse near a plane rigid wall, interaction of three synchronized and nonsynchronized bubbles located in line, and interaction between sixteen bubbles with the centers located on a plane surface) are considered for validation. It is shown that in all cases the proposed model results well agree with the corresponding available experimental data and numerical solutions. For demonstration of the applicability, the proposed model is used to analyze radial oscillations, spatial movements and surface deformations of the bubbles in rectilinear, planar and spatial clusters under harmonic forcing. Initially 125, 121 and 125 bubbles of the clusters are located at the nodes of 1D, 2D and 3D uniform grids, respectively. It is found that the monotonically decaying radial oscillations of the 1D cluster bubbles during one period of forcing are quite similar to those of a single bubble. Unlike this, the attenuation of the radial oscillations of the central bubbles in the 2D and 3D clusters is essentially non-monotonic. The maximum pressures in the central bubbles of the 2D and 3D clusters are significantly higher than in the single and 1D cluster bubbles. In all clusters the maximum deformations of bubbles are in the form of the second spherical harmonics, except for the central bubble of the 3D cluster, the maximum deformations of which are in the form of the fourth harmonics.