Help with logarithmic definite integral: $\int_0^1\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\, dx$ I'm look for a closed form evaluation of the following improper definite integral involving logarithms:

$$\begin{align}
I:&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\
&\approx 0.579307.
\end{align}$$


My attempt:
My first idea was to transform the integral using the substitution $u=-\ln{(x)}$, but this didn't really result in anything recognizably easier than the original integral:
$$\begin{align}
I&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\
&=\int_{0}^{\infty}(-u)\ln^3{(1-e^{-u})}\,\mathrm{d}u\\
&=-\int_{0}^{\infty}u\ln^3{(1-e^{-u})}\,\mathrm{d}u.
\end{align}$$
The second thing I tried was a substitution I've known to be successful with similar problems: 
$$\begin{align}
u&=-\ln{(1-x)},\\
-u&=\ln{(1-x)},\\
e^{-u}&=1-x,\\
x&=1-e^{-u},\\
\mathrm{d}x&=e^{-u}\,\mathrm{d}u.
\end{align}$$
Applying the substitution, the integral becomes:
$$\begin{align}
I&=\int_{0}^{1}\frac{1}{x}\ln{(x)}\ln^3{(1-x)}\,\mathrm{d}x\\
&=\int_{0}^{\infty}\frac{(-u)^3}{1-e^{-u}}\ln{(1-e^{-u})}\,e^{-u}\,\mathrm{d}u\\
&=-\int_{0}^{\infty}\frac{u^3\,e^{-u}}{1-e^{-u}}\ln{(1-e^{-u})}\,\mathrm{d}u\\
&=\int_{0}^{\infty}\frac{u^3}{1-e^{u}}\ln{(1-e^{-u})}\,\mathrm{d}u.
\end{align}$$
This looks like it has a shot at being evaluated in terms of series expansions, but I'm stuck at this point.
Suggestions welcome.

Update:
Using the series expansions $\ln{(1-e^{-x})}=-\sum_{n=1}^{\infty}\frac{1}{n}e^{-nx}$ and $\sum_{k=1}^{\infty}e^{-kx}=\frac{1}{e^{x}-1}$, the integral becomes:
$$\begin{align}
I&=\int_{0}^{\infty}\frac{x^3}{1-e^{x}}\ln{(1-e^{-x})}\,\mathrm{d}x\\
&=\int_{0}^{\infty}\frac{x^3}{e^{x}-1}\sum_{n=1}^{\infty}\frac{1}{n}e^{-nx}\,\mathrm{d}x\\
&=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{\infty}\frac{x^3\,e^{-nx}}{e^{x}-1}\,\mathrm{d}x\\
&=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{\infty}x^3\,e^{-nx}\sum_{k=1}^{\infty}e^{-kx}\,\mathrm{d}x\\
&=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=1}^{\infty}\int_{0}^{\infty}x^3\,e^{-(k+n)x}\,\mathrm{d}x\\
&=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=1}^{\infty}\frac{6}{(k+n)^4}\\
&=\sum_{n=1}^{\infty}\frac{1}{n}\psi^{(3)}(n+1),
\end{align}$$
where $\psi^{(n)}(x)$ is the $n$-th derivative of the digamma function. Unfortunately, I don't know how to evaluate this last series. But this feels like progress.
 A: If you're interested in evaluating that $\psi'''(n+1)$ series you obtained in your original evaluation, a fun way to go about it is to use contours. A method popularized by Flajolet and Salvy in their famous paper on the topic.
$$\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$
Begin with the kernel $$f(z)=\pi \cot(\pi z)\psi'''(-z)$$
This has poles at the positive integers, n;  negative integers, -n, and 0.
The series at the positive integers is:
$$\frac{6}{(z-n)^{5}}-\frac{2\pi^{2}}{(z-n)^{3}}-\frac{\psi'''(n+1)}{z-n}+\cdot\cdot\cdot $$
Thus, the sum of the residues is:
$$Res(f,n)=\frac{6}{4!}\frac{d^{4}}{dz^{4}}\left[\frac{1}{(z-n)^{5}}\cdot \frac{1}{n}\right]-\frac{2\pi^{2}}{2!}\frac{d^{2}}{dz^{2}}\left[\frac{1}{(z-n)^{3}}\cdot \frac{1}{n}\right]-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$
$$=6\sum_{n=1}^{\infty}\frac{1}{n^{5}}-2\pi^{2}\sum_{n=1}^{\infty}\frac{1}{n^{3}}-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$
$$=6\zeta(5)-2\pi^{2}\zeta(3)-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}$$
the series at the negative integers has a simple pole.
$$Res(f, -n)=\sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}$$$
The residue at z=0 can be found by the Laurent expansion and noting the coefficient of the 1/z term:
$$Res(f,0)=24\zeta(5)$$
Now, put it altogether and set equal to 0 due to there being no poles in the contour:
$$6\zeta(5)-2\pi^{2}\zeta(3)+24\zeta(5)+\sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}-\sum_{n=1}^{\infty}\frac{\psi'''(n+1)}{n}=0$$
solve for the series in question. Note, since $\frac{\psi'''(n)-\psi'''(n+1)}{n}=\frac{6}{n^{5}}$, this means that $\displaystyle \sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}$ is equal to the series in question plus $6\zeta(5)$.  In other words by calling our series S:
$\displaystyle \sum_{n=1}^{\infty}\frac{\psi'''(n)}{n}=S+6\zeta(5)$
$$-2S+24\zeta(5)-2\pi^{2}\zeta(3)=0$$
$$S=12\zeta(5)-\pi^{2}\zeta(3)$$
A: We could use the series expansion $$\frac{\ln(1-x)}{1-x}=-\sum_{n=1}^\infty H_nx^n,$$ where $H_n=1+\frac12+\frac13+\ldots+\frac1n$ is the harmonic number. Then we get $$I=\int_0^1\frac{\ln(1-x)\ln^3x}{1-x}dx=-\sum_{n=1}^\infty H_n\int_0^1 x^n\ln^3 x\,dx=6\sum_{n=1}^\infty \frac{H_n}{(n+1)^4}.$$
The last series could be evaluated using Euler's formula $$\sum_{n=1}^\infty \frac{H_n}{n^4}=3\zeta(5)-\zeta(2)\zeta(3)$$ (see formula (20) here: http://mathworld.wolfram.com/HarmonicNumber.html): $$I=6\sum_{n=1}^\infty \frac{H_n}{(n+1)^4}=6\sum_{n=1}^\infty \left(\frac{H_{n+1}}{(n+1)^4}-\frac{1}{(n+1)^5}\right)=6(2\zeta(5)-\zeta(2)\zeta(3)).$$
A: Introduce a two-parameter deformation
$$\mathcal{I}(a,b)=\int_0^1 x^{a-1}\left[(1-x)^b-1\right] dx=\frac{\Gamma(a)\Gamma(b+1)}{\Gamma(a+b+1)}-\frac1a.$$
The integral we are looking for is obtained as
\begin{align*}
I&=\frac{\partial^3}{\partial b^3}\left[\frac{\partial\mathcal{I}(a,b)}{\partial a}\biggl|_{a=0}\right]_{b=0}=\\
&=\frac{\partial^3}{\partial b^3}\left[\frac{\gamma^2}{2}+\frac{\pi^2}{12}+\gamma\,\psi(1+b)+\frac{\psi^2(1+b)-\psi'(1+b)}{2}\right]_{b=0}=\\
&=12\zeta(5)-6\zeta(2)\zeta(3).
\end{align*}
At the last step, one needs to use that $\psi^{(n)}(1)=(-1)^{n+1}n!\,\zeta(n+1)$.
