A direct numerical simulation study has been carried out to better understand the dynamics of deformable bubbles in bubble-induced turbulence. A mathematical model for bubble collision, bounce and coalescence is developed based on the Volume-of-Fluid (VoF) method. A surface-tension-holding film thickness δc is introduced to determine whether the bubbles will coalesce or bounce after collision. Bubbly flows with the gas volume fractions αg=5.6% (low), 11.2% (medium) and 22.1% (high) are considered in the study. The DNS results show that the dissipation of turbulent kinetic energy in the cases under consideration takes place when the wavenumber κ is smaller than twice the wavenumber for the bubble diameter 2κd. The bubble-induced turbulence is anisotropic at the length scales close to the bubble diameter due to the elongated turbulent structures in the bubble-rising direction. The bubbles collide intensively with each other when they are in a cluster surrounded by the liquid with low pressure. Parallel clustering of bubbles is found when the gas volume fraction has a low or medium value. Three-dimensional clustering is found when the bubbles are densely populated. Two different mechanisms of bubble collision have been identified from the DNS results, termed parallel collision and vertical collision. Parallel collision is often observed when bubbles are sparsely populated. In a parallel collision, the relative velocity of the bubbles slows down as two bubbles approach each other due to the jet flow between them, while the relative velocity increases sharply after the bounce. In a vertical collision, by contrast, the relative velocity of bubbles accelerates as two bubbles approach each other, while it slows down during the collision. Vertical collisions occur when the bubbles are more densely populated. The numerical results also show the significant effects of δc on bubble coalescence.