Many systems, including biological tissues and foams, are made of highly packed units having high deformability but low compressibility. At two dimensions, these systems offer natural tesselations of plane with fixed density, in which transitions from ordered to disordered patterns are often observed, in both directions. Using a modified Cellular Potts Model algorithm that allows rapid thermalization of extensive systems, we numerically explore the order-disorder transition of monodisperse, two-dimensional cellular systems driven by thermal agitation. We show that the transition follows most of the predictions of Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory developed for melting of 2D solids, extending the validity of this theory to systems with many-body interactions. In particular, we show the existence of an intermediate hexatic phase, which preserves the orientational order of the regular hexagonal tiling, but looses its positional order. In addition to shedding light on the structural changes observed in experimental systems, our study shows that soft cellular systems offer macroscopic systems in which KTHNY melting scenario can be explored, in the continuation of Bragg's experiments on bubble rafts.