I've conjectured the following identity for $n\geq0$ integers: $$ \int_0^1 \ln(x)^n \operatorname{Ei}(x) \, dx = (-1)^{n+1}n! \cdot \left(-\operatorname{Ei}(1)+\sum_{k=1}^{n+1} {_kF_k}\left(\begin{array}c1,1,\dots,1\\2,2,\dots,2\end{array}\middle|\,1\right)\right), $$ where $\ln$ is the natural logarithm, $\operatorname{Ei}$ is the exponential integral and ${_pF_q}$ is the generalized hypergeometric function.
For $n=0$ it is $1-e+\operatorname{Ei}(1)$ and for $n=1$ it is $e-1-\gamma$, where $e$ is base of the natural logarithm, and $\gamma$ is the Euler-Mascheroni constant.
How could we prove this identity?