What is the most relevant choice for the factorization scale of the Drell-Yan process? Can it be the invariant mass at the reaction threshold?
The choice of factorisation scale for any process mimics the typical scales that are present in the problem, whether that be a hard scattering virtuality, mass, C.O.M energy or indeed transverse momentum. It may also be assigned the same value as the renormalisation scale at which the strong coupling is evaluated.
As far as I know, there is no definite answer as to what should be chosen as the factorisation scale - one chooses the scale such to incorporate higher order corrections through elimination of large logarithms that would otherwise spoil the perturbative expansion. The caveat is that with many scales in the problem, there is typically more than one type of logarithm so they can’t all be nullified. So, it becomes an optimisation exercise into what is the most efficient way of suppressing or softening the effect of large logarithms.
This is all under the premise of resummation and is the thing to do because an all order expansion (which is by no means achievable currently) will be void of all such artifacts ($\mu_R^2$, $\mu_F^2$ etc) arising in the choice of mathematical regularisation/factorisation framework.
There is a work here https://arxiv.org/abs/hep-ph/0703156v1 in which the authors find $\mu_F \sim Q/2$, where $Q$ is a final state invariant mass, to be an optimal choice of factorisation scale in the Drell-Yan process you mention.