Expressing Zeta function using Gamma series

My first thought is to give an integral representation for the general term of the series, $$\frac{1}{n^{1-s}}-\frac{B(n+s,1-s)}{\Gamma(1-s)}=\frac{1}{\Gamma(1-s)}\left(\int_{0}^{+\infty}x^{-s}e^{-n x}\,dx-\int_{0}^{1}x^{-s}(1-x)^{n+s-1}\,dx\right)$$

$$\frac{1}{n^{1-s}}-\frac{B(n+s,1-s)}{\Gamma(1-s)}=\frac{1}{\Gamma(1-s)}\left(\int_{0}^{+\infty}x^{-s}e^{-n x}\,dx-\int_{0}^{+\infty}(1-e^{-x})^{-s}e^{-(n+s)x}\,dx\right)$$ Summing over $n\geq 1$ we get:

$$\sum_{n\geq 1}\left(\frac{n^s}{n^1}-\frac{\Gamma(n+s)}{\Gamma(n+1)}\right)=\frac{1}{\Gamma(1-s)}\int_{0}^{+\infty}\left(\frac{1}{x^s}-\frac{1}{(e^x-1)^s}\right)\frac{dx}{e^x-1}$$ for every $s$ with real part $\in(0,1)$. The explicit computation of the last integral as $\,\frac{\pi}{\sin(\pi s)}+\Gamma(1-s)\,\zeta(1-s)\,$ proves OP's identity $(4)$. For the computation we may use, for instance, the classical application of Ramanujan's master theorem to Bernoulli polynomials.