If you take a matrix of standard normal draws and apply a QR decomposition, the distribution induced over Q will be uniform over all orthogonal matrices. You can get both rotations (determinant = 1) or reflections (determinant = -1), so if you really want to stick to the special orthogonal manifold of just rotations, you have to somehow eliminate the determinant = -1 alternatives.
If you then try to impose a different distribution on top of that uniform, you run into two problems. First, you have to do a change of variables adjustment before you can apply the different distribution, and second, you induce a lot of constraint geometry and dependence through the orthonormal constraint (e.g., if one component gets bigger, the sum of the squares of the other ones has to go down by an equal amount, inducing a lot of correlation).
Alternatively, if you start with a skew symmetric matrix with random normal entries and put that through a matrix exponential, you also get uniformly distributed orthonormal matrices. But it has the same problems and I believe is not as stable arithmetically as QR.