``everybody knows'' that an integral orthogonal matrix is a signed permutation matrix, so there are exactly $2^n n!$ such matrices in $O(n).$ Now, what if we ask for the enumeration of elements of $O(n$ with rational elements, where all the denominators are equal to $q$ (where $q$ is an arbitrary integer, though if the question is much easier for $q$ prime, that's fine too...)? I would suspect this has been studied... (one can ask the same question for arbitrary algebraic groups/number fields, of course).

Boris Venkov asked this question, once, in the coffea room of the Maths Department of the university of Geneva. Here are some things that relate to the case $q=2^m$. Let $I_n$ be the 'euclidean' $\mathbf Z$module of rank $n$, and $A\mapsto G_n(A)$ the group it defines. Let us write $a(n,m)$ for the number of matrices in $G_n(\mathbf Z[\frac{1}{2}])$ all of whose coordinates have $2$valuation at least $m$. First note that $G_n(\mathbf Z[\frac{1}{2}])$ is finite for $n<5$. In fact, one even has $G_n(\mathbf Z[\frac{1}{2}])=G_n(\mathbf Z)\simeq\mathbf Z/2^n\rtimes S_n$ for $n<4$ and $G_4(\mathbf Z[\frac{1}{2}])=Aut(D_4)=G_4(\mathbf Z)\rtimes \mathbf Z/3$, where $D_4$ denotes the sublattice of elements of even length in $I_4$. The $\mathbf Z/3$ is generated by $$ \begin{bmatrix} 1/2 & 1/2 & 1/2 & 1/2\\\ 1/2 & 1/2 & 1/2 & 1/2\\\ 1/2 & 1/2 & 1/2 & 1/2\\\ 1/2 & 1/2 & 1/2 & 1/2 \end{bmatrix}$$ Thus in this case, you have $a(4,0)=384=2^7*3$ and $a(4,1)=1152=2^7*3^2$. Next, note that $G_5(\mathbf Z[\frac{1}{2}])$ is an amalgamated sum : $A\star_C B$ with $A=G_5(\mathbf Z)$, $B=Aut(I_1)\times Aut(D_4)$ and $C=A\cap B\simeq \mathbf Z/2\times G_4(\mathbf Z)$. This decomposition comes from the action of $G_5(\mathbf Z[\frac{1}{2}])$ on the BruhatTits Building of $G_5(\mathbf Q_2)$ : a tree whose vertices are 5valent and 3valent. The number $\frac{a(5,m)}{\vert G_5(\mathbf Z)\vert}$ counts the $5$valent vertices that are at (combinatorial) distance less than $2m$ from a fixed one. Thus $\frac{a(5,m)}{\vert G_5(\mathbf Z)\vert}=10*\frac{8^m1}{7}+1$. Things are becoming more complicated in higher dimensions since the BruhatTits building is not a tree anymore ($n\geq 6$) and the action of $G_n(\mathbf Z[\frac{1}{2}])$ on the vertices of a same type (say corresponding to unimodular lattices) is not transitive anymore ($n\geq 9$). You may adapt the same method to obtain similar results for small values of $n$ and $q=p^m$ : the BruhatTits building is a tree for $n=3$, the action of $G_3(\mathbf Z[\frac{1}{p}])$ is transitive on one type of vertices (but rarely on maximal simplices) ... 


This is not a complete answer, but it's a start. In the $2 \times 2$ case, the act of choosing the first column and clearing denominators describes a twotoone map from orthogonal matrices to primitive vectors in the lattice $\mathbb{Z} \oplus \mathbb{Z}$. In particular, there are twice as many orthogonal matrices with denominator exactly $q^2$ as there are primitive vectors of length $q$, and there are twice as many orthogonal matrices with denominator dividing $q^2$ as there are vectors of length $q$. The latter quantity is enumerated by the theta function of the lattice, which is a modular form of weight 1, and its coefficients grow logarithmically. In higher dimension, I guess the generating function comes from successively choosing vectors in orthogonal complements, so it should be at least related to modular forms. 

