This equation: $$\int_0^\infty \frac{x^{s-1}}{e^x-1} dx=\Pi(s-1)\zeta(s) $$ was used by Riemann in his famous paper from 1859. Seemingly it follows from: $$ \int_0^\infty e^{-nx}x^{s-1}dx= \frac{\Pi(s-1)}{n^s}$$
a result achieved by repeated integration by parts. How does the first equation follows from the second?
A faster way to find the second equation would be to simply make the substitution $u = nx$ which leads to :
$$ \int_0^\infty e^{-u} \left( \frac{u}{n}\right)^{s-1} \frac{du}{n} = \frac{\Pi(s-1)}{n^s} \quad \quad (*)$$
You can also write :
$$ \int_0^\infty \frac{x^{s-1}}{e^x - 1}dx = \int_0^\infty e^{-x} x^{s-1}\frac{1}{1 - e^{-x}}dx$$
And, you can write (geometric series) :
$$ \forall x >0, \ \frac{1}{1 - e^{-x}} = \sum_{n=0}^\infty e^{-nx}$$
which leads to
$$ \int_0^\infty \frac{x^{s-1}}{e^x - 1}dx = \int_0^\infty e^{-x} x^{s-1}\sum_{n=0}^\infty e^{-nx} dx $$
We can easily switch sum and integral, indeed :
$$ \forall n \in \mathbb N,\ \sum_{k=0}^n e^{-(n+1)x} x^{s-1} \leqslant e^{-x} x^{s-1} = \varphi(x) $$
so we can apply the dominated convergence theorem with $\varphi \in L^1(\mathbb R_+^*)$. Using $(*)$ :
$$ \int_0^\infty \frac{x^{s-1}}{e^x - 1}dx = \sum_{n=0}^\infty \frac{\Pi(s-1)}{(n+1)^s} $$
Which is exactly :
$$ \int_0^\infty \frac{x^{s-1}}{e^x - 1}dx = \Pi(s-1)\zeta(s) $$
First we have that $$\int_0^\infty \frac{x^{s-1}}{e^x-1}dx = \int_0^\infty x^{s-1}\sum _{n=0}^{\infty}e^{-nx}dx$$ Beliveing in uniform convergence we swap the sum and the integral to get $$\sum _{n=0}^{\infty}\int_0^\infty x^{s-1}e^{-nx}dx$$ with a simple change of variable, we set $nx = t$ so $dt = ndx$, then $$\sum _{n=0}^{\infty}\int_0^\infty x^{s-1}e^{-nx}dx = \sum _{n=0}^{\infty}\int_0^\infty \frac{1}{n}\left(\frac{t}{n}\right)^{s-1}e^{-t}dt = \underbrace{\sum_{n=0}^\infty n^{-s}}_{\text{Def. of }\zeta(s)}\underbrace{\int_0^\infty t^{s-1}e^{-t}dt}_{\text{Def. of }\Gamma(s)} = \zeta(s)\Gamma(s)$$ Just remember now the relation between the pi function and the Euler gamma function $$\Pi(z) = \Gamma(z+1)$$
