The shape of groups of galaxies has been since now a disregarded topic in the literature, mainly due to the fact that individual ellipticity estimates of poor systems are strongly biased toward higher values. Such a bias turns out to depend on both the ellipticity itself and on the richness of the group. In this paper we have devised a modified version of the classical Lucy's iterative method, which is able to perform a statistical rectification of the observed ellipticity distribution of groups. Since this method needs an estimate of the average radial profile of the groups, we have faced also this problem. Our analysis has been applied to a sample of 54 groups belonging to the Tully (1987) catalog. The sample is selected in order to reduce the well-known, distance-dependent observational biases, which are present in all available group catalogs (Pisani et al. 1992). We found statistical evidence that galaxy groups share a common law for the radial distribution of members. This law can be quite satisfactorily represented by a simple Gaussian model, whereas a de Vaucouleurs profile (for instance) is definitely not consistent with the data. It is important to stress that the rectified distribution of apparent ellipticities looks very different from the original (unrectified) one. This implies that any statistical analysis which does not take into account the above mentioned bias, has to give meaningless results. The final aim of the paper is to produce the distribution of intrinsic ellipticities of groups, by deprojection of their apparent (rectified) ellipticity distribution. To do this we have to assume a biaxial (oblate or prolate) intrinsic shape of the groups are some model for their spatial orientation. Because of this reason, before the deprojection, we have searched for possible alignments of the group major axes along the supergalactic plane or toward the Local Supercluster center. We found a weak indication in favor of the latter hypothesis, but more sizable samples are needed to give statistical significance to this conclusion. At this stage, we have been forced to perform the deprojection by assuming a random orientation of the groups in the space. The apparent (rectified) ellipticity distribution which is obtained from our sample of galaxy groups, turned out to be consistent with both the oblate and the prolate intrinsic shapes. However, the oblate hypothesis should imply that a great fraction of groups has a very flattened, disklike structure. This circumstance is rather difficult to be explained within the available cosmological scenarios. Moreover, the convergence of the deprojection algorithm in the oblate case is hard-earned, and the final agreement with the data is less satisfactory than that obtained in the prolate case. All these facts led us to prefer the prolate solution with respect to the oblate one. Obviously a mixture of two different populations (oblate and prolate), as well as a triaxial solution, cannot be excluded.