In this paper, we develop a quantile functional regression modeling framework that models the distribution of a set of common repeated observations from a subject through the quantile function, which is regressed on a set of covariates to determine how these factors affect various aspects of the underlying subject-specific distribution. To account for smoothness in the quantile functions, we introduce custom basis functions we call textit{quantlets} that are sparse, regularized, near-lossless, and empirically defined, adapting to the features of a given data set and containing a Gaussian subspace so {non-Gaussianness} can be assessed. While these quantlets could be used within various functional regression frameworks, we build a Bayesian framework that uses nonlinear shrinkage of quantlet coefficients to regularize the functional regression coefficients and allows fully Bayesian inferences after fitting a Markov chain Monte Carlo. Specifically, we apply global tests to assess which covariates have any effect on the distribution at all, followed by local tests to identify at which specific quantiles the differences lie while adjusting for multiple testing, and to assess whether the covariate affects certain major aspects of the distribution, including location, scale, skewness, Gaussianness, or tails. If the difference lies in these commonly-used summaries, our approach can still detect them, but our systematic modeling strategy can also detect effects on other aspects of the distribution that might be missed if one restricted attention to pre-chosen summaries. We demonstrate the benefit of the basis space modeling through simulation studies, and illustrate the method using a biomedical imaging data set in which we relate the distribution of pixel intensities from a tumor image to various demographic, clinical, and genetic characteristics.