This is a question that I wish more undergrads would ask! It’s a really important calculation, because it really requires you to “go back to the fundamentals” when you answer it.
First, let me remove an unnecessary complication you are adding: the Lorentz boost is linear and so if we solve this problem for one particle we will solve it by superposition for $N$ particles. You only need to handle one at a time. Furthermore we know that we can choose our coordinates arbitrarily so let us set our one particle moving in the xy-plane through the origin and then boost in the x-direction, that makes this essentially only a 2-dimensional problem. (If you do not like that multiply through everything by a $\delta(z).$)
Let me do as is my custom and define the relativistic time coordinate $w = ct$ as well, with the charged particle's velocity vector $\mathbf u$ transformed into some vector $\vec\eta = (\eta_x, \eta_y) = \mathbf u/c.$ In this language we would write $$\rho(w, x, y) = q~\delta(x - \eta_x w)\delta(y - \eta_y w),\\\mathbf J(w, x, y) = \vec\eta~c\rho(w, x, y).$$ We have our 2D setup in the unprimed coordinates, great.
Now comes the part where we “go back to the fundamentals” and ask, what does the transformed charge density $\rho'(w', x', y')$ really mean? We have to go back to the definition, it is the charge contained in the box formed by the product of open intervals $(x', x' + \delta x') \times (y', y' + \delta y')$ at the time $w'$, divided by $\delta x'~\delta y'$. We can then transform that primed-stationary box back to the unprimed coordinates and we see that it is moving across space, and it is length-contracted. The length contraction turns out to be kind of gorgeously irrelevant here due to the fact that we're dealing with a point charge, as is the fact that “simultaneous” in the primed frame is not “simultaneous” across this square in the unprimed frame: we anticipate that the charge “should be” a Lorentz-scalar so we just need to find the time when the box hits our point charge.
So the box is a set of points $(x, y) \in (\gamma~(x' + \beta w'), \gamma~(x' + \beta w') + \delta x'/\gamma) \times (y', y' + \delta y')$ at time $w = \gamma~(w' + \beta x')$. And so the point particle described by $\rho$ is inside this volume only when $$x' + \beta w' = \eta_x (w' + \beta x')\\
y' = \eta_y ~\gamma~(w' + \beta x').
$$
This in turn implies that the point charge is seen traveling in the primed coordinates at the point $$x' =
\frac{\eta_x -\beta}{1 - \eta_x \beta}~w',\\
y' = \frac{\gamma^{-1}~\eta_y} {1 - \eta_x \beta}~ w',$$ and so if it is to still be a point charge with magnitude $q$ (that is, if we assume that charge is a Lorentz scalar) then the proper expression for $\rho', \mathbf J'$ must be just, $$\rho'(w', x', y') = q~\delta\left(x' - \frac{\eta_x -\beta}{1 - \eta_x \beta}~w'\right) ~\delta\left(y' - \frac{\gamma^{-1}~\eta_y} {1 - \eta_x \beta}~ w'\right),\\J_x' = \frac{\eta_x - \beta}{1 - \eta_x \beta} ~c\rho',\\
J_y' = \frac{\gamma^{-1}~\eta_y} {1 - \eta_x \beta}~ c\rho'.$$ Now compare this with the naive transformation of the fields we have, $$
\begin{align}
\rho(w', x', y') &= q~\delta\big(\gamma (x' + \beta w') - \eta_x \gamma (w' + \beta x')\big) \delta\big(y' - \eta_y \gamma (w' + \beta x')\big)\\
&=\frac{q}{\gamma~(1-\beta\eta_x)} \delta\left(x' - \frac{\eta_x - \beta}{1 - \beta\eta_x}\right)\delta\left(y' - \eta_y \frac{\gamma^{-1}~\eta_y} {1 - \eta_x \beta}~ w'\right)\\
&= \frac{\rho'}{\gamma~(1 - \beta \eta_x)},
\end{align}$$and note how this special pre-factor comes in because $\delta(a x) = \delta(x)/|a|.$ Once you see this pre-factor, I claim you have solved the entire problem. At long last you find that in fact, $$ c \rho' = \gamma ~(c \rho - \beta J_x),\\
J_x' = \gamma~ (J_x - \beta c \rho),\\
J_y' = J_y,$$ and therefore the single particle current density has transformed like a 4-vector under boosts, with transformation matrix $$\begin{bmatrix}\gamma&-\gamma\beta&0\\-\gamma\beta&\gamma&0\\0&0&1\end{bmatrix}.$$ Combine this with the fact that $\mathbf J$ transforms like an ordinary vector under rotations, and you have that arbitrary single-particle $(c\rho, \mathbf J)$s are 4-vectors, and thus their superpositions are, too.