Finding Laurent expansion for $\zeta(s)$ is equivalent to finding a power series representation for
$$
F(s)=\zeta(s)-{1\over s-1}
$$
at $s=1$. This means that we need to develop strategies allowing us to deduce a formula for $F^{(k)}(s)$. This can be done by plugging the Dirichlet series representation of $\zeta(s)$ into Euler-Maclaurin formula:
$$
\zeta(s)=\sum_{n=1}^\infty{1\over n^s}={1\over s-1}+\frac12-s\int_1^\infty{\overline B_1(x)\over x^{s+1}}\mathrm dx
$$
Consequently we have for all $k>1$ that
$$
F^{(k)}(1)={\mathrm d^k\over\mathrm ds^k}\left[-s\int_1^\infty{\overline B_1(x)\over x^{s+1}}\mathrm dx\right]_{s=1}
$$
After simplifications, we can observe that
$$
\begin{aligned}
{\partial^k\over\partial s^k}[-sx^{-s-1}]
&=(-s)(-\log x)^kx^{-s-1}-k(-\log x)^{k-1}x^{-s-1} \\
&=(-1)^kx^{-s-1}[(-s)\log^kx+k\log^{k-1}x] \\
\end{aligned}
$$
As a result, we have
$$
\begin{aligned}
F^{(k)}(1)
&=(-1)^k\int_1^\infty{\overline B_1(x)[k\log^{k-1}x-\log^kx]\over x^2}\mathrm dx \\
&=(-1)^k\int_1^\infty\overline B_1(x)\mathrm d\left(\log^kx\over x\right)
\end{aligned}
$$
One can verify that this quantity converges
To eliminate $\overline B_1(x)$ in the above integral, we apply Euler-Maclaurin formula to $\log^kx/x$ (for $k>1$):
\begin{aligned}
\sum_{n=1}^N{\log^kn\over n}
&={\log^{k+1}N\over k+1}+{\log^kN\over2N}+\int_1^N\overline B_1(x)\mathrm d\left(\log^kx\over x\right) \\
&={\log^{k+1}N\over k+1}+\gamma_k+o(1)
\end{aligned}
where $\gamma_k$ is the Stieltjes constants:
$$
\gamma_k=\int_1^\infty\overline B_1(x)\mathrm d\left(\log^kx\over x\right)
$$
As a result, we can plug Stieltjes constants back into $F(s)$ to get
$$
F(s)=F(1)+\sum_{k=1}^\infty{(-1)^k\gamma_k\over k!}(s-1)^k
$$
Now it remains to determine $F(1)$, and using Euler-Maclaurin again on harmonic series allows us to determine $F(1)=\gamma$. Consequently, the Laurent expansion of $\zeta(s)$ at $s=1$ is as follows
$$
\zeta(s)={1\over s-1}+\sum_{k=0}^\infty{(-1)^k\gamma_k\over k!}(s-1)^k
$$
where
$$
\gamma_k=\lim_{N\to\infty}\left\{\sum_{n=1}^N{\log^kn\over n}-{\log^{k+1}N\over k+1}\right\}
$$
$\gamma_0=\gamma$ makes the above expression valid for $k\ge0$.