In this paper, we give an estimate of sub-Laplacian of Riemannian distance functions in pseudo-Hermitian geometry which plays a similar role as Laplacian comparison theorem in Riemannian geometry, and deduce a prior horizontal gradient estimate of pseudo-harmonic maps from pseudo-Hermitian manifolds to regular balls of Riemannian manifolds. As an application, Liouville theorem is established under the conditions of nonnegative pseudo-Hermitian Ricci curvature and vanishing pseudo-Hermitian torsion. Moreover, we obtain the existence of pseudo-harmonic maps from complete noncompact pseudo-Hermitian manifolds to regular balls of Riemannian manifolds.