Based on thermodynamics, we discuss the galactic clustering of expanding Universe by assuming the gravitational interaction through the modified Newtons potential given by $f(R)$ gravity. We compute the corrected $N$-particle partition function analytically. The corrected partition function leads to more exact equations of states of the system. By assuming that system follows quasi-equilibrium, we derive the exact distribution function which exhibits the $f(R)$ correction. Moreover, we evaluate the critical temperature and discuss the stability of the system. We observe the effects of correction of $f(R)$ gravity on the power law behavior of particle-particle correlation function also. In order to check feasibility of an $f(R)$ gravity approach to the clustering of galaxies, we compare our results with an observational galaxy cluster catalog.