The first stages of the path instability phenomenon affecting the buoyancy-driven motion of gas bubbles rising in weakly or moderately viscous liquids are examined using a recently developed numerical approach designed to assess the global linear stability of incompressible flows involving freely evolving interfaces. Predictions for the critical bubble size and frequency of the most unstable mode are found to agree well with reference data obtained in ultrapure water and in several silicone oils. By varying the bubble size, stability diagrams are built for several specific fluids, revealing three distinct regimes with different bifurcation sequences. The spatial structure of the unstable modes is analysed, together with the variations of the bubble shape, position and orientation. For this purpose, displacements of the bubble surface are split into rigid-body components and volume-preserving deformations, allowing us to determine how the relative magnitude of the latter varies with the fluid properties and bubble size. Predictions obtained with freely deformable bubbles are compared with those found by maintaining the bubble shape determined in the base state frozen during the stability analysis. This comparison reveals that deformations leave the phenomenology of the first bifurcations unchanged in low-viscosity fluids, especially water. Hence, in such fluids, bubbles behave essentially as freely moving rigid bodies submitted to constant-force and zero-torque constraints, at the surface of which the fluid obeys a shear-free condition. In contrast, deformations change the nature of the primary bifurcation in oils slightly more viscous than water, whereas, somewhat surprisingly, they leave the near-threshold phenomenology unchanged in more viscous oils.