regularization of a divergent integral
Edit
The formula proposed by Tom may be obtained with the classical method : $$\int_0^{\infty} \frac{x^{s-1}}{e^x(e^x-1)} dx $$ $$ = \sum_{n=2}^{\infty} \int_0^{\infty} x^{s-1} e^{-nx} dx$$ $$ = \sum_{n=2}^{\infty} \int_0^{\infty} \left(\frac{t}{n}\right)^{s-1} e^{-t} dt/n $$ $$ = \left(\sum_{n=2}^{\infty} \frac{1}{n^s}\right) \int_0^{\infty} t^{s-1} e^{-t} dt $$ $$ = \left(\zeta(s)-1\right) \Gamma(s)$$
As $s\to 0$ we have $\displaystyle \left(\zeta(s)-1\right) \Gamma(s)\sim -\frac 3{2s} +\frac 32\gamma-\frac{\ln(2\pi)}2$ (with $\gamma=\lim_{s\to 0} \frac 1s-\Gamma(s)$ the Euler constant $0.5772156649\cdots$ and $\frac{-\ln(2\pi)}2=\zeta'(0)$) which is still divergent...
I am not sure that we can avoid this because of the equivalence $\int \frac 1{x^2} dx$ near $0$ but perhaps that your 'regularization' could consist in the limit : $\displaystyle \lim_{\epsilon\to 0} \frac{\left(\zeta(-\epsilon)-1\right) \Gamma(-\epsilon)+\left(\zeta(\epsilon)-1\right) \Gamma(\epsilon)}2=\frac 32\gamma-\frac{\ln(2\pi)}2$
or a symmetrical pondered integral centered at $0$ with the same result...
Corresponding series
Let's substitute this definition of the Mittag-Leffler function : $$\operatorname{E}_s(z)=\sum_{k=0}^\infty \frac{z^k}{\Gamma(sk+1)}$$ in your integral : $$I(s)=s\int_{0}^{\infty}\frac{E_{s}(x^{s})-1}{xe^{x}(e^{x}-1)}dx$$ $$I(s)=s\int_{0}^{\infty}\frac{\sum_{k=1}^\infty \frac{x^{sk}}{\Gamma(sk+1)}}{xe^{x}(e^{x}-1)}dx=s\sum_{k=1}^\infty \frac 1{\Gamma(sk+1)}\int_{0}^{\infty}\frac{ x^{sk-1}}{e^{x}(e^{x}-1)}dx$$
We may use our previous result to rewrite $I(s)$ as : $$I(s)=s\sum_{k=1}^\infty \frac {\left(\zeta(sk)-1\right) \Gamma(sk)}{\Gamma(sk+1)}$$ this becomes the simple (and convergent if $\Re(s)>0$ and $\frac 1s$ not integer at least...) : $$\boxed{I(s)=\sum_{k=1}^\infty \frac {\zeta(sk)-1}{k}}$$
Approximation and Visualization
We may remove the simple poles in this series at $s=\frac 1n$ for $n\in \mathbb{N}^*$ by subtracting $\frac 1{sk-1}$ at the numerator (Stieltjes expansion) with the hope of getting something simpler.
Indeed for $s$ small I found this rather accurate numerical approximation (relative error $<\frac 1{20000}$ in the rectangle $\Re(s)\in (10^{-5},2),\ |\Im(s)|<1$) :
$$\sum_{k=1}^\infty \frac {\zeta(sk)-1-\frac 1{sk-1}}{k}\approx\frac {\ln(s)}2-0.95973512990747991661376-0.04053073339766363s+0.00026485233s^2$$
I think that the $\frac {\ln(s)}2$ term is exact (from numerical 'scaling properties') as well as the provided digits of the last three constants.
Since $\displaystyle \sum_{k=1}^\infty \frac 1{k(sk-1)}=-\gamma -\psi\left(1-\frac 1s\right)$ with $\psi$ the digamma function we have the approximation : $$I(s)\approx -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2-0.95973512990747991661376-0.04053073339766363s+0.00026485233s^2$$
that allows us to produce this picture of the imaginary part of $I(s)$ for $\Re(s)>0$ :
We observe a jump of amplitude $\pi$ (near the line $x=\Re(s)=0$ at the right around the line $y=\Im(s)=0$) coming from $\frac {\ln(s)}2$.
Let's observe that our expansion is analytic everywhere except at the simple poles $s=\frac 1n$ and the branch point $s=0$ allowing us to show the pictures of the imaginary and real part of the (possible) Riemann surface of the analytic extension of $I(s)$ on the whole plane (imagine an infinity of sheets superposed around the branch point $0$ generated by the $\frac {\ln(s)}2$ term) :
Asymptotic Expansion (direct method)
We want to rewrite following variant of $I(s)$ as a Maclaurin series $P(s)$ : $$J(s):=\sum_{k=1}^\infty \frac {\zeta(sk)-1-\frac 1{sk-1}}{k}-\dfrac{\ln(s)}2$$
Let's split the $-1$ at the numerator in two $-\frac 12$ terms to study the asymptotic behavior of this series (and obtain convergence of the individual 'power $\log$ integrals' we will see later...) : $$J(s)= \lim_{N\to\infty} J_N(s)\ \ \text{with}\ \ J_N(s)=\sum_{k=1}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
Let's first show that $\ \lim_{N\to\infty} J_N(s)$ exists for any $s>0$ such that $\frac 1s \notin \mathbb{N}$ (so that all the terms of the sum are defined).
Since $\ \sum_{k=1}^N \frac 1k= H_n=\gamma +\ln\bigl(N+\frac 12\bigr) + \operatorname{O}\bigl(\frac 1N\bigr)\ $ we have $\frac 12 \sum_{k=1}^N \frac 1k +\frac{\ln(s)}2=\frac {\gamma}2+\frac 12 \ln\bigl(s\bigl(N+\frac 12\bigr)\bigr)+\operatorname{O}\bigl(\frac 1N\bigr)=\frac {\gamma}2+\frac {\ln(sN)}2+\operatorname{O}\bigl(\frac 1N\bigr)$
For $sN \gg 1$ we have (setting $x:=sk$) :
$$s\sum_{k=N}^M \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}\sim \int_{sN}^{sM} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}{x}\,dx$$
the error is of order $\frac 1N$ because of the $H_M-H_N$ terms and because $\zeta(x)-1 \sim 2^{-x}$ for $x \gg 1$.
A precise evaluation of the last integral is : $\left[\frac 32 \ln(x) -\ln(x-1)+\sum_{k>1} \operatorname{Ei}(-x\ln(k))\right]^{sM}_{sN}$ (with $\operatorname{Ei}$ the exponential integral and $\operatorname{Ei}(-u)\sim \frac {e^{-u}}{-u}$ for $u\gg 1$) but this isn't needed since the integral is clearly equivalent to $\frac 12 \left[\ln(x)\right]^{sM}_{sN}\ $ with again an error of order $\frac 1N$ for $sN\gg 1$.
Both $\frac {\ln(sN)}2$ terms cancel in $J_N(s)$ for large $sN$ so that the limit $J(s)$ exists for $s>0$ (for $\frac 1s \in \mathbb{N}$ $J(s)$ may be defined as a limit since we removed the poles).
Let's use the Euler Maclaurin asymptotic formula with $f(k)=\frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}$ to expand directly $J_N(s)$ : $$\sum_{k=0}^N f(k) \sim \int_0^N f(x) dx +\frac {f(0)+f(N)}2+\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} \bigl(f^{(2k-1)}(N)-f^{(2k-1)}(0)\bigr)$$ becomes : $$ s\sum_{k=0}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}\sim \int_0^{sN} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}{x}\,dx + \frac {f(0)+f(N)}2+\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} \left(f^{(2k-1)}(N)-f^{(2k-1)}(0)\right)$$
Note that we started at $k=0$ with $\displaystyle\ f(0)=\lim_{k\to 0}\ s\frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{sk}=(1+\zeta'(0))s\ $
and that, for any $n\in \mathbb{N}^*$, we have :
$f^{(n-1)}(0)=\lim_{N\to 0}\ s\left(\frac d{dN}\right)^{n-1}\ \left(\frac {\zeta(sN)-\frac 12-\frac 1{sN-1}}{sN}\right)=\lim_{N\to 0}\ s\left(\frac {s\cdot d}{d(sN)}\right)^{n-1} (\cdots)=\lim_{x\to 0}\ s^n\left(\frac {\zeta(x)-\frac 12-\frac 1{x-1}}x\right)^{(n-1)}=s^n\lim_{x\to 0} \frac{\left(\zeta(x)-\frac 1{x-1}\right)^{(n)}}n$
or simply $\ \displaystyle f^{(n-1)}(0)=\frac{n!+\zeta^{(n)}(0)}n s^n\ $ for $n>0$
For $n\ge 0$ and $N \gg 1$ we may show that $f^{(n)}(N)\sim \frac {(-1)^n n!}{2 N^{n+1}}\ $ so that the Euler Maclaurin asymptotic formula becomes (after removing $\frac{H_N+\ln(s)}2$ on both sides to get $J_N(s)$ at the left) :
$$J_N(s)= s\sum_{k=0}^N \frac {\zeta(sk)-1-\frac 1{sk-1}}{sk}-\ln(s)\sim \int_0^{sN} \frac {\zeta(x)-\frac 12-\frac 1{x-1}}x\,du - \frac {\gamma}2-\frac {\ln(sN)}2+\operatorname{O}\bigl(\frac 1N\bigr)+ \frac {f(0)}2-\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} f^{(2k-1)}(0)$$
At this point let's rewrite $\ln(sN)=\int_0^{sN} \frac 1{x+1}\,du+\operatorname{O}\bigl(\frac 1N\bigr)$ and go at the limit $N\to \infty$ : $$J(s)\sim - \frac {\gamma}2+\int_0^\infty \frac {\zeta(x)-\frac 12-\frac 1{x-1}}x-\frac 1{2(x+1)}\,dx + \frac {f(0)}2-\sum_{k=1}^\infty \frac{B_{2k}}{(2k)!} f^{(2k-1)}(0)$$
We recognize a Maclaurin series in $s$ : $$J(s)\sim K_0 -K_1 \frac s2-\sum_{k=1}^\infty K_{2k} \frac{B_{2k}}{(2k)!} s^{2k} $$
with the constant term given by $\displaystyle\ \boxed{K_0:= -\frac{\gamma}2 +\int_0^\infty \frac {\zeta(x)-\frac 1{x-1} -\frac 12 }x -\frac 1{2(1+x)}\,dx}$
(the numerical value is $K_0=-0.95973512990\cdots$ as wished ; I don't know a closed form...)
and the $K_n$ coefficients of $s^n$ by $\boxed{K_n:=\frac{n!+\zeta^{(n)}(0)}n}\ \text{for} \ n>0$
The wished analytic extension of $I(s)$ will be :
$$\boxed{I(s)\sim -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2+K_0 -K_1 \frac s2-\sum_{k=1}^\infty K_{2k} \frac{B_{2k}}{(2k)!} s^{2k}}$$
Short table of $K_n$ values :
($\gamma_n$ is the nth Stieltjes constant : $\gamma_0= \gamma,\ \gamma_1= 0.07281584548367672486058\cdots$)
$\displaystyle K_1=1+\zeta'(0)=1-\frac{\ln(2\pi)}2=0.0810614667953272582\cdots$
$\displaystyle K_2=1+\frac{\zeta''(0)}2=1-\frac{\ln(2\pi)^2}4+\frac{\gamma^2}4-\frac{\pi^2}{48}-\frac{\gamma_1}2=
-0.0031782279542924256\cdots$
$K_3=\frac{3!+\zeta'''(0)}3=-0.00157038895408481592\cdots$
$K_4=\frac{4!+\zeta^{(4)}(0)}4=0.0007242029965730102531\cdots$
($-\frac {K_1}2$ and $-\frac {K_2}{12}$ may be compared to the coefficients for $s$ and $s^2$ obtained previously : $-0.040530733397\cdots$ and $0.00026485233\cdots$).
(it is interesting to observe that $\zeta^{(n)}(0)\approx -n!$ as well as $\zeta^{(n)}\left(\frac 12\right)\approx -2^{n+1}n!$)
Asymptotic Expansions (the same result after an excursion with Mittag-Leffler, another expression for $K_0$)
To get the exact terms of the expansion of $I(s)$ let's provide the Maclaurin expansion in $s$ of $s\left(E_s(x^s)-1\right)$ (the coefficients are functions of $x$ that we will have to integrate later over $\mathbb{R}^+$) :
$$s\left(E_s(x^s)-1\right)=C_0(x) e^x-\frac s2-\frac{s^2}{12}\left[\ln(x)+\gamma\right]+\frac{s^4}{720}\left[\ln(x)^3+3\gamma\ln(x)^2+3\left(\gamma^2-\zeta(2)\right)\ln(x)+\gamma^3-3\gamma\zeta(2)+2\zeta(3)\right]-\frac{s^6}{30240}\left[\ln(x)^5+5\gamma\ln(x)^4+10\left(\gamma^2-\zeta(2)\right)\ln(x)^3+10\left(\gamma^3-3\gamma \zeta(2)+2\zeta(3)\right)\ln(x)^2+5\left(\gamma^4-6\gamma^2\zeta(2)+8\gamma \zeta(3)+\pi^4/60\right)\ln(x)+\gamma^5-10\gamma^3\zeta(2)+20\gamma^2\zeta(3)+\gamma\pi^4/12-20\zeta(2)\zeta(3)+24\zeta(5) \right]+P_8(\ln(x))\operatorname{O}\bigl(s^8\bigr)$$
The term independent of $s$ : $C_0(x) e^x$ with $\displaystyle C_0(x)=\int_0^1 \frac{\gamma(t,x)}{\Gamma(t)}\,dt\ $ (and $\gamma(t,x)$ the lower 'incomplete gamma function') was obtained by considering the limit as $n\to \infty$ of formula (3.4) from Haubold's paper (with $s=\frac 1n$) : $$E_{\frac 1n}\left(x^{\frac 1n}\right)=e^x\left[1+\sum_{r=1}^{n-1}\frac{\gamma\left(1-\frac rn ,x\right)}{\Gamma\left(1-\frac rn\right)}\right]$$
A shorter expression for the Maclaurin expansion (that could again be obtained with Euler-Maclaurin) should be :
$$\boxed{s\left(E_s(x^s)-1\right)\sim C_0(x) e^x-\frac s2-\sum_{n=2}^\infty s^n\frac{B_n}{n!}\ \sum_{k=0}^{n-1} \binom{n-1}{k} S_k \ln(x)^{n-1-k}} $$
with $B_n$ the 'Bernoulli numbers' ($B_{2n+1}=0$ for $n>0$) and $S_k$ the coefficients appearing in the expansion of the reciprocal gamma function at $z=0$ (and at $z=1$) : $$\frac 1{\Gamma(s)}=\sum_{k=0}^\infty \frac {S_k}{k!} s^{k+1}\ \ \text{so that}\ \ S_k=\lim_{s\to 0} \left(\frac 1{\Gamma(s+1)}\right)^{(k)}$$
$S_0=1$
$S_1=\gamma$
$S_2=\gamma^2-\zeta(2)$
$S_3=\gamma^3-3\gamma \zeta(2)+2\zeta(3)$
$S_4=\gamma^4-6\gamma^2\zeta(2)+8\gamma \zeta(3)+\pi^4/60$
$S_5=\gamma^5-10\gamma^3\zeta(2)+20\gamma^2\zeta(3)+\gamma\pi^4/12-20\zeta(2)\zeta(3)+24\zeta(5)$
$C_0$ too may be rewritten with the $S_k$ terms as $\displaystyle C_0(x)=\sum_{k=0}^\infty \frac{S_k}{(-\ln(x))^{k+1}}$.
We will again start with : $\displaystyle J(s)= \lim_{N\to\infty} J_N(s)$ with $\ \displaystyle J_N(s)=\sum_{k=1}^N \frac {\zeta(sk)-\frac 12-\frac 1{sk-1}}{k}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$
but will use Mittag-Leffler after that : $$J_N(s)=s\sum_{k=1}^N \frac {\zeta(sk)\Gamma(sk)-\frac {\Gamma(sk)}2-\frac {(sk-1)\Gamma(sk-1)}{sk-1}}{(sk)\Gamma(sk)}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
$$J_N(s)=s\sum_{k=1}^N \frac {\zeta(sk)\Gamma(sk)-\frac {\Gamma(sk)}2-\Gamma(sk-1)}{\Gamma(sk+1)}-\frac 12 \sum_{k=1}^N \frac 1k -\frac{\ln(s)}2$$
Let's rewrite everything in integral form (using $\ \sum_{k=1}^N \frac 1k= \gamma +\int_0^N \frac 1{1+x} + \operatorname{O}\bigl(\frac 1N\bigr)\ $, $\Gamma(u)=\int_0^\infty t^{u-1} e^{-t} dt\ $ and considering the limit as $N\to \infty$) :
$$J_N(s)=\sum_{k=1}^N \left( \frac s{\Gamma(sk+1)}\int_{0}^{\infty}\frac{x^{sk-1}}{e^x-1}-\frac {x^{sk-1}}{2e^x}-\frac {x^{sk-1}}{x e^x}\,dx\right)-\frac{\gamma}2 -\frac{\ln(N s)}2+\operatorname{O}\bigl(\frac 1N\bigr)$$
$$J_N(s)=-\frac{\gamma}2+\int_{0}^{\infty}\frac{\sum_{k=1}^N \frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)\,dx -\frac{\ln(N s)}2+\operatorname{O}\bigl(\frac 1N\bigr)$$
$$J_N(s)=-\frac{\gamma}2 +\int_{0}^{\infty}\frac{\sum_{k=1}^N \frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)+\frac{e^{-N s x}-e^{-x}}{2x}\,dx +\operatorname{O}\bigl(\frac 1N\bigr)$$
that becomes at the limit $N \to \infty$ :
$$J(s)\sim -\frac{\gamma}2 +\int_{0}^{\infty}\frac{\sum_{k=1}^\infty\frac {sx^{sk}}{\Gamma(sk+1)}}{xe^{x}}\left(\frac {e^x}{e^{x}-1}-\frac 12-\frac 1x\right)-\frac 1{2(1+x)}\,dx$$
$$\boxed{J(s)\sim -\frac{\gamma}2 +\int_{0}^{\infty}\frac{s \left(E_{s}(x^{s})-1\right)}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)-\frac 1{2(1+x)}\,dx}$$
From the MacLaurin expansion of $s \left(E_{s}(x^{s})-1\right)$ in powers of $s$, we get :
$$J(s)\sim K_0 -K_1 \frac s2-\sum_{n=2}^\infty K_n \frac{B_n}{n!} s^n $$
and the wished analytic extension of $I(s)$ :
$$\boxed{I(s)\sim -\gamma -\psi\left(1-\frac 1s\right) +\frac {\ln(s)}2+K_0 -K_1 \frac s2-\sum_{n=2}^\infty K_n \frac{B_n}{n!} s^n}$$
with the constant term given by : $$K_0:= -\frac{\gamma}2 +\int_0^\infty \frac{C_0(x)}x\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)-\frac 1{2(1+x)}\,dx$$
with the same numerical value and the $K_n$ coefficients of $s^n$ for $n>0$ by :
$$K_n:= \sum_{m=0}^{n-1} f_m \binom{n-1}{m} S_{n-1-m}$$
with :
$$f_m:=\int_{0}^{\infty}\frac{\ln(x)^m}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)\,dx$$
Let's use (with a derivation similar to the previous one!) :
$$Z(s):=\left(\zeta(s)-\frac 1{s-1}-\frac 12\right) \Gamma(s)=\int_0^{\infty} \frac{x^s}{xe^{x}}\left(\frac 1{e^{x}-1}-\frac 1x+\frac 12\right)\,dx$$ (the idea of the additional term $-\frac 1x+\frac 12$ appeared in the Hermite-Stieltjes letters of 1885 allowing to extend the integral giving zeta near $x=0$ by removing the $\frac 1{x^2}$ and $\frac 1x$ singularities ; just some pages later we find the famous claim of Stieltjes of a proof of R.H. !)
since $x^s=e^{s\ln(x)}$ we get : $$f_m= Z^{(m)}(0)=\lim_{s\to 0}\left(\frac {d}{ds}\right)^m \left[\left(\zeta(s)-\frac 1{s-1}-\frac 12\right)\Gamma(s)\right]$$
From the binomial theorem it is clear that : $$K_{n+1}= \sum_{m=0}^n \binom{n}{m} f_m S_{n-m}= \lim_{s\to 0}\sum_{m=0}^n \binom{n}{m} Z^{(m)}\left(\frac 1{\Gamma(s+1)}\right)^{(n-m)}= \lim_{s\to 0} \left(\frac {Z(s)}{\Gamma(s+1)}\right)^{(n)}$$
allowing some simplifications : $$K_n=\lim_{s\to 0} \Bigl(\frac{\zeta(s)-\frac 1{s-1}-\frac 12}s\Bigr)^{(n-1)}=\lim_{s\to 0} \frac{\left(\zeta(s)-\frac 1{s-1}\right)^{(n)}}n\ \text{for} \ n>0$$ $$\text{i.e.}\ \ \boxed{K_n=\frac{n!+\zeta^{(n)}(0)}n}\ \text{for} \ n>0$$ i.e. the same result as previously
We could continue this very interesting subject and study the multiplicative partition (“Factorisatio Numerorum”) obtained by rewriting $\displaystyle I(s)=\sum_{k=1}^\infty \frac {\zeta(sk)-1}{k}$ as the logarithm of $\displaystyle \prod_{n>1} \frac 1{1-n^{-s}}$ getting this Dirichlet series or... but I'll have to make all this visible one day so...
Sorry if all this is not exactly what you hoped...