Here is another solution. We start with @Omran Kouba's simplification:
$$ I = -\frac{7}{24}J + \frac{1}{3}K, $$
where
$$ J := \int_{0}^{1} \frac{\log(1-x)\log^3 x}{1-x} \, \mathrm{d}x
\qquad\text{and}\qquad
K := \int_{0}^{1} \frac{\log^3 x}{1-x^2} \log\left(\frac{1-x}{1+x}\right) \, \mathrm{d}x. $$
Step 1. To make use of this decomposition, we first introduce two auxiliary integrals:
$$ A := \int_{0}^{1} \frac{\log (1-x)\log^3 x}{(1-x)x} \, \mathrm{d}x
\qquad\text{and}\qquad
B := \int_{0}^{1} \frac{\log (1+x)\log^3 x}{(1+x)x} \, \mathrm{d}x.
$$
Plugging the power series $\frac{\log(1-x)}{1-x} = -\sum_{n=1}^{\infty} H_n x^n$ into the integral and performing term-wise integration, we easily find that
$$ \color{red}{ A = 6 \sum_{n=1}^{\infty} \frac{H_n}{n^4} }
\qquad\text{and}\qquad
\color{red}{ B = 6 \sum_{n=1}^{\infty} (-1)^n \frac{H_n}{n^4} }. $$
We will demonstrate how to compute these sums in Step 3.
Step 2. Next, we relate $A$ and $B$ to $J$ and $K$. We achieve this by decomposing the respective integrand of $A$ and $B$ algebraically. Indeed, for $J$, we can write
$$\begin{align*}
\color{red}{J}
&= \int_{0}^{1} \left[ \color{blue}{ \frac{\log(1-x)}{1-x} } \right] \log^3 x \, \mathrm{d}x \\
&= \int_{0}^{1} \left[ \color{blue}{ \frac{\log(1-x)}{(1-x)x} - \frac{\log(1-x)}{x} } \right] \log^3 x \, \mathrm{d}x \\
&= \color{red}{ A - 6\zeta(5)}.
\end{align*}$$
The blue part in the first and second line are equal algebraically. Similarly, for $K$, we have
$$\begin{align*}
\color{red}{K}
&= \int_{0}^{1} \left[ \color{blue}{ \frac{1}{1-x^2} \log\left(\frac{1-x}{1+x}\right) } \right] \log^3 x \, \mathrm{d}x \\
&= \int_{0}^{1} \left[ \color{blue}{ \frac{\log (1-x)}{(1-x)x} + \frac{\log (1+x)}{(1+x)x} - \frac{\log (1-x^2)}{(1-x^2)x} } \right] \log^3 x \, \mathrm{d}x \\
&= \int_{0}^{1} \frac{\log (1-x)\log^3 x}{(1-x)x} \, \mathrm{d}x
+ \int_{0}^{1} \frac{\log (1+x)\log^3 x}{(1+x)x} \, \mathrm{d}x
- \smash[b]{ \underbrace{ \int_{0}^{1} \frac{\log (1-x^2)\log^3 x}{(1-x^2) x} \, \mathrm{d}x }_{x^2 \mapsto x} } \\
&= \color{red}{\frac{15}{16} A + B}.
\end{align*}$$
Again, the blue part in the first and second line are identical as functions. Therefore we can express $I$ as a linear combination of $A$, $B$, and $\zeta(5)$:
$$ \color{red}{ I = \frac{1}{48} A + \frac{1}{3} B + \frac{7}{4}\zeta(5). } $$
Step 3. So it remains to compute $A$ and $B$. To this end, we extend the harmonic numbers $H_n$ to the meromorphic function on $\mathbb{C}$ via the identity $H_t = \psi(t+1) - \psi(1)$, where $\psi$ is the digamma function. Then by the Euler's reflection formula, we have
$$ H_t = H_{-1-t} - \frac{\pi}{\tan(\pi t)}. \tag{*} $$
This relation will be used critically in our computation.
Now fix a number $0 < \sigma < 1$. Then by the residue calculation followed by the identity $\text{(*)}$, we obtain
$$\begin{align*}
A
&= \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{H_t}{t^4} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t \\
&= \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{H_{-1-t}}{t^4} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t
- \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{1}{t^4} \frac{\pi^2}{\tan^2(\pi t)} \, \mathrm{d}t
\end{align*}$$
To proceed, we apply the substitution $t \mapsto -t$ to the first term. This leads to:
$$\begin{align*}
&\frac{6}{2\pi i}\int_{\sigma+i\infty}^{\sigma-i\infty} \frac{H_{-1-t}}{t^4} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t \\
&= - \frac{6}{2\pi i}\int_{-\sigma+i\infty}^{-\sigma-i\infty} \frac{H_{t-1}}{t^4} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t \\
&= -6 \, \underset{t=0}{\mathrm{Res}} \, \left[ \frac{H_{t-1}}{t^4} \frac{\pi}{\tan(\pi t)} \right] - \frac{6}{2\pi i}\int_{\sigma+i\infty}^{\sigma-i\infty} \frac{H_t - t^{-1}}{t^4} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t \\
&= -6 \left[ 2\zeta(2)\zeta(3) - \zeta(5) \right] - A + 6\zeta(5).
\end{align*}$$
For the second term, note that we have $\int_{\sigma+i\infty}^{\sigma-i\infty} t^{-4} \, \mathrm{d}t = 0$. So we can add a constant to $\frac{\pi^2}{\tan^2(\pi t)}$ without changing the value of the second integral. In particular, we can use $\frac{\pi^2}{\sin^2(\pi t)} = \pi^2 + \frac{\pi^2}{\tan^2(\pi t)}$ instead. Since $(\frac{\pi}{\tan(\pi t)})' = -\frac{\pi^2}{\sin^2(\pi t)}$, we can then perform integration by parts to obtain
$$\begin{align*}
- \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{1}{t^4} \frac{\pi^2}{\tan^2(\pi t)} \, \mathrm{d}t
&= - \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{1}{t^4} \frac{\pi^2}{\sin^2(\pi t)} \, \mathrm{d}t \\
&= \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{4}{t^5} \frac{\pi}{\tan(\pi t)} \, \mathrm{d}t \\
&= 24\zeta(5).
\end{align*}$$
Combining altogether and solving the resulting equation for $A$, we obtain
$$ \color{red}{ A = -6\zeta(2)\zeta(3) + 18\zeta(5) } $$
A similar technique shows that we also have
$$\begin{align*}
\color{red}{B}
&= \frac{6}{2\pi i} \int_{\sigma+i\infty}^{\sigma-i\infty} \frac{H_t}{t^4} \frac{\pi}{\sin(\pi t)} \, \mathrm{d}t
= \color{red}{ 3\zeta(2)\zeta(3) - \frac{177}{16}\zeta(5) }.
\end{align*}$$
Plugging these back into $I$, the desired identity is established.