When you calculate work, you do so along a given path. Here, that path has tangent vector $d\mathbf s$. This is a vector with direction; the minus sign will ultimately come from choosing the path's orientation--inward or outward.
Edit: Aha, I think I've found the unintuitive part. The key is in the use of the coordinate $r$ to parameterize the path, in that $r$ is larger at the start of the path and smaller at the end. This runs counter to what you would usually do when parameterizing such a path with an arbitrary parameter.
Let $\mathbf s_0$ and $\mathbf s_1$ be the starting and ending points of a path $\mathbf s(\lambda) = \mathbf s_0 + (\mathbf s_1 -\mathbf s_0)\lambda$. The work integral is then
$$W = \int_{\mathbf s_0}^{\mathbf s_1} \mathbf F(\mathbf s) \cdot d\mathbf s= \int_0^1 \mathbf F(\mathbf s(\lambda)) \cdot \frac{d\mathbf s}{d\lambda} \, d\lambda = \int_0^1 \mathbf F (\mathbf s(\lambda))\cdot (\mathbf s_1 - \mathbf s_0) \, d\lambda$$
For two finite points, the basic approach is sound, but it breaks down when you have a point at infinity involved. This is the reason that the problem of assembling a configuration is usually attacked with a different basic parameterization.
Instead, set $\mathbf s(\lambda) = \lambda \hat{\mathbf a}$ for some unit vector $\hat{\mathbf a}$ and set the bounds of the integral as being from $[\infty, R)$. This is the important point: even though the path is being traversed coming in from infinity, the parameterization means that $d\mathbf s/d\lambda = + \hat{\mathbf a}$, not minus as I originally thought. The path's still oriented outward; we're just traversing it backwards.
Here's how that integral looks:
$$W = \int_{\infty}^R \mathbf F(\lambda \hat{\mathbf a}) \cdot \hat{\mathbf a} \, d\lambda$$
Of course, we know the expression for the electric force:
$$\mathbf F(\mathbf r) = k\frac{qq_0 \mathbf r}{|r|^3}$$
Plug in $\mathbf r = \lambda \hat{\mathbf a}$ to get
$$\mathbf F(\lambda \hat{\mathbf a}) = k \frac{qq_0 \lambda \hat{\mathbf a}}{\lambda^3} = k \frac{q q_0 \hat{\mathbf a}}{\lambda^2}$$
We find that the integrand is then
$$W = \int_{\infty}^R k \frac{q q_0}{\lambda^2} \hat{\mathbf a} \cdot \hat{\mathbf a} \, d\lambda = \int_\infty^R k \frac{q q_0}{\lambda^2} \, d\lambda = - k \frac{q q_0}{R} < 0$$
The work is negative, so the change in potential energy $\Delta U = - W$ is positive as required.
So where is the problem then? As we've seen, there actually shouldn't be an extra negative sign coming in on line 4 (as posted in the OP's question). This is somewhat obscured because an explicit parameterization of the path is never written down in the first place--usually, you don't have to, but this problem is tricky enough that it helps immensely.
Best Answer
In mechanics, we love energy: indeed, it often really simplifies solving problems. Often, we want to know how the speed of a particle varies, or, in other words, what its kinetic energy is. Luckily, for some forces, knowing the variation in kinetic energy is really easy: these forces are the conservative forces. A force $\textbf{F}$ is conservative when its work along a path $\Gamma$ only depends on the initial and final positions: then, we can introduce a potential energy, say $E_p$, so that
$$ W(\textbf{F}) = \int_\Gamma \textbf{F}\cdot\textbf{dl} = E_{pi} - E_{pf}$$
However, the kinetic energy theorem states that the total work done on a particle equals its kinetic energy variation, so
$$ \Delta KE = W(\textrm{F}) = -\Delta E_p $$
Now, let's focus on Coulomb's law. Consider a charge $Q$ in $P$, and another one, say $q$ in $M$. It appears that $Q$ creates a force $\textbf{F}_e$ on $q$, with
$$\textbf{F}_e = \frac{Qq}{4\pi\epsilon \|\textbf{PM}\|^3}\textbf{PM}$$
Since $\textbf{F}_e$ only depends on the position of $P$ and $M$, it is derived from a potential energy $E_p$, and we find that
$$ E_p = \frac{Qq}{4\pi\epsilon\|\textbf{PM}\|}$$
However, even if the notation is usefull for a few particles, it really becomes difficul when studying a complex distribution of charges. Then, we prefer to see a single charge as the source of an electric field $\textbf{E}$, which exists everywhere, with $$\textbf{E}(M) = \frac{Q}{4\pi\epsilon \|\textbf{PM}\|^3}\textbf{PM}$$
Now, multiplying $\textbf{E}$ by the charge $q$ gives the electric force created by $Q$. So, $\textbf{E}$ and $\textbf{F}_e$ are almost the same: we introduce $U$ which is the potential created by $Q$, so that $E_p = qU$.