Quarkonium production in hadron collisions at low transverse momentum q(perpendicular to) << M with M as the quarkonium mass can be used for probing transverse momentum dependent (TMD) gluon distributions. For this purpose, one needs to establish the TMD factorization for the process. We examine the factorization at the one-loop level for the production of eta(c) or eta(b). The perturbative coefficient in the factorization is determined at one-loop accuracy. Comparing the factorization derived at tree level and that beyond the tree level, a soft factor is, in general, needed to completely cancel soft divergences. We have also discussed possible complications of TMD factorization of p-wave quarkonium production.