We analyze the transverse momentum ( p T ) spectra of lepton pairs ( ℓ ℓ ¯ ) generated in the Drell-Yan process, as detected in proton-nucleus (pion-nucleus) and proton-(anti)proton collisions by ten collaborations over a center-of-mass energy s N N or s if in a simplified form) range from ~ 20 GeV to above 10 TeV. Three types of probability density functions (the convolution of two Lévy-Tsallis functions, the two-component Erlang distribution, and the convolution of two Hagedorn functions) are utilized to fit and analyze the p T spectra. The fit results are approximately in agreement with the collected experimental data. Consecutively, we obtained the variation law of related parameters as a function of s and invariant mass Q . In the fit procedure, a given Lévy-Tsallis (or Hagedorn) function can be regarded as the probability density function of transverse momenta contributed by a single quark ( q ) or anti-quark ( q ¯ ). The Drell-Yan process is then described by the statistical method.