Define integral for $\gamma,\zeta(i) i\in\mathbb{N}$ and Stirling numbers of the first kind Consider the integral
$$\int\limits_0^{\infty}e^{-x}x^k\ln(x)^n\dfrac{dx}x$$
For $n=3$ we have
$$(-\gamma^2-2\zeta(3)-3\zeta(2)\gamma)\genfrac{[}{]}{0pt}{}{k}{1}+3(\gamma^2+\zeta(2))\genfrac{[}{]}{0pt}{}{k}{2}-6\gamma\genfrac{[}{]}{0pt}{}{k}{3}+6\genfrac{[}{]}{0pt}{}{k}{4}$$
where $\genfrac{[}{]}{0pt}{}{n}{k}$ is Stirling number of the first kind. For different n there are simmilar formulas. In the general case
$$\int\limits_0^{\infty}e^{-x}x^k\ln(x)^n\dfrac{dx}x\in\bigoplus_{j=1}^{n+1}\genfrac{[}{]}{0pt}{}{k}{j}\mathbb{Z}[\gamma,\zeta(i)]_{n+1-j}$$
where $$\mathbb{Z}[\gamma,\zeta(i)]=\bigoplus_{j=0}^{\infty}\mathbb{Z}[\gamma,\zeta(i)]_{j}$$ is graded ring.
Is ir possible to generalize this formula to multiple zeta values?
 A: This is not a complete answer, rather two observations:
Let $\mathcal{I}_{k,n}$ denote the integral. Integration by parts, assuming $k \geqslant  2$ and $n \geqslant 1$ gives the recurrence equation:
$$
    \mathcal{I}_{k+1,n} = n  \mathcal{I}_{k,n-1} + k \mathcal{I}_{k,n}  \tag{1}
$$
which we can use to reduce the problem to $k=1$. Let $\mathcal{J}_{k,n} = \frac{1}{n!}\mathcal{I}_{k,n}$. It is readily seen that $\mathcal{J}_{k,n}$ satisfies the recurrence relation of the unsigned Stirling numbers of the first kind:
$$
  \mathcal{J}_{k+1,n} =  \mathcal{J}_{k,n-1} + k \mathcal{J}_{k,n}  \tag{2}
$$
One the other hand:
$$
   \frac{\mathrm{d}^n}{\mathrm{d} s^n} x^{s-1} = x^{s-1} \log(x)^n
$$
Thus:
$$
   \int_0^\infty \mathrm{e}^{-x} x^{k-1} \log(x)^n \mathrm{d} x = \left. \frac{\mathrm{d}^n}{\mathrm{d} s^n} x^{s-1} \int_0^\infty \mathrm{e}^{-x} x^{s-1}  \mathrm{d} x \right|_{s=k} = \left. \frac{\mathrm{d}^n}{\mathrm{d} s^n} \Gamma(s) \right|_{s=k} = \frac{\mathrm{d}^n}{\mathrm{d} k^n} \exp(\log\Gamma(k) )
$$
Now use Faà di Bruno formula:
$$
  \frac{\mathrm{d}^n}{\mathrm{d} k^n} \exp(\log\Gamma(k) ) = \Gamma(k) \sum_{m=0}^n B_{n,m}\left( \psi(k), \psi^{(1)}(k), \ldots, \psi^{(n-1)}(k)\right)
$$
where $\psi(x)$ is the digamma function, and $\psi^{(m)}(x)$ is the polygamma function, and where $B_{n,k}$ stand for the Faà di Bruno polynomials, which are homogeneous:
$$
   B_{n,k}\left(\lambda x_1, \lambda x_2, \ldots, \lambda x_n\right) =  \lambda^k B_{n,k}\left(x_1, x_2, \ldots, x_n\right) 
$$
$$
   B_{n,k}\left(\lambda x_1, \lambda^2 x_2, \ldots, \lambda^n x_n\right) =  \lambda^n B_{n,k}\left(x_1, x_2, \ldots, x_n\right) 
$$
Assume that recurrence relation $(1)$ was used to reduce to $k=1$ case. For this $k$, values of polygammas are related to $\zeta$-values:
$$
   \psi(1) = -\gamma \qquad \psi^{(m)}(1) = \left. \frac{\mathrm{d}^m}{\mathrm{d} s^m} \psi(s) \right|_{s=1} = (-1)^{m+1} m! \sum_{n=1}^\infty \frac{1}{n^{m+1}} = (-1)^{m+1} m! \zeta(m+1)
$$
Thus
$$
  \int_0^\infty \mathrm{e}^{-s} \log(x)^n \mathrm{d} x = (-1)^n \sum_{m=0}^n B_{n,m} \left(\gamma, \zeta(2), 2! \zeta(3), \ldots, n! \zeta(n+1)\right)
$$
and this explains the grading.
