Decomposition of a single 4D rotation I have a $4\times 4$ matrix $M$ which represents a general 4-dimensional rotation.
$$
M = \pmatrix{a_{11} &amp;a_{12} &amp;a_{13} &amp;a_{14}\\a_{21} &amp;a_{22} &amp;a_{23} &amp;a_{24}\\a_{31} &amp;a_{32} &amp;a_{33} &amp;a_{34}\\a_{41} &amp;a_{42} &amp;a_{43} &amp;a_{44}} 
$$
There are also six 4D basis rotation matrices $M_{XY}$, $M_{YZ}$, $M_{ZX}$, $M_{XW}$, $M_{YW}$, $M_{ZW}$.

How to decompose 4D rotation $M$ into basis rotations (i.e., to find corresponding angles $T_{XY}$, $T_{YZ}$, $T_{ZX}$, $T_{XW}$, $T_{YW}$, $T_{ZW}$)?

Wikipedia says that any 4D rotation given by its matrix can be decomposed into two isoclinic rotations. However, I do not know if it is possible to decompose them further into basis rotations. Also, there may be another (simpler) way.
You should be able to build your basis matrices by doing something similar to the Gram-Schmidt process: by multiplying by suitably-chosen rotation matrices you can set each of the off-diagonal coefficients in turn equal to 0.  For instance, consider applying a $M_{ZW}$ rotation by some angle $\theta$ to the matrix:
$$M' = M_{ZW}(\theta)\ M = \pmatrix{1&amp;0&amp;0&amp;0\\0&amp;1&amp;0&amp;0\\0&amp;0&amp;\cos\theta&amp;-\sin\theta\\0&amp;0&amp;\sin\theta&amp;\cos\theta}\pmatrix{a_{11}&amp;a_{12}&amp;a_{13}&amp;a_{14}\\a_{21}&amp;a_{22}&amp;a_{23}&amp;a_{24}\\a_{31}&amp;a_{32}&amp;a_{33}&amp;a_{34}\\a_{41}&amp;a_{42}&amp;a_{43}&amp;a_{44}}
=\pmatrix{a_{11}&amp;a_{12}&amp;a_{13}&amp;a_{14}\\a_{21}&amp;a_{22}&amp;a_{23}&amp;a_{24}\\*&amp;*&amp;a_{33}\cos\theta-a_{43}\sin\theta&amp;a_{34}\cos\theta-a_{44}\sin\theta\\*&amp;*&amp;a_{33}\sin\theta+a_{43}\cos\theta&amp;a_{34}\sin\theta+a_{44}\cos\theta}$$
(Here the *s are 'don't-cares' - they have a value, but it's essentially unimportant.)  Now, by choosing $\theta$ correctly, you can make $M'_{34} = a_{34}\cos\theta-a_{44}\sin\theta$ equal to zero without affecting any other upper-triangular elements.  Once you've got that matrix, you can continue in a similar fashion by finding a suitable $M_{YW}(\phi)$ to zero the $M_{24}$ element, then an $M_{XW}$ to zero the $M_{14}$ element, etc.; repeating this procedure for each combination of unequal coordinates (in some canonically-defined order) will eventually lead you to an orthogonal (by hypothesis) lower-triangular matrix, which must therefore be the identity.
You can think of this process in terms of applying successive coordinate-plane rotations, each one on the right of what has gone before, so that ultimately the resulting rotation is the identity.  

Then multiplying the identity matrix by the inverses of these coordinate-plane rotations, each one also on the right of what has gone before, you will end up with an equation showing how to express the original rotation as a product of coordinate-plane rotations.

This actually works in any dimension. Since the dimension of SO(n) is n-choose-2 = n(n-1)/2, that should be the maximum number of coordinate-plane rotations needed.

In a given dimension, selecting a fixed ordering for performing this procedure defines a continuous map from the n(n-1)/2 dimensional torus (since each angle may be thought of as a point on circle) to SO(n+1):

Tn(n-1)/2  —>  SO(n)

that is in fact surjective (since each rotation can be expressed as a product of coordinate-plane rotations).  The circle factors of Tn(n-1)/2 are indexed by pairs $(i, j)$ with $1 \le i \lt j \le n$ (since each one represents one coordinate plane).

The 2-plane rotations used in this procedure are sometimes called Householder rotations.