Jets constructed via clustering algorithms (e.g., anti-$k_T$, soft-drop) have been proposed for many precision measurements, such as the strong coupling $alpha_s$ and the nucleon intrinsic dynamics. However, the theoretical accuracy is affected by missing QCD corrections at higher orders for the jet functions in the associated factorization theorems. Their calculation is complicated by the jet clustering procedure. In this work, we propose a method to evaluate jet functions at higher orders in QCD. The calculation involves the phase space sector decomposition with suitable soft subtractions. As a concrete example, we present the quark-jet function using the anti-$k_T$ algorithm with E-scheme recombination at next-to-next-to-leading order.