There are several equivalent ways of implementing Dirichlet boundary conditions with the finite element method. I'll give a brief overview, but you'll probably need more details, which can be found in the book Understanding and Implementing the Finite Element Method.
It may be easier to analyze this for the Poisson problem $-\nabla^2u = f$ than for the heat equation; the techniques are the same in both cases, and if you understand them for the steady-state problem then it shouldn't be a big leap to adapt them to the time-dependent one. For the initial value, that amounts to setting $u^0_i = \mathbf{u}_0(x_i)$. You'll then use this as part of the right-hand side in the linear system you've written down when you go to solve for the next time-step $\mathbf{u}^1$. The boundary conditions on $\partial\Omega$ are trickier.
One way to implement Dirichlet boundary conditions is to alter the matrices $A$ and $B$ and the load vector in order to guarantee that the boundary conditions are satisfied exactly. Suppose that $i$ is some unknown of the linear system corresponding to a Dirichlet boundary point; the modification you have to make is to zero out row $i$ of $B$, set row $i$ of $A$ to be the corresponding row in the identity matrix and set $\mathbf{l}_i$ equal to $g(x_i)$.
To illustrate this approach, suppose you had a $3 \times 3$ chunk of the linear system that looked like this:
$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33} \end{matrix}\right]
\left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right]
= \left[\begin{matrix} l_1 \\ l_2 \\ l_3\end{matrix}\right]$
and that $u_2$ is a degree of freedom corresponding to a boundary node. The modified system would look like this:
$\left[\begin{matrix} a_{11} & a_{12} & a_{13} \\
0 & 1 & 0 \\
a_{31} & a_{32} & a_{33} \end{matrix}\right]
\left[\begin{matrix} u_1 \\ u_2 \\ u_3\end{matrix}\right]
= \left[\begin{matrix} l_1 \\ g_2 \\ l_3\end{matrix}\right]$
You can see how the modified system isn't symmetric anymore. Nonetheless, $u_1$ and $u_3$ are still coupled to $u_2$. The asymmetry of the system isn't usually an issue in practice provided that your initial guess for $u$ exactly satisfies the Dirichlet boundary conditions; try throwing the problem into a naive CG solver and see what happens. (Can you come up with a reason why this might be so?)
Another approach is to reduce the system size. Instead of solving for all of the coefficients $u_i$ in the expansion $\mathbf{u} = \sum_iu_i\phi_i$, you'll only solve for the coefficients corresponding to interior nodes. This means extracting sub-matrices out of $A$ and $B$ for those unknowns, and adding an extra component to the load vector $\mathbf{l}$ consisting of the effective forcing on the interior from boundary conditions. Referring back to the example matrices from before, we're left with the following reduced system for $u_1$ and $u_3$:
$\left[\begin{matrix} a_{11} & a_{13} \\
a_{31} & a_{33} \end{matrix}\right]
\left[\begin{matrix} u_1 \\ u_3\end{matrix}\right] =
\left[\begin{matrix} l_1 - a_{12}g_2 \\ l_3 - a_{32}g_2\end{matrix}\right]$
Changing the linear system or deflating it are, from what I've observed, the most common approaches. You can also add on boundary conditions later as a constraint, or you can use penalty methods, but these are less common. Of the two, I tend to favor the first approach.