I have a question about denoising fMRI data.
I know that some software or packages do that automatically.
So we can regress-out motion/wm/csf/global signals or use a_comp_cor and ica_aroma.
My exact question is about the order of scrubbing. Is that the last step of denoising? Can we do that at the first step? What reason is behind the decision?

In a lot of cases, motion outliers are added as one-hot encoded vectors into your matrix of all confounds (e.g., mean compartment signals, acompcor etc). So all that denoising happens simultaneously in a general linear model.

If we are talking scrubbing a la power et al 2012, volumes are actually suppressed.
Then it makes sense to start with that, and then regress other confounds.
The reason is that motion parameters for example are highly contaminated by motion peaks.
So if you are regressing them out before scrubbing, these particular time points drive the regression.
It makes more sense to regress out the residual smaller motion, than to regress out large motion effects which will get discarded anyway.

Now, as @Steven pointed out, lots of groups implement scrubbing as a big regression model where impulse binary vectors are used to model time points with excessive motions. In that case denoising and scrubbing are done simultaneously. This leaves you with weirdly interpolated time series, and I don’t think it’s as clean as previous method but you keep a regular sampling grid and you don’t have to change anything in your software to do scrubbing, as long as you can add regressors. Which is why I believe most software ended up with that approach.