Limitation on the shift theorem of DFT due to frequency resolution?
But what if I want to shift different frequency components by a different number of samples? In that case, you can't find FFT of the whole data set and multiply by the same exponential. You would have to treat different frequencies separately right?


What you're describing is a system with a non-constant (i.e. frequency-dependent) group delay.

Group delay is the derivative of phase, hence,  $\frac{\mathrm d\phantom f}{\mathrm df}\phi(f)$ of a system. 

In filter design, you (most of the time) want a constant group delay, so you design your filter with linear phase (which means you need to use FIR filters, and you need to make them symmetric); the derivative of a linear function is a constant!

Now, you want a given "group delay profile". Sooo, what about actually implementing the following idea:

Essentially, you approach this via the window design method that one typically uses to design linear-phase filters of a wanted amplitude response. You, instead, want a to define a group delay, and hence, a phase response over frequency.


Convert a definition of your group delay to a phase definition by numerically integrating that profile (i.e. do a cumulative sum). (the integration is the inverse operation to the differentiation that delivers the group delay from the phase response; so, the result should be that each value is a $e^{j2\varphi}$, with $2\pi\varphi$ being the relative delay you want). Do that definition at a high "oversampling", i.e. if you just need an actual "resolution" for your frequency-dependent delay of sampling rate / 10, define that function at 200 points. Note down that you have a "frequency-sampling-rate" of $r_\text{freq}=\frac{f_\text{sample}}{200}$.
Convolve the result with a low-pass linear filter impulse response. Why will become clearer in the next step.
Convert the result to time domain by doing an IDFT (IFFT). The effect of the low-pass filter above is now, a lot of the resulting values are (close enough to) zero. So, you can omit these. 


If I'm not mistaken, you now have a filter of desired group delay profile plus a constant (the constant stems from the low-pass filtering in frequency domain).



Another route if you only care about non-overlapping bands and these should have very different group delays, making it hard to properly build that filter via the filter design method:


Find the "hardest" (i.e. narrowest) bands; design a FIR band pass filter for these using your favorite math toolkit's filter design methodology (I propose designing an equivalent-width low pass and shifting it in frequency). The length of that filter will be the longest you'll need for all the band pass filters. Sufficiently over-sample, i.e. re-design that filter, but use let's say 20 times the length. Use that length for all other filters, too.
The group delay of a linear-phase phase filter is half its length. Now, all your bandpasses have the same group delay $\tau_\text{offset}$.
Adjust the group delay of the individual filters by transforming (FFT), multiplying by a complex sinusoid of appropriate frequency, and transforming back (IFFT). If your frequency isn't (close enough) to $f_s / \text{length}(filter)$, you'll get leakage / Gibb's phenomenon. You can trade accuracy for oscillations by applying a window. If your filters are long, it's probably not even necessary to window.
Sum up all the filter impulse responses.
(optional) if it turns out that a lot of the resulting taps are close to zero, truncate.