Clustering is widely observed in many turbulent flows, where it results from the inability of finite inertia particles to comply with the different time scales, which characterize a turbulent field. Depending on their inertia, particles are found to be instantaneously organized in clusters, whose size depends on the Kolmogorov-Stokes number and which presumably form as a consequence of particle ejection from persistent vortical structures. In reacting flows, the abrupt acceleration of the fluid across the thin flame front due to combustion adds new and unexpected features. The particles follow such acceleration with a certain time lag which, coupled with the flame front fluctuations, gives rise to an entirely different mechanism of cluster formation. As suggested in previous studies, a possible indicator of this preferential localization is the so-called clustering index, quantifying the departure of the actual particle arrangement from the Poissonian distribution. Most of the clustering is found in the flame brush region, where it cannot be explained by the standard arguments used in cold flows. Actually, the effect is significant also for very light particles, where the simple model we propose, based on the Bray-Moss-Libby formalism, is able to account for most of the deviation from the Poissonian. When the particle inertia increases, the effect becomes larger and it is found to persist well within the region of the burned gases. The observed clustering is confirmed by a more precise analysis in terms of a generalization of the radial distribution function to inhomogeneous, anisotropic flows. The results taken from a direct numerical simulation with single step kinetics favorably compare with experiments on a premixed Bunsen turbulent flame. The present findings are expected to be of some relevance for the plenty of applications dealing with particles in presence of combustion, e.g., liquid droplet swarms for combustion temperature control, soot dynamics, or combustion-oriented particle image velocimetry.