[guided]We want to prove more than set-theoretic existence: the orthonormal frames must vary smoothly with the base point. The right way to see this is to use an arbitrary smooth local frame and then correct it by a smooth positive-definite matrix.
Let $U \subset M$ be an open set over which $E$ has a smooth local frame $e_1,\dots,e_n$. This means that each $e_i: U \to E$ is a smooth section and that $e_1(x),\dots,e_n(x)$ is a basis of $E_x$ for every $x \in U$. The metric $g$ is smooth, so the matrix of its fibrewise inner product in this frame is the smooth map
\begin{align*}
G: U \to \operatorname{Sym}^+_n(\mathbb{R})
\end{align*}
defined by
\begin{align*}
G(x)_{ij} := g_x(e_i(x),e_j(x)).
\end{align*}
The target is $\operatorname{Sym}^+_n(\mathbb{R})$ because each $g_x$ is a positive-definite symmetric [bilinear form](/page/Bilinear%20Form).
The problem is now a finite-dimensional linear algebra problem depending smoothly on $x$: given the positive-definite matrix $G(x)$, find all matrices $B \in GL_n(\mathbb{R})$ such that the frame with columns $e_i(x)B_{ij}$ is orthonormal. That condition is exactly
\begin{align*}
B^\top G(x)B = I_n.
\end{align*}
Let $S: \operatorname{Sym}^+_n(\mathbb{R}) \to \operatorname{Sym}^+_n(\mathbb{R})$ denote the positive-definite square-root map, so $S(B)^2=B$ for every positive-definite matrix $B$. We use smooth functional calculus for symmetric matrices. Its hypotheses are satisfied here because $t \mapsto \sqrt{t}$ is smooth on the open set $(0,\infty)$ and every matrix in $\operatorname{Sym}^+_n(\mathbb{R})$ is symmetric with spectrum contained in $(0,\infty)$. Therefore $S$ depends smoothly on the positive-definite matrix. Define
\begin{align*}
R: U \to GL_n(\mathbb{R})
\end{align*}
by $R(x):=S(G(x))^{-1}$. Since $S(G(x))$ is symmetric and positive-definite, $R(x)^\top G(x)R(x)=I_n$.
Now take $A \in O(n)$ and define
\begin{align*}
u_{x,A}: \mathbb{R}^n \to E_x
\end{align*}
by
\begin{align*}
u_{x,A}(a) := \sum_{i=1}^n e_i(x)\,(R(x)A a)_i.
\end{align*}
This is a linear isomorphism because $R(x)A \in GL_n(\mathbb{R})$. For $a,b \in \mathbb{R}^n$, the matrix formula for the metric in the frame $e_1,\dots,e_n$ gives
\begin{align*}
g_x(u_{x,A}(a),u_{x,A}(b)) = a^\top A^\top R(x)^\top G(x)R(x)A b.
\end{align*}
Using $R(x)^\top G(x)R(x)=I_n$ and $A^\top A=I_n$, this becomes
\begin{align*}
g_x(u_{x,A}(a),u_{x,A}(b)) = a^\top b = \langle a,b\rangle_0.
\end{align*}
Thus $u_{x,A}$ is a $g$-orthonormal frame.
Conversely, suppose $u \in \operatorname{Fr}_g(E)_x$. Since $e_1(x),\dots,e_n(x)$ is a basis of $E_x$, there is a unique matrix $B \in GL_n(\mathbb{R})$ such that
\begin{align*}
u(a)=\sum_{i=1}^n e_i(x)(Ba)_i.
\end{align*}
The assumption that $u$ is $g$-orthonormal says $B^\top G(x)B=I_n$. Multiplying this identity by $S(G(x))$ on the left and right shows that $S(G(x))B$ is orthogonal. If we define $A:=S(G(x))B$, then $A \in O(n)$ and $B=R(x)A$. Therefore every orthonormal frame over $U$ has the form $u_{x,A}$ for a unique $A \in O(n)$.
This gives a local trivialization
\begin{align*}
U \times O(n) \to \operatorname{Fr}_g(E)|_U, \quad (x,A) \mapsto u_{x,A}.
\end{align*}
The map is smooth because $e_i$, $G$, the square-root map $S$, matrix inversion, and matrix multiplication are smooth on their domains. Its inverse is smooth because it is obtained by writing the frame in the smooth local frame $e_1,\dots,e_n$ to get $B$, then forming $A=S(G(x))B$. The right action of $O(n)$ corresponds to multiplication in the second factor, so this local model proves that $\operatorname{Fr}_g(E)$ is a smooth principal $O(n)$-subbundle of $\operatorname{Fr}(E)$.[/guided]