Proof of a few equations involving $\int_{\alpha}^{\infty}\frac{1}{t\left(e^{t}\pm1\right)}dt$ I derived these formulas with the Laurent series and Euler-Maclaurin summation formula. I can demonstrate this later if anyone's curious.  I'm wondering if there is another way.  I'm also interested in finding generalized formulas.
$$\lim\limits_{\alpha\to 0}\left[\frac{\ln\left(\alpha\right)}{2}+\int_{\alpha}^{\infty}\frac{1}{t\left(e^{t}+1\right)}dt\right]=\frac{1}{2}\left(\ln\left(\pi\right)-\ln\left(2\right)-\gamma\right)$$
$$\lim\limits_{\alpha\to 0}\left[\frac1\alpha+\frac{\ln\left(\alpha\right)}{2}-\int_{\alpha}^{\infty}\frac{1}{t\left(e^{t}-1\right)}dt\right]=\frac{1}{2}\left(\ln\left(\pi\right)+\ln\left(2\right)-\gamma\right)$$
A manipulation of these equations yields 
$$\lim\limits_{s\to -1}\left[\frac{1}{\ln|s|}+\left(-\frac{1}{s+1}+\frac{1}{2}\right)\ln|\ln|s||+\int_{s}^{\infty}\frac{\ln\left(\ln\left(u\right)\right)}{\left(u+1\right)^{2}}du\right]=\frac{1}{2}\left(\ln\left(\pi\right)-3\ln\left(2\right)-\gamma\right)$$
$$\lim\limits_{s\to 1}\left[\frac{1}{\ln\left(s\right)}+\left(\frac{1}{s-1}+\frac{1}{2}\right)\ln\left(\ln\left(s\right)\right)-\int_{s}^{\infty}\frac{\ln\left(\ln\left(u\right)\right)}{\left(u-1\right)^{2}}du\right]=\frac{1}{2}\left(\ln\left(\pi\right)+\ln\left(2\right)-\gamma\right)$$

Here is my elementary method of deriving these:
Begin with the Euler-Maclaurin summation formula:
$$
\begin{align}
\frac{1}{h}\int_a^b f(t)dt
&=\sum_{k=0}^n f(kh+a)-\left(\frac{f(a)+f(b)}{2}\right)
\\&-\sum_{k=1}^n \frac{h^{2k-1}B_{2k}}{(2k)!} \left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right)
\\&-R
\end{align}
$$
where $h=\frac{b-a}{n}$ and $R$ is the remainder term. Letting $n=\frac{b-a}{x}$ and rearranging we get
$$
\begin{align}
\sum_{k=1}^{(b-a)/x} \frac{x^{2k-1}B_{2k}}{(2k)!} \left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right)
&=\sum_{k=0}^{(b-a)/x} f(kx+a)-\left(\frac{f(a)+f(b)}{2}\right)
\\&-\frac{1}{x}\int_a^b f(t)dt
\\&-R
\end{align}
$$ 
Limiting $b\to 0$ and $a\to -\infty$, we have
$$
\begin{align}
\lim\limits_{\substack{%
     a \to -\infty\\
     b \to 0}}
\sum_{k=1}^{-a/x} \frac{x^{2k-1}B_{2k}}{(2k)!} \left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right)
&=\lim\limits_{\substack{%
     a \to -\infty\\
     b \to 0}}\left(\sum_{k=0}^{-a/x} f(kx+a)-\left(\frac{f(a)+f(b)}{2}\right)-\frac{1}{x}\int_a^b f(t)dt\right)
