Does $\int_0^{\infty} \left( p + q W \left( r e^{- s x + t} \right) + u x \right) e^{- x} d x$ have a closed-form expression?

Does $\int_0^{\infty} \left( p + q W \left( r e^{- s x + t} \right) + u x \right) e^{- x} d x$ (with 6 variables) where W is the Lambert W function (also known as ProductLog in Mathematica) have a closed-form expression? If we drop the variable $s$ from the expression Maple is able calculate $$ \int_0^{\infty} \left( p + q W \left( r e^{- x + t} \right) + u x \right) e^{- x} d x = q W \left( r e^t \right) + \frac{q}{W \left( r e^t \right)} - q + u + p - \frac{q}{r e^t} $$ which agrees with numerical calculations, so my gut feeling is that such an expression should exist for the 6-variable expression as well, but Maple nor Mathematica are able to compute it.

To simplify the problem a bit let's consider a simpler integral, $$\int_0^{\infty} W \left( e^{- s x} \right) e^{- x} d x$$ If $s=1$ then we have $$\begin{array}{ll} \int_0^{\infty} W \left( e^{- x} \right) e^{- x} d x & = \frac{1 - 2 W \left( 1 \right) + W \left( 1 \right)^2}{W \left( 1 \right)}\\ & = 0.330366124761680583225170439162 \end{array}$$ which is the solution to $$\left\{ y : y - \frac{1}{W \left( 1 \right)} = W \left( 1 \right) - 2 \right\}$$

If $s=1/2$ then Maple is not able to immediately find an answer but it is seen to be true numerically that we have(thanks to the amazing and wonderful RIES program) $$\begin{array}{ll} \int_0^{\infty} W \left( e^{- \frac{x}{2}} \right) e^{- x} d x & = \frac{- 1 + 2 W \left( 1 \right) - W \left( 1 \right)^2 + 4 W \left( 1 \right)^3}{4 W \left( 1 \right)^2}\\ & = 0.421516016690748181742333199330 \end{array}$$ which is the solution to $$\left\{ y : \cos \left( \pi \sqrt{W \left( 1 \right) - y} \right) = \sin \left( \frac{\pi}{2 W \left( 1 \right)} \right) \right\}$$

Any ideas would be greatly appreciated!

$$\int_0^{\infty} W \left( e^{- s x} \right) e^{- x} d x = W(1)+\left(-\frac{1}{s}\right)^{-\frac{1}{s}}\left(\Gamma\left(\frac{1}{s}\right)-s\;\Gamma\left(1+\frac{1}{s},-\frac{W(1)}{s}\right)\right),$$

where $\Gamma(z)$ is the Gamma function and $\Gamma(a,z)$ is the incomplete Gamma function.

To calculate this, change variable $y=e^{-x}$, then transpose the function graph to replace $W(y)$ with its inverse function $e^z z$.