In a Galerkin approach for stationary PDEs, you write your unknown solution as a linear combination of "nice" basis functions $\{\phi_j\}_{j=1,\dots,n}$,
$$u(x) \approx \sum_{j=1}^n u_j \phi_j(x),$$
and project the equation onto the space spanned by (these or other) functions, exploiting the linearity of the differential operators and the inner product. In your case, you get
$$\sum_{j=1}^n u_j (c\nabla^2 \phi_j(x)+a\nabla \phi_j(x) + b,\psi_i(x))=0\qquad\text{for all }1\leq i\leq n.$$
This yields a set of linear equations for the unknown coefficients $u_j$.
In the finite element method, boundary conditions are implemented differently for Dirichlet and for Neumann conditions. Homogeneous Dirichlet conditions are built into the choice of the basis: $\phi_j(0)=\psi_i(0)=0$ for all $1\leq i,j\leq n$. The Neumann conditions, on the other hand, become part of the operator:
You integrate by parts in the first term of the sum to obtain (I am leaving out the lower order terms since they do not play such a role here)
$$\sum_{j=1}^n u_j \left(\phi_j'(1)\psi_i(1)-(c\nabla \phi_j(x),\nabla\psi_i(x)) \right)\qquad\text{for all }1\leq i\leq n$$
(using that all boundary terms for $x=0$ vanish due to the choice of $\psi$).
Of course, then you use that
$$\sum_{j=1}^n u_j \phi_j'(1)\psi_i(1) = u'(1)\psi_i(1)$$
so that becomes part of your right-hand side. In your case, $u'(1)=0$ so the whole term vanishes.
In the method of lines, you approximate your unknown solution $u(t,x)$ as
$$u(t,x)\approx \sum_{j=1}^n u_j(t) \phi_j(x),$$
i.e., the whole space dependence is in the basis functions, while the whole time dependence is in the scalar coefficients. Your projected equation is now
$$\sum_{j=1}^n u_j(t) (c\nabla^2 \phi_j(x)+a\nabla \phi_j(x) + b,\psi_i(x)) = \sum_{j=1}^n \frac{d}{dt}u_j(t)(\phi_j(x),\psi_i(x))\qquad\text{for all }1\leq i\leq n.$$
Now you do the same steps; since the inner products on the left-hand side are exactly the same as before, you get exactly the same matrices. On the right-hand side, you have a term just like a standard mass matrix $M$ with entries $M_{ij}=(\phi_j(x),\psi_i(x))$, where you respect the Dirichlet conditions by the choice of basis functions, but you don't have to worry about Neumann conditions. Note that the time dependence is purely in the coefficients, so for the (spatial) mass matrix it makes no difference whether this comes from approximating $u$ or $\partial_t u$.
The way inhomogeneous Dirichlet conditions are treated mathematically is to reduce the problem to the homogeneous case: You write
$$u(t,x) = g(t,x) + \tilde u(t,x)$$
for some fixed $g$ with $g(t,0)=1$ and $g'(t,1)=0$. Then, $\tilde u$ satisfies
$$c\nabla^2 \tilde u + a\nabla \tilde u = -b -c\nabla^2 g -a\nabla g,\qquad \tilde u(0) = \tilde u'(1) = 0,$$
and you can proceed as above to compute $\tilde u(t,x)$ and then $u(t,x) = \tilde u(t,x)+g(t,x)$. If you work with the standard hat functions, you could simply take $g(t,x) = \phi_1(x)$.
If the initial condition (which you have not stated) satisfies $u_0(0)=1$, you could also replace the first equation of the system of ODEs by $\partial_t u_1 = 0$, i.e., replace the first row in the mass matrix by $(1,0,\dots,0)$, the first row of the remaining matrices by $(0,\dots,0)$, and the first component of the right-hand side by $0$.
If memory serves correctly, the method of lines is discussed in Ern, Guermond: Theory and Practice of Finite Elements, Springer, 2004, Chapter 6.1.
