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.