You need to be careful with the way you handle unconstrained sources in the context of statistical analysis. See the tutorial Difference: https://neuroimage.usc.edu/brainstorm/Tutorials/Difference#Unconstrained_sources
Now, do I compute within-subject source average first and then run Time-frequency process?
No, never compute time-frequency decompositions on averages. Average the TF of single trials:
https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#MEG_recordings:_Single_trials
Or do I take the sources directly and compute time frequency to obtain one averaged source time frequency file, for each subject?
Source maps (links) of individual trials => Time-frequency => Average (can be the online average of the Morlet wavelets or Hilbert processes, it avoids creating hundreds of useless files).
One option is to process only some specific scouts instead of the full source maps, in order to keep the size of the data generated (and the number of multiple comparisons to correct for!) reasonable.
After that, I'll project on default anatomy right?
If processing full source maps: Projecting to the template and then doing the TF analysis might be easier to handle, but I'm not sure. If you want to project on the template the TF values, you need to keep the complex values instead of saving only the power.
Both approaches should give the same results, the two operations are permutable.
If you use TF of scouts: no need to project to a standard space, you can directly compare the values between subjects.
And then what pipeline do I follow to perform statistics between two conditions?
This can be tricky, especially with unconstrained sources.
Non-parametric permutation tests should work with any kind of data in input, but the effects of interest still have to present in the data you are testing.