This is a matrix algebra heavy solution. I think there might be a better way to do this, but here goes.
Take from Denis Andrienko's notes:
$$Q_{ij} = \frac{1}{N}\sum_\alpha(u_i^\alpha u_j^\alpha - \frac{1}{3}\delta_{ij})$$
where $\mathbf{u}^\alpha$ is the unit vector pointing along the long axis of molecule $\alpha$ located at $\mathbf{x}^\alpha$, and $N$ the number of molecules.
Or if you prefer matrix notation,
$$Q = \frac{1}{N}\sum_\alpha((\mathbf{u}^\alpha)(\mathbf{u}^\alpha)^T - \frac{1}{3}I)$$
We write out $\langle Q \rangle$ in an orthonormal basis $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$, where one of the axes is set by the director $\mathbf{e}_1 = \langle\sum_\alpha \mathbf{u}^\alpha\rangle \big/ \|\langle\sum_\alpha \mathbf{u}^\alpha\rangle\|$:
$$\langle Q\rangle = U^T \frac{1}{N}\sum_\alpha \langle (\mathbf{u}^\alpha)(\mathbf{u}^\alpha)^T \rangle U - \frac{1}{3}I$$
where the matrix $U = [\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3]$. This can be written as
$$\langle Q\rangle = \frac{1}{N}\sum_\alpha \langle ((\mathbf{u}^\alpha)^TU)^T((\mathbf{u}^\alpha)^TU) \rangle - \frac{1}{3}I$$
Writing this out in its full glory:
$$\langle Q\rangle =
\frac{1}{N}\sum_\alpha
\begin{pmatrix}
\langle((\mathbf{u}^\alpha)^T\mathbf{e}_1)^2\rangle & \langle(\mathbf{u}^\alpha)^T\mathbf{e}_1(\mathbf{u}^\alpha)^T\mathbf{e}_2\rangle & \langle(\mathbf{u}^\alpha)^T\mathbf{e}_1(\mathbf{u}^\alpha)^T\mathbf{e}_3\rangle \\
\langle(\mathbf{u}^\alpha)^T\mathbf{e}_1(\mathbf{u}^\alpha)^T\mathbf{e}_2\rangle & \langle((\mathbf{u}^\alpha)^T\mathbf{e}_2)^2\rangle & \langle(\mathbf{u}^\alpha)^T\mathbf{e}_2(\mathbf{u}^\alpha)^T\mathbf{e}_3\rangle \\
\langle(\mathbf{u}^\alpha)^T\mathbf{e}_1(\mathbf{u}^\alpha)^T\mathbf{e}_3\rangle & \langle(\mathbf{u}^\alpha)^T\mathbf{e}_2(\mathbf{u}^\alpha)^T\mathbf{e}_3\rangle & \langle((\mathbf{u}^\alpha)^T\mathbf{e}_3)^2\rangle
\end{pmatrix}
- \frac{1}{3}I
$$
Writing $\mathbf{e}_1$ out, we find that $\langle((\mathbf{u}^\alpha)^T\mathbf{e}_1)^2\rangle$ = $\langle\cos^2\theta^\alpha\rangle$. If the molecule is uniaxial, we have $\langle((\mathbf{u}^\alpha)^T\mathbf{e}_2)^2\rangle = \langle((\mathbf{u}^\alpha)^T\mathbf{e}_3)^2\rangle = \langle\frac{1}{2}\cos^2(\frac{\pi}{2}-\theta^\alpha)\rangle$ (if its not uniaxial, one gets more weight while the other the same amount less). The offdiagonals vanish and the sum is just an average over all molecules.
Thus we have
$$
\langle Q \rangle =
\begin{pmatrix}
\frac{2}{3}S & 0 & 0 \\
0 & -\frac{1}{3}S & 0 \\
0 & 0 & -\frac{1}{3}S
\end{pmatrix}
$$
where $S = \frac{1}{2}\langle 3\cos^2\theta^\alpha - 1 \rangle$.
That is to say that we can write
$$\langle Q_{ij}\rangle = S(n_i n_j - \frac{1}{3}\delta_{ij})$$
where the unit vector $\mathbf{n}$ specifies the direction of the principal axis of $\langle Q_{ij}\rangle$ in the transformed coordinates (e.g. $\textbf{n} = (1, 0, 0)$ points towards $\textbf{e}_1$ etc.).