\end{align}
$$
The remainder disappears as $n\to\infty$. Now make a variable substitution in the limit $a\to -ax$.
$$
\begin{align}
\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\sum_{k=1}^{a} \frac{x^{2k-1}B_{2k}}{(2k)!} \left(f^{(2k-1)}(b)-f^{(2k-1)}(-ax)\right)
&=\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\left(\sum_{k=0}^{a} f((k-a)x)-\left(\frac{f(-ax)+f(b)}{2}\right)-\frac{1}{x}\int_{-ax}^b f(t)dt\right)
\\&=\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\left(\sum_{k=0}^{a} f(-kx)-\left(\frac{f(-ax)+f(b)}{2}\right)-\frac{1}{x}\int_{-ax}^b f(t)dt\right)
\end{align}
$$
Now use the following hint.
$$\frac{1}{z(e^z-1)}=\frac{1}{z^2}-\frac{1}{2z}+\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!}z^{2k-2}$$
which when we integrate we get
$$\begin{align}
\int_x^\infty\frac{1}{z(e^z-1)}dz
&=K-\left(-\frac{1}{x}-\frac{\ln{|x|}}{2}+\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!(2k-1)}x^{2k-1}\right)
\end{align}$$
Where $K$ stands for the integral evaluated at $\infty$. 
Let $f(x)=\text{Ei}\left(x\right)-\ln\left|x\right|-\gamma$. Note that $\lim\limits_{t\to 0}f^{(m)}(t)=\frac1m$ and $\lim\limits_{t\to -\infty}f^{(m)}(t)=0$ for $m\ge1$. Further note that $\lim\limits_{t\to 0} f(t) = 0$ so we may substitute the sum: $\sum_{k=0}^{a} f(-kx)=\sum_{k=1}^{a} f(-kx)$. Now we have
$$
\begin{align}
\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!(2k-1)}x^{2k-1}
&=\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\left(\sum_{k=1}^{a} f(-kx)-\left(\frac{f(-ax)+f(b)}{2}\right)-\frac{1}{x}\int_{-ax}^b f(t)dt\right)
\\&=\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\left(\sum_{k=1}^{a}\text{Ei}(-kx)-\sum_{k=1}^{a}\ln\left|-kx\right|-\sum_{k=1}^{a}\gamma-\left(\frac{f(-ax)+f(b)}{2}\right)-\frac{1}{x}\int_{-ax}^b f(t)dt\right)
\\&=-\int_x^\infty\frac{1}{z(e^z-1)}dz+\lim\limits_{\substack{%
     a \to \infty\\
     b \to 0}}\left(-a\ln|x|-\ln|a!|-a\gamma-\left(\frac{f(-ax)+f(b)}{2}\right)-\frac{1}{x}\int_{-ax}^b f(t)dt\right)
\\&=-\int_x^\infty\frac{1}{z(e^z-1)}dz+\frac{1}{x}+\frac{\ln\left|x\right|}{2}+\frac{1}{2}\left(\gamma-\ln\left(2\pi\right)\right)
\end{align}
$$
(The limit is tricky which is why I left out some steps). Therefore
$$\begin{align}
\int_x^\infty\frac{1}{z(e^z-1)}dz=\frac12 (\gamma-\ln(2\pi))+\frac{1}{x}+\frac{\ln{|x|}}{2}-\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!(2k-1)}x^{2k-1}
\end{align}$$
We can derive the other equation with $f(x)=-\text{Ei}(x)+2\text{Ei}(2x)-\ln|4x|-\gamma$.
 A: Too long for a comment: I join to @clathratus and would love to see derivation for the following reasons: I work on an explicit formula for the Laplace transform of the complex valued digamma function. The real one is discussed in  Dixit in the context of the famous integral OLIVIER OLOA
\begin{equation}
\int_0^{\frac{\pi }{2}} \frac{\theta^2}{\theta^2+\log (2\, \cos (\theta ))} \, d\theta = \frac{1}{8} \pi \, (1+\log (2 \,\pi )-\gamma )
\end{equation}
Since your formulas have some similarity your derivations could lead to a new ansatz for the proof.
EDIT
We start with the substitution of $y=t-\alpha$. Now we see, that:
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+\int_0^{\infty } \frac{1}{\left(e^{y}+1\right) (\alpha +y)} \, dy\right)=\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+\int_0^{\infty } \frac{1}{\left(e^{\alpha+y}+1\right) (\alpha +y)} \, dy\right)$$
The integral:
$$\mathcal{I}\left( \alpha \right) =\int_0^{\infty } \frac{1}{\left(e^{y}+1\right) (\alpha +y)} \, dy$$
is already discussed here Yuriy S. Jack D'Aurizio uses the "evaluating integrals over the positive real axis property of the Laplace transform" WIKIPEDIA to transform the integral to 
$$\mathcal{I}\left( \alpha \right) =\frac{1}{2} \int_0^{\infty } e^{-\alpha \, s} \left(\psi \left(\frac{s+2}{2}\right)-\psi \left(\frac{s+1}{2}\right)\right) \, ds$$
This can further written as the Laplace transform of 
$$\mathcal{I}\left( \alpha \right) =\mathcal{L\texttt{$\left\{\psi \left(\frac{s+2}{2}\right)\right\}$}}-\mathcal{L\texttt{$\left\{\psi \left(\frac{s+1}{2}\right)\right\}$}}$$
In Moll  the relationship between the Laplace transform of the diagamma function and the integral of  OLOA is shown. This ansatz could be used to proof the expression above.
EDIT-COMPLETE PROOF
Now, let's start with:
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+\frac{1}{2} \int_0^{\infty } e^{-\alpha\, s} \left(\psi \left(1+\frac{s}{2}\right)-\psi \left(\frac{1}{2}+\frac{s}{2}\right)\right)\, ds\right)$$
Transformation with $z\,=\,\frac{s}{2}$ leads to: 
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+\int_0^{\infty } e^{-2 \,\alpha \, z} \psi(1 + z) \, dz-\int_0^{\infty } e^{-2\, \alpha\, z} \psi\left(\frac{1}{2}+z\right) \, dz\right)$$
Now look at the first integral: An explicit formula for the Laplace transform of the digamma function is discussed in Dixit. For the readability, here the result:

