Since you only impose mild assumptions on $f$, $g$, the proof is somewhat technical, e.g. we cannot work with Riemann sums because $f^2$ might not be Riemann integrable.
Without loss of generality, I will assume that $a=0$ and $X_0=0$. Write $X_t = M_t+A_t$ where $$M_t := \int_0^t f(s) \, dB_s \qquad A_t := \int_0^t g(s) \, ds.$$ If we denote by $\langle \cdot,\cdot \rangle$ the quadratic (co)variation, then $$\langle X,X \rangle_t = \langle M+A,M+A \rangle_t = \langle M,M \rangle_t + 2 \langle M,A \rangle_t + \langle A,A \rangle_t. \tag{1}$$ This follows by a straight-forward computation similar to that in your question. We are going to show that \begin{align*} \langle M,M \rangle_t &= \int_0^t f(s)^2 \, ds \tag{2} \\ \langle M,A \rangle_t &= 0 \tag{3} \\ \langle A,A \rangle_t &= 0. \tag{4} \end{align*}
Proof of $(4)$:
Let $g=g(t,\omega)$ be a measurable function such that $g(\cdot,\omega) \in L^1([0,T])$ for any $T>0$. Then $t \mapsto A_t(\omega) = \int_0^t g(s,\omega) \, ds$ is a continuous function for each $\omega$ and so $A_{\bullet}(\omega)$ is uniformly continuous on $[0,T]$. If $\Pi=\{0=t_0<\ldots<t_n=T\}$ is a partition of $[0,T]$ with mesh size $|\Pi|$, then \begin{align*} \sum_i |A_{t_{i+1}}-A_{t_i}|^2 &\leq \sup_{|s-t| \leq |\Pi|, s,t \in [0,T]} |A_{s}-A_t| \sum_{i=1}^n |A_{t_{i+1}}-A_{t_i}| \\ &\leq \sup_{|s-t| \leq |\Pi|, s,t \in [0,T]} |A_{s}-A_t| \int_0^T |g(s)| \, ds. \end{align*} Because of the uniform continuity on $[0,T]$, the right-hand side converges a.s. to $0$ as the mesh size $|\Pi|$ tends to zero. This proves $\langle A,A \rangle_T=0$.
Proof of $(3)$:
This is quite similar to the previous proof. Take measurable $f,g$ such that $f(\cdot,\omega) \in L^2([0,T])$ and $g(\cdot,\omega) \in L^1([0,T])$ for $T>0$. The stochastic integral $M_t = \int_0^t f(s) \, dB_s$ has continuous sample paths with probability $1$. Exactly as in the previous part, we get $$\sum_i |M_{t_{i+1}}-M_{t_i}| \, |A_{t_{i+1}}-A_{t_i}| \leq \sup_{|s-t| \leq |\Pi|, s,t \in [0,T]} |M_s-M_t| \int_0^T |g(s)| \,ds.$$ Because of the uniform continuity on compact time intervals, the right-hand side converges to $0$ as $|\Pi| \to 0$. Hence, $\langle M,A \rangle_T=0$ for all $T>0$.
Proof of $(2)$:
For simple functions this is a straight-forward calculation, see this question. To extend $(2)$ to a larger class of functions, we need to use approximation techniques. For brevity of notation set $$S_{\Pi}(Y,Z) := \sum_{i} (Y_{t_{i+1}}-Y_{t_i})(Z_{t_{i+1}}-Z_{t_i})$$ and $S_{\Pi}(Y) =: S_{\Pi}(Y,Y)$.
Case 1: $f$ satisfies $\mathbb{E}\int_0^T f(s)^2 \, ds < \infty$ for each $T>0$.
Since $f$ is progressively measurable and satisfies the above integrability condition, there exists a sequence of simple functions $(f_n)_{n \in \mathbb{N}}$ such that $$\mathbb{E}\int_0^T |f(s)-f_n(s)|^2 \, ds \to 0, \qquad T>0 \tag{5}$$ and $$\mathbb{E} \left| \int_0^T f_n(s) \, dB_s - \int_0^T f(s) \, dB_s \right|^2 \to 0, \qquad T>0. \tag{6}$$ Set $M_n(t):=\int_0^t f_n(s) \, dB_s$ and fix $T>0$. We have $$\langle M,M \rangle_T = \langle M-M_n,M-M_n \rangle_T + 2 \langle M-M_n,M_n \rangle_T + \langle M_n,M_n \rangle_T. \tag{7}$$ Let $\Pi$ be a partition of $[0,T]$ . Taking expectation and applying Itô's isometry, we find \begin{align*} \mathbb{E}(S_{\Pi}(M-M_n)) &= \sum_i \mathbb{E}\int_{t_i}^{t_{i+1}} (f_n(s)-f(s))^2 \, ds \\ &= \mathbb{E}\int_0^T (f_n(s)-f(s))^2 \, ds \end{align*} Letting $|\Pi| \to 0$ using Fatou's lemma, we get
$$\mathbb{E}(\langle M-M_n \rangle_T) \leq \mathbb{E}\int_0^T (f_n(s)-f(s))^2 \, ds \xrightarrow[n \to \infty]{(5)} 0, $$ which shows that $\langle M-M_n\rangle_T \to 0$ in $L^1$. Similarily, an application of Itô's isometry (combined with the polarization identity, see here) shows that $$\mathbb{E}(S_{\Pi}(M-M_n,M_n)) = \mathbb{E}\int_0^T (f_n(s)-f(s)) f_n(s) \, ds.$$ Applying the Cauchy-Schwarz inequality and using $(5)$, it follows that the right-hand side converges to $0$ as $n \to \infty$ (uniformly in $\Pi$), and so $\langle M-M_n,M_n \rangle_T \to 0$ in $L^1$. Finally we already know that $\langle M_n,M_n \rangle_T=\int_0^T f_n(s)^2$ and so $$\lim_{n \to \infty} \langle M_n,M_n \rangle_T = \int_0^T f(s)^2 \, ds.$$ Letting $n \to \infty$ in $(7)$ proves the assertion.
Case 2: $\int_0^t f(s)^2 \, ds < \infty$ with probability $1$.
In order to extend $(2)$ such function we need to truncate $f$, e.g. consider $f_n := (-n) \vee f \wedge n$. Each $f_n$ satisfies the integrability assumption from Case 1, and so we know the quadratic variation of $\int_0^t f_n(s) \, dB_s$. Now, similar as in the previous part, we can use this knowledge to compute the quadratic varion of $\int_0^T f(s) \, dB_s$. Let me know in case that you really want to see all the details.
