We report a calculation of the perturbative matching coefficients for the transverse-momentum-dependent parton distribution functions for quark at the next-to-next-to-next-to-leading order in QCD, which involves calculation of non-standard Feynman integrals with rapidity divergence. We introduce a set of generalized Integration-By-Parts equations, which allows an algorithmic evaluation of such integrals using the machinery of modern Feynman integral calculation.