Propositon:
  Let $\alpha > 0$,  $\gamma$  the Euler-Mascheroni constant and $\psi\left(x\right)$ the digamma function, then: 

\begin{align}
L_{C}(\alpha)= \int_0^{\infty } e^{-\alpha z} \psi (1+x) \, dx = 2 \alpha  \sum _{n=1}^{\infty } \frac{\log (n)}{\alpha ^2+4 \pi^2 n^2}+
\end{align}
\begin{equation*}
\left(\frac{1}{e^{\alpha }-1}-\frac{1}{\alpha }+1\right) \log \left(\frac{2 \pi}{\alpha }\right)+
\frac{1}{4} \left(\psi \left(\frac{i \alpha }{2 \pi }\right)+\psi \left(-\frac{i \alpha}{2 \pi }\right)\right)-\frac{\log (\alpha )+\gamma }{\alpha }
\end{equation*}
Inserting $L_{C}$ we get:
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+L_{C}(2\, \alpha)-\int_0^{\infty } e^{-2\, \alpha\,  z} \psi \left(\frac{1}{2}+z\right) \, dz\right)$$
The digamma function of the second integral can be transformed with The Wolfram Functions Site:
$$\psi (2\, z)=\frac{1}{2} \left(\psi (z)+\psi \left(z+\frac{1}{2}\right)\right)+\log (2)$$
Now we perform the integration:
$$2\, \log (2)\, \int_0^{\infty } e^{-2\, \alpha \, z} \, dz=\frac{\log (2)}{\alpha}$$ 
and use the well known relationship:
$$\psi(1 + z)\,=\,\psi(z)+\frac{1}{z} $$
to further simplify the expression:
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+L_{C}(2 \,\alpha)-\int_0^{\infty } e^{-2\, \alpha\,  z} (2 \psi (2\, z+1)-\psi (z+1)) \, dz\right)$$
Last but not least, we do again the transformation $y =2\,z$, to write the original expression as a function of $L_{C}(\alpha)$: 
$$\underset{\alpha \rightarrow 0}{\lim }\left(\frac{\log (\alpha )}{2}+\frac{\log (2)}{\alpha }+2 L_{C}(2 \,\alpha)-L_{C}(\alpha)\right)$$
Taking the limit $\alpha \rightarrow 0$ and considering that:
$$\underset{\alpha \rightarrow 0}{\lim }\left(2 \alpha  \sum _{n=1}^{\infty } \frac{\log (n)}{\alpha ^2+4\, \pi^2 \,n^2}\right)=0$$
leads to the expected resultat. This completes the proof.
A: For the first Limit, write ($s\rightarrow 1_-$, $a\rightarrow 0_+$ in this order)
$$
I(s,a)=J_1(s)-J_2(s,a)=\int_0^{\infty}dx\frac{x^{-s}}{(e^x+1)}-\frac12\int_0^{\alpha}x^{-s}(1+O(x))dx
$$
the first integral is a integral repesentation of the Dirichlet eta function (Proof: Taylor Expansion of the denonominator), which is related to the Riemann Zeta function as follows $\eta(s)=(2^s-1)\zeta(s)$, which is easily seen from the respective series representation. We therefore have
$$
I(s,a)=\Gamma(1-s)\eta(s)+\frac{a^{-s+1}}{2(s-1)}+O(a^{2-s}(s-2)^{-1})=\\ \color{red}{\Gamma(1-s)(2^s-1)\zeta(s)}+\color{blue}{\frac{a^{-s+1}}{2(s-1)
}}+\color{green}{O(a^{2-s}(s-2)^{-1})}\quad (\star)
$$
Taking the limit in $s$ we get (Proof:  expand both sides of the functional equation of the Zetafunction around $s=1$ and use $\Gamma(s)=1-\gamma(s-1)+O((s-1)^2)$ ),
$$
\color{red}{\Gamma(1-s)(2^s-1)\zeta(s)}=\color{red}{-\frac 1{2 (s-1)}+\frac{1}2 \left(-\gamma+\log(\pi/2)\right)+O(s-1)}
$$
and furthermore 
$$
\color{blue}{\frac{a^{-s+1}}{2(s-1)}}=\color{blue}{\frac 1{2 (s-1)}-\frac{1}2\log(a)+O(s-1)}
$$
as well as 
$$
\color{green}{O(a^{2-s}(s-2)^{-1})}=\color{green}{a+O(s-1)}
$$
the green part is clearly negligible as $a\rightarrow 0+$
so, as we subsitute back into $(\star)$ we get in total
$$
\lim_{a\rightarrow 0+} (I(1,a)+\color{blue}{\frac{1}2\log(a)})=\color{red}{\frac{1}2 \left(-\gamma+\log(\pi/2)\right)}
$$
as expected (note how the singular term exactly cancels, this was the motivation for the inital split of the integral).
The second Limit should be computable by the same method, taking one more term in the expansions of the different terms.
Edit: the singular part of the (exponentials) in the second Limit is $1/t+1/2$ which has to be substituted into $J_2(s,a)$
A: The result can also be found using a Mellin transform technique:
\begin{align}
I^+&=\int_{\alpha}^{\infty}\frac{1}{t\left(e^{t}+1\right)}\,dt\\
&=\frac{1}{2}\int_1^\infty \frac{e^{-\alpha u/2}}{\cosh (\alpha u/2)} \,du\\
&=\frac{1}{2}\int_0^\infty f(u)h\left( \frac{\alpha u}{2} \right)\,du
\end{align} 
where
\begin{align}
f(z)&=\begin{cases}
     0\text{ if } z\le 1\\
     z^{-1}\text{ if } z> 1
     \end{cases}\\
