Multizeta function values
The elements of $S$ are conjectured to be $\mathbb{Q}$-linearly independent, and so a basis for the $\mathbb{Q}$-linear span of the multiple zeta values.
This is what Francis Brown accomplished at the motivic level. Actually, his proof via motives goes through the direction that is currently missing on the level of numbers, which is why it doesn't lead to an algorithm for writing down the (conjecturally unique) expression. In Theorem 7.4 of [Mixed Tate motives over $\mathbb{Z}$ (Ann. Math., 2012)], he proved that the motivic enrichments $\zeta^{\mathbb{m}}(\mathbf{s})$, $\mathbf{s} \in \{2,3\}^{\times}$ of the $\{2,3\}^{\times}$-zeta values are $\mathbb{Q}$-linearly independent. The dimension $d_n$ of the $\mathbb{Q}$-linear span $\mathcal{M}_n$ of the motivic MZV of a given weight $n$ is well known to be equal to the $z^n$ coefficient of $1/(1-z^2-z^3)$, which is also the number of expressions of $n$ as a sum of $\{2,3\}$-elements, and so the number of $\{2,3\}^{\times}$ zeta values of weight $n$. The linear independence of the $\zeta^{\mathbb{m}}(\mathbf{s}), \, \mathbf{s} \in \{2,3\}^{\times}$ then implies the equality: these elements form a basis for $\bigoplus_n \mathcal{M}_n$. There is a realization homomorphism $\mathcal{M} \to \mathbb{R}$ sending $\zeta^{\mathbb{m}}(\mathbf{s}) \to \zeta(\mathbf{s})$, and it follows (non-constructively) that there is a representation on the level of numbers.
However, proving $\mathbb{Q}$-linear independence remains completely out of reach on the level of numbers. We know very little here: even the irrationality of $\zeta(5)$ has remained an unsolved problem (for instance, see this MO question and answer).
What does follow immediately from the motivic formalism is the upper bound $\dim{M_n} \leq d_n$ of the conjectured dimension of the $\mathbb{Q}$-linear span of the weight $n$ MZV. This was proved independently by Terasoma and Goncharov; it means that there are many linear relations among MZV of a given (high enough) weight. Given that, identities like the weight $12$ one you quote should not come as a surprise. Explicit such identities were obtained by Gangl, Kaneko and Zagier. They are linked to modular forms of weight $n$, via Eisenstein series.
Vesselin has already given an excellent answer to your first question, so let me just elaborate a little bit on the relation
$$ 28\zeta(3,9)+150\zeta(5,7)+168\zeta(7,5)=\frac{5197}{691}\zeta(12)\tag{1} $$
and its connection with modular forms for $\operatorname{SL}_2(\mathbb{Z})$. First of all, it is a rather exotic and surprising relation (I will qualify this statement below). In a nutshell, the coefficients $28,150,168$ on the left-hand side of (1) come from period integrals of the cusp form $\Delta$, while the denominator $691$ of the right-hand side comes from the constant term of the Eisenstein series $G_{12}$.
To any cuspidal Hecke eigenform $f$ of weight $k$, one can associate an even (i.e. invariant under $Y\mapsto -Y$) homogeneous polynomial $r_f(X,Y)^{+} \in K_f[X,Y]_{k-2}$ of degree $k-2$, where $K_f$ is the number field generated by the Fourier coefficients of $f$. Explicitly,
$$ r_f(X,Y)^{+}=\frac{\operatorname{Re}\left( (2\pi i)^{k-1}\int_0^{i\infty}f(z)(X-zY)^{k-2}dz \right)}{\omega_f^{+}} $$ where the integral is along the geodesic path from $0$ to $i\infty$ on the upper half-plane, and $\omega_f^{+} \in \mathbb{R}$ is a suitable ''period''. This construction is part of the very well-known Eichler--Shimura--Manin theory, expositions of which can be found in many places; see e.g. Lang's Introduction to Modular Forms for an introduction, or Brown's Multiple Modular Values and the relative completion of the fundamental group of $M_{1,1}$, Section 7 for a more condensed presentation.
Now in weight $12$ there is a unique cuspidal Hecke eigenform, namely Ramanujan's $\Delta$, and one has (up to scaling by an element of $\mathbb{Q}^{\times}$)
$$ r_{\Delta}(X,Y)^+=\frac{36}{691}(X^{10}-Y^{10})-X^2Y^2(X^2-Y^2)^3. $$
The essential part is the second summand $p(X,Y):=X^2Y^2(X^2-Y^2)^3$, the first one is a coboundary in the sense of group cohomology for $\operatorname{SL}_2(\mathbb{Z})$, and hence can be ignored (for the time being, but see below). Now Gangl--Kaneko--Zagier's result (Theorem 3 in their paper) says in essence that if numbers $q_{r,s}$ are defined by $p(X+Y,Y)=\sum_{r,s \geq 1}\binom{r+s-2}{r-1}q_{r,s}X^{r-1}Y^{s-1}$, then
$$ \sum_{r+s=12, \, r,s \, odd}q_{r,s}\zeta(r,s) \equiv 0 \mod \zeta(12) $$
which, after rescaling, yields (1) modulo $\zeta(12)$. The upshot is that the coefficients $28,150,168$ come from integrals involving the cusp form $\Delta$, and therefore the left-hand side of (1) has a rather clean, arithmetic-geometric interpretation.
To the best of my knowledge, no such interpretation is known for the right-hand prefactor $\frac{5197}{691}$ (in Gangl--Kaneko--Zagier's original approach, they need to compute some kind of beta-integral). However, the occurrence of the denominator $691$ seems to be an artifact of the well-known congruence $G_{12} \equiv \Delta \mod 691$ on the level of Fourier expansions, where $G_{12}$ is the Hecke Eisenstein series of weight $12$ (basically because the numerator of the constant term of G_{12} is divisible by 691). The precise relation should be explained by recent work of Hain--Matsumoto on Universal Mixed Elliptic Motives (in particular, the previously neglected term $\frac{36}{691}(X^{10}-Y^{10})$ resurfaces as the dual of the cocycle of $G_{12}$). As far as I can see, this does not give an explanation about the numerator $5197$, though.
Finally, I would like to point out that there is another connection of (1) to cusp forms, which is not (yet?) arithmetic-geometric but in my opinion very nice: namely, via the theory of double Eisenstein series, which are to usual Eisenstein series what double zeta values are to single zeta values. Since my answer is already much too long, I refer to the original paper of Gangl--Kaneko--Zagier, Section 7 for the details. Some of this has also been generalized by Henrik Bachmann in his PhD thesis (see in particular the bottom of p.34 for details on how to obtain (1) via this theory).
So why is the relation (1) surprising? It is more generally believed that all algebraic relations between multiple zeta values are implied by the (regularized) double shuffle equations, which are a class of quadratic relations between multiple zeta values. Specialized to the depth two case, they read
$$ \zeta(r)\zeta(s)=\zeta(r,s)+\zeta(s,r)+\zeta(r+s)=\sum_{m+n=r+s}\left( \binom{m-1}{r-1}+\binom{m-1}{s-1}\right)\zeta(m,n). \tag{2} $$ One of the main results of Gangl--Kaneko--Zagier is that (1) cannot be deduced from (2), and one needs to resort to the double shuffle equations in higher depth. If one were able to deduce from (2) that $28\zeta(4)\zeta(8)+\frac{95}{3}\zeta(6)^2 \in \mathbb{Q}\zeta(12)$ (which is of course a consequence of Euler's theorem on even zeta values), this ''depth-defect'' would go away, alas one cannot. This is one of the very perplexing features of the depth filtration on multiple zeta values, and that it can be resolved using modular forms points at some deep relation between the two objects, which we do not completely understand yet.