You would have an answer to your question if you could classify pairs of elements in $GL_n(\mathbf F_q)$ up to simultaneous conjugacy. This is the notorious matrix pair problem, which is the quintessential wild problem in representation theory (see When is a classification problem "wild"? , Are wild problems related to undecidable ones? and How can classifying irreducible representations be a "wild" problem?).
So you are unlikely to get a nice answer for general $n$.
However, this may be doable (treating $q$ as a formal variable) for small $n$, maybe up to $n=4$; for each $n$, there are finitely many similarity class types of matrices in $GL_n(\mathbf F_q)$ (see http://www.sagemath.org/doc/reference/combinat/sage/combinat/similarity_class_type.html). A similarity class type is basically a union of conjugacy classes which have conjugate centralizers. To classify pairs of matrices up to similarity, you put the first one in canonical form, and then classify the conjugacy classes in its centralizer.
This has essentially been done for $n\leq 4$ in my paper with Singla and Spallone: Similarity of matrices over local rings of length two, at least the number of classes is known - see Remark 1.1 and Table 7.
In order to answer your question, a more delicate analysis will be needed, for going from pairs to subgroups.