h(z)&=\frac{e^{-z}}{\cosh z}
\end{align}
Mellin transforms of these functions are (see Ederlyi TI 6.6.6, for example)
\begin{align}
\mathcal{M}\left[f(z)\right]&=\frac{1}{1-s} \text{ for }\Re s<1\\
\mathcal{M}\left[h(z)\right]&=2^{1-s}\left( 1-2^{1-s} \right)\Gamma(s)\zeta(s) \text{ for }\Re s>0
\end{align} 
with the Mellin convolution transform (DLMF),
we can choose $c>0$ to express
\begin{align}
I^+&=\frac{1}{4i\pi}\int_{c-i\infty}^{c+i\infty} 2^{1-s}\left( 1-2^{1-s} \right)\Gamma(s)\zeta(s)\left( \frac{\alpha}{2} \right)^{-s}\,\frac{ds}{s}
\end{align}
The oles of the function lie at $s =0,-1,-3,-5,\cdots$, with $s=0$ being double. Corresponding residues are $\ln\pi -2\ln2-\gamma-\ln\frac{\alpha}{2},\frac{\alpha}{2},-\frac{\alpha^3}{72},\cdots$. By closing the integral with the left half-circle, we obtain the expansion 
\begin{equation}
I^+\sim \frac{1}{2}\left( \ln\frac{\pi}{2\alpha} -\gamma\right)+\frac{\alpha}{4}-\frac{\alpha^3}{144}+O\left( \alpha^5 \right)
\end{equation} 
The same method applies for the second integral, with the tabulated transform
\begin{equation}
\mathcal{M}\left[\frac{e^{-z}}{\sinh z}\right]=2^{1-s}\Gamma(s)\zeta(s) \text{ for }\Re s>1
\end{equation} 
As $\zeta(-2n)=0$, the poles lie  at $s=1,0,-1,-3,\cdots$ ($s=0$ being double) with the corresponding residues $\frac{2}{\alpha},\ln\left( \frac{\alpha}{2\pi} \right)+\gamma,\frac{2}{\alpha},-\frac{\alpha}{6},\cdots$. This leads to the result
\begin{equation}
I^{-}\sim \frac{1}{2}\left( \ln\left( \frac{\alpha}{2\pi} \right)+\gamma \right)+\frac{1}{\alpha}-\frac{\alpha}{12}+\cdots
\end{equation} 
In both cases, the proposed expression in terms of the Bernoulli numbers can be retrieved as $\zeta(-2n-1)=-\frac{B_{2n+2}}{2n+2}$.
