Closed-form of $\int_0^1 B_n(x)\psi(x+1)\,dx$ Is there a closed-form of the following integral?
$$I_n = \int_0^1 B_n(x)\psi(x+1)\,dx,$$
where $B_n(x)$ are the Bernoulli polynomials and $\psi(x)$ is the digamma function.
The motivation of the problem was this question. In this answer achille hui showed that for all $n>0$
$$\int_0^1 B_n\left(\left\{ \frac1x \right\}\right) \frac{dx}{x}
= -\int_0^1 B_n(u) \psi(1+u) du,$$
where $\{x\}$ is the fractional part of $x$.
By looking at specific $n$ values we could observe interesting patterns.
$$\begin{align}
I_0 & =  0\\
I_1 & =  1-\frac 12 \ln(2\pi) \\
I_2 & =  -\frac 12 + 2 \ln A \\
I_3 & = \frac{1}{12} - \frac{3}{4\pi^2}\zeta(3) = \frac{1}{12} - \frac{1}{8}\frac{\zeta(3)}{\zeta(2)}\\
I_4 & = \frac{1}{45} - 4\zeta'(-3)\\
I_5 & = -\frac{13}{360} + \frac{15}{4\pi^4}\zeta(5) = -\frac{13}{360} + \frac{1}{24}\frac{\zeta(5)}{\zeta(4)}\\
I_6 & = -\frac{1}{252}-6\zeta'(-5)\\
I_7 & = \frac{47}{1260} - \frac{315}{8\pi^6}\zeta(7) = \frac{47}{1260} - \frac{1}{24}\frac{\zeta(7)}{\zeta(6)}\\
I_8 & = -\frac{8}{1575}-8\zeta'(-7)\\
I_9 & = -\frac{1703}{25200} -\frac{2835}{4\pi^8}\zeta(9) = -\frac{1703}{25200} -\frac{3}{40}\frac{\zeta(9)}{\zeta(8)}\\
I_{10} & = \frac{2461}{83160} - 10\zeta'(-9)\\
\dots
\end{align}$$
where $A$ is the Glaisher–Kinkelin constant, $\zeta(3)$ is the Apéry's constant and in general $\zeta$ is the Riemann zeta function and $\zeta'$ is its derivative.
 A: This is an answer for the odd case.
Proposition. Let $k=1,2,3,\ldots$. Then

$$
\int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx=(-1)^{k+1}\frac{(2k+1)!}{(2\pi)^{2k+1}}\pi \: \zeta(2k+1)-\sum_{j=0}^{2k}\!\frac{ {{2k+1}\choose j} B_j}{2k+1-j} \quad (*)
$$

The proof is here, observing that
$$\begin{align}
\int_0^1 B_{2k+1}(x)\: \psi (x+1) \:dx & =  \int_0^{1} \frac{B_{2k+1}\left(\left\{ 1/x \right\}\right)}{x}dx. \tag1
\end{align}$$ A proof of $(1)$ may be found here.
A: The complete answer for the question with $\,n\ge 1\,$ is the union of the answers of Start wearing purple and Olivier Oloa (the sign of the right side of the equation $(*)$ is wrong there), we get:  

$$\int\limits_0^1 B_n(x)\psi(x+1)dx = (-1)^{n-1} \left(n \,\zeta’(1-n)- B_n H_{n-1}\right) + \sum\limits_{k=1}^n {\binom n k} \frac{B_{n-k}}{k}$$  

A well-functioning trick is to invert $(1)$ of Start wearing purple’s answer by using: 
$$b_{m-1}=\sum\limits_{v=0}^{m-1} {\binom m v} a_v \enspace \Leftrightarrow \enspace a_m=\frac{1}{m+1}\sum\limits_{v=0}^m {\binom {m+1} {v+1}} B_{m-v}b_v $$  
Note 1:
$$\sum\limits_{k=1}^n {\binom n k} \frac{B_{n-k}}{k}=\int\limits_0^1\frac{B_n(x)-B_n(0)}{x}dx$$
Note 2:
$$n\zeta’(1-n)- B_n H_{n-1}=\frac{d}{dt}B_t(1)|_{t=n}-H_n B_n(1)$$
$\hspace{4cm}$ Informations about $\,B_t(x)\,$ can be found here .
A: The paper

V. Adamchik, Polygamma functions of negative order, J. Comp. Appl. Math. 100, (1998), 191-199) 

contains the following statement (Proposition 3, Eq. (17)):

Let $n\in\mathbb{Z}_{\ge0}$ and $\Re z>0$, then
  \begin{align}\int_0^z x^n\psi\left(x\right) dx=&(-1)^{n-1}\zeta'(-n)+\frac{(-1)^n}{n+1}B_{n+1}H_n+\\ &
\sum_{k=0}^n(-1)^kz^{n-k}{n \choose k}\left[\zeta'(-k,z)-\frac{B_{k+1}(z)H_k}{k+1}\right],\end{align}
  where $B_n$ and $B_n(z)$ are Bernoulli numbers and polynomials, $H_n$ are harmonic numbers, and $\zeta(s,z)$ denotes Hurwitz zeta function.

Specializing to $z=1$, we get
$$\int_0^1 x^n\psi\left(x\right) dx=\sum_{k=0}^{n-1}(-1)^k{n \choose k}\left[\zeta'(-k)-\frac{B_{k+1}H_k}{k+1}\right].\tag{1}$$
Using (1) and the recursion relation $\psi(x+1)-\psi(x)=x^{-1}$, the integral $\int_0^1 P\left(x\right)\psi\left(x+1\right) dx$ can be computed for any polynomial $P(x)$.
Edit: There remains the task of explaining why Bernoulli polynomials provide a better basis than monomials $x^n$ (instead of linear combination of $\zeta'(-k)$ we get only one of these values). I think this could be understood after careful reading of the paper.
