Prove $\int^\infty_0 \frac{e^{-t}}{t}\left[\frac1{t^2}-\frac1{(1-e^{-t})^2}+\frac1{1-e^{-t}}-\frac1{12}\right]dt=\frac34-\zeta'(-1)+\zeta'(0)$ How to prove
$$\int^\infty_0 \frac{e^{-t}}{t}\left[\frac1{t^2}-\frac1{(1-e^{-t})^2}+\frac1{1-e^{-t}}-\frac1{12}\right]dt=\frac34-\zeta'(-1)+\zeta'(0)$$
?
This integral appeared in my answer, and according to an arXiv paper and the OP's conjecture this equality is very likely to be true. This is also supported by numerical evidence.
You cannot find the integral in the arXiv paper, as the integral arises in my lengthy proof of a statement (whose proof is omitted) in the paper. Thus I think it is not very useful to provide the link.
Real or complex approaches are welcomed. Thanks in advance.
 A: I will first evaluate $$I(s) = \int_{0}^{\infty} t^{s-1}e^{-t} \left(\frac{1}{t^{2}} - \frac{1}{(1-e^{-t})^{2}} + \frac{1}{1-e^{-t}} - \frac{1}{12} \right) \, \mathrm dt \tag{1}$$ for $\operatorname{Re}(s)>2$.  
Then to find $I(0)$, I'll take the limit as $s$ approaches $0$.
(The Mellin transform, like the Laplace transform, is an analytic function where the integral converges absolutely. Since the integral on the right side of $(1)$ behaves like $t^{s+1}$ near $t=0$, the integral defines an analytic function for $\operatorname{Re}(s) >-2$.)
With the restriction that $\text{Re}(s)>2$, we can break up the integral into four separate convergent integrals and evaluate each integral separately.
The first integral is just $\Gamma(s-2)$, the third integral is $\Gamma(s) \zeta(s)$, and the fourth integral is $\frac{1}{12} \, \Gamma(s)$. 
The second integral is $\Gamma(s) \zeta(s-1)$, which can be derived by differentiating the integral of the Bose–Einstein distribution and using the polylogarithm property $\frac{\mathrm d}{\mathrm dz} \operatorname{Li}_{s}(z) = \frac{\operatorname{Li}_{s-1}(z)}{z}$ along with the fact that $\operatorname{Li}_{s}(1) = \zeta(s)$.
Combing all 4 integrals, we get   $$I(s) = \left(\Gamma(s-2) -\Gamma(s) \zeta(s-1) + \Gamma(s)\zeta(s) - \frac{1}{12} \,\Gamma(s)\right) $$ for $\operatorname{Re}(s) >2$.
The above expression is an analytic function for $\operatorname{Re}(s) >-2$. (The singularities at $s=2$, $s=1$, $s=0$, and $s=-1$ are removable.) Combining this with the property of the Mellin transform mentioned previously, it follows by the identity theorem that $$I(s) =\left(\Gamma(s-2) -\Gamma(s) \zeta(s-1) + \Gamma(s)\zeta(s) - \frac{1}{12} \,\Gamma(s)\right) $$ for $\operatorname{Re}(s) >-2$.  
To find $I(0)$, we need to expand the terms in Laurent series at $s=0$.
The gamma function has simple poles at zero and the negative integers with residue $\frac{(-1)^{n}}{n!}$.
See here.
So at $s=-2$, $\Gamma(s) = \frac{1}{2(s+2)} + \mathcal{O}(1)$.
The constant term of Laurent series of $\Gamma(s)$ at $s=-2$ is then $$ \begin{align} \lim_{s \to -2} \left(\Gamma(s) - \frac{1}{2(s+2)} \right) &= \lim_{s \to -2} \left(\frac{\Gamma(s+3)}{s(s+1)(s+2)} - \frac{1}{2(s+2)}\right) \\ &= \lim_{s \to -2} \frac{2 \Gamma(s+3)-s(s+1)}{2s(s+1)(s+2)} \\ &= \frac{1}{2}\lim_{ s \to -2} \frac{2 \Gamma'(s+3)-2s-1}{(s+1)(s+2)+s(s+2)+s(s+1)}\\  &= \frac{2\Gamma'(1)+3}{4} \end{align}$$
Since the Laurent series of $\Gamma(s)$ at $s=-2$ has the same coefficients as the Laurent series of $\Gamma(s-2)$ at $s=0$, we get $$ \Gamma(s-2) = \frac{1}{2s} +\frac{\Gamma'(1)}{2} + \frac{3}{4} + \mathcal{O}(s)$$
Similarly, $$\Gamma(s) =  \frac{1}{s} + \Gamma'(1) + \mathcal{O}(s)$$
Therefore, $$\small I(0) = \lim_{s \to 0} \left[\frac{1}{2s} + \frac{\Gamma'(1)}{2} + \frac{3}{4} + \mathcal{O}(s) - \left(\frac{1}{s} + \Gamma'(1) + \mathcal{O}(s) \right) \left(\zeta(-1) - \zeta(0) + \zeta'(-1)s - \zeta'(0)s+ \mathcal{O}(s^{2}) \right) - \frac{1}{12} \left(\frac{1}{s} + \Gamma'(1) + \mathcal{O}(s) \right)\right] $$
where $\zeta(0) = - \frac{1}{2}$ and $\zeta(-1) = -\frac{1}{12}$.
This leads to mass cancellation, and we end up with $$I(0) = \lim_{ s \to 0} \left( \frac{3}{4} - \zeta'(1) + \zeta'(0) + \mathcal{O}(s)\right) = \frac{3}{4} - \zeta'(1)  +\zeta'(0) $